2.1 Obtention des données

La première étape consiste à créer de la requête “overpass”, qui est l’API couramment utilisée pour télécharger des données du serveur OpenStreetMap.

gp_osmdata_l8 <- opq("Guadeloupe") %>%
  add_osm_feature(key="admin_level", value="8") %>%
  osmdata_sf()

Ici, on précise admin_level=8 qui correspond aux données du niveau communal, et on demande à enregistrer le résultat de la requête dans un objet (gp_osmdata_sf) qui est une liste de dataframes au format sf (plus d’informations sur ce format). A noter que la requête overpass peut nécessiter plusieurs essais avant d’aboutir, si le serveur ne répond pas (cela dépend de la charge du serveur).

Pour enregistrer les données sur le disque, et ne pas avoir à les télécharger à chaque nouvelle session, il suffit d’entrer la commande write_rds(gp_osmdata_l8, "data/gp_osmdata_l8.rds")

write_rds(gp_osmdata_l8, "gp_osmdata_l8.rds", compress="xz")

Et pour recharger les données à partir de ce fichier :

gp_osmdata_l8 <- read_rds("gp_osmdata_l8.rds")

Regardons de quoi est constitué notre premier objet :

gp_osmdata_l8

## Object of class 'osmdata' with:
##                  $bbox : 15.8320085,-61.809764,16.5144801,-61.0003663
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 56852 points
##             $osm_lines : 'sf' Simple Features Collection with 375 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 59 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : 'sf' Simple Features Collection with 32 multipolygons

Nous voyons que que cette liste contient un data-frame de classe sf appelé osm_multipolygons formé de 32 multipolygones. C’est ce data-frame qui contient les données correpondant aux 32 communes que compte la Guadeloupe.

Regardons sa structure en listant les noms des variables qu’il contient, telles que le nombre de population, le code postal, le code insee, le nom (de la commune), etc…

names(gp_osmdata_l8$osm_multipolygons)

##  [1] "osm_id"             "name"               "addr.postcode"     
##  [4] "admin_level"        "alt_name"           "alt_name.fr"       
##  [7] "alt_name.gcf"       "boundary"           "capital"           
## [10] "label.FR.history"   "loc_name"           "loc_name.fr"       
## [13] "loc_name.gcf"       "name.be"            "name.be.tarask"    
## [16] "name.bg"            "name.el"            "name.fr"           
## [19] "name.frp"           "name.gcf"           "name.he"           
## [22] "name.ja"            "name.ka"            "name.ko"           
## [25] "name.lt"            "name.mr"            "name.ru"           
## [28] "name.sr"            "name.uk"            "name.zh"           
## [31] "population"         "population.date"    "postal_code"       
## [34] "ref.FR.SIREN"       "ref.INSEE"          "ref.LOCODE"        
## [37] "source"             "source.population"  "source.postal_code"
## [40] "type"               "wikidata"           "wikipedia"         
## [43] "geometry"

La dernière variable geometry contient les données spatiales (positions spatiales des multipolygones) et les autre contiennent les métadonnées rattachées à ces mutltipolygones.

Nous pouvons sélectionner les variables pour ne garder que celles qui nous semblent intéressantes :

gp_com_lim_sf <- gp_osmdata_l8$osm_multipolygons %>%
  dplyr::select(name, ref.INSEE, population) %>%
  arrange(name) %>%
  st_make_valid()

On obtient un data-frame de classe sf avec les 4 variables d’intérêt (nom de la commune, numéro INSEE, taille de la population, géométrie), qui donne les limites des 32 communes de la Guadeloupe.

L’avantage d’avoir un data-frame de classe sf est que l’on peut effectuer des calculs et des transformations sur les variables qu’il contient. Par exemple, on pourrait vouloir calculer la densité de population au km² : on dispose en effet de la taille de la population de chaque commune, et à partir des géométries on peut calculer les surfaces des multipolygones.

Il faut commencer par convertir la variable population en classe numeric car pour l’instant elle est stockée sous forme character puis on utilise la fonction st_area pour calculer la superficie des communes et la transformer en km² (car l’unité de base est le mètre) :

# conversion de la classe de la variable population
gp_com_lim_sf$population <- as.numeric(gp_com_lim_sf$population)

# calcul de la superficie et conversion en km²
gp_com_lim_sf$superficie <- as.numeric(units::set_units(st_area(gp_com_lim_sf), km^2))

Puis nous pouvons calculer les densités de population et ajouter les valeurs à une nouvelle variable population_density :

gp_com_lim_sf$population_density <- gp_com_lim_sf$population/gp_com_lim_sf$superficie

Pour l’instant, on dispose donc d’un objet contenant des mutlipolygones avec le contour des communes et quelques variables d’intérêt. Si par exemple on veut afficher la localisation des centres des communes (bourgs) il faut chercher ces données dans un autre data-frame.

En fait ces données sont présentes dans l’objet initialement créé gp_osmdata_l8, mais comme il s’agit de localisation précises, ponctuelles, il ne faut plus les chercher dans osm_multipolygons mais dans osm_points. Or, ce dta-frame contient de mutliples points d’intérêts en plus des centre-bourgs. Explorons le data-frame gp_osmdata_l8$osm_points :

names(gp_osmdata_l8$osm_points)

##  [1] "osm_id"                   "name"                    
##  [3] "addr.postcode"            "admin_level"             
##  [5] "alt_name"                 "alt_name.fr"             
##  [7] "alt_name.gcf"             "amenity"                 
##  [9] "bus"                      "capital"                 
## [11] "capital_ISO3166.1"        "cargo"                   
## [13] "crossing"                 "description"             
## [15] "direction"                "ele"                     
## [17] "ferry"                    "ford"                    
## [19] "height"                   "highway"                 
## [21] "intermittent"             "junction"                
## [23] "loc_name"                 "loc_name.fr"             
## [25] "loc_name.gcf"             "name.be"                 
## [27] "name.be.tarask"           "name.bg"                 
## [29] "name.el"                  "name.fr"                 
## [31] "name.frp"                 "name.gcf"                
## [33] "name.he"                  "name.ja"                 
## [35] "name.ka"                  "name.ko"                 
## [37] "name.lt"                  "name.mr"                 
## [39] "name.ru"                  "name.sr"                 
## [41] "name.uk"                  "name.zh"                 
## [43] "natural"                  "place"                   
## [45] "population"               "public_transport"        
## [47] "ref.FR.SIREN"             "ref.LOCODE"              
## [49] "seamark.harbour.category" "source"                  
## [51] "tourism"                  "traffic_calming"         
## [53] "waterway"                 "waterway.sign"           
## [55] "wikidata"                 "wikipedia"               
## [57] "geometry"

On trouve une variable nommée addr.postcode dont on peut afficher le contenu :

gp_osmdata_l8$osm_points$addr.postcode

##   [1] NA      NA      NA      NA      NA      NA      NA      NA      NA     
##  [10] NA      NA      NA      NA      NA      NA      NA      NA      NA     
##  [19] NA      NA      NA      NA      NA      NA      NA      NA      NA     
##  [28] NA      NA      NA      NA      NA      NA      NA      NA      "97100"
##  [37] "97110" NA      "97160" NA      NA      NA      NA      NA      NA     
##  [46] NA      NA      NA      NA      NA      "97114" "97113" NA      NA     
##  [55] NA      NA      NA      NA      NA      NA      NA      NA      NA     
##  [64] NA      NA      NA      NA      NA      NA      NA      NA      NA     
##  [73] NA      NA      NA      NA      NA      NA      NA      NA      NA     
##  [82] NA      NA      NA      NA      NA      NA      NA      "97121" NA     
##  [91] NA      NA      NA      NA      NA      NA      NA      NA      NA     
## [100] NA     
##  [ reached getOption("max.print") -- omitted 56752 entries ]

On constate qu’il y a beaucoup de NA et seulement quelques valeurs réelles correspondant à un code postal, donc probablement rattachées aux centre-bourgs. Si tel est le cas, vérifions s’il y en a bien 32 :

length(which(!is.na(gp_osmdata_l8$osm_points$addr.postcode)))

## [1] 32

Parfait, nous avons bien 32 points auxquels sont rattachés un code postal. Construisons un nouvel objet sf avec ces points et 2 variables d’intérêt, name (noms des communes) et addr.postcode (code postaux) :

gp_com_cen_sf <- gp_osmdata_l8$osm_points %>%
  dplyr::filter(!is.na(addr.postcode)) %>%
  dplyr::select(name, addr.postcode) %>%
  arrange(name) %>%
  st_make_valid()

2.2 Visualisation des cartes

Nous pouvons faire nos premières cartes statiques, en utilisant ggplot2 qui permet d’afficher les données au format sf grâce à la commande geom_sf() :

p1 <- ggplot() +
  geom_sf(data=gp_com_lim_sf) +
  labs(title = "Contours des communes de la Guadeloupe", x = "Longitude", y = "Latitude") +
  annotation_scale(style="ticks", location="br") +
  theme_bw()
p1

Nous pouvons aussi choisir de colorer chaque commune avec une couleur arbitraire en ajoutant l’esthétique adéquate (aes(fill=name) et un échelle qualitative de couleurs (scale_fill_discrete):

p2 <- p1 +
  geom_sf(data = gp_com_lim_sf, aes(fill = name)) +
  scale_fill_discrete() +
  guides(fill="none")
p2

Vous aurez remarqué que nous sommes repartis du premier graphique (objet p1) ce qui permet de ne changer que le strict nécessaire (ce qui a été défini dans l’objet p1 et qui ne change pas, tels que le thème, les labels des axes et le titre, n’ont pas besoin d’être spécifiés à nouveau)

A cette carte, nous pouvons ajouter des points chaque centre-bourg et des labels avec les noms des communes correspondantes (on utilisera ici la fonction geom_text_repel du paquet ggrepel du paquet du même nom, qui permet d’éviter aux labels de se superposer) :

p3 <- p2 +
  geom_sf(data = gp_com_cen_sf, size = 2, shape=21, color="black", fill="white") +
  geom_text_repel(data=gp_com_cen_sf,
                            aes(label = name, geometry = geometry), fontface="bold",
                            stat = "sf_coordinates", size=3, min.segment.length = 0) +
  labs(title = "Communes de la Guadeloupe")
p3

La pleine puissance du format sf associé à ggplot2 s’exprime quand on veut des affichages condtionnels. Par exmple, on peut faire en sorte que les couleurs des multipolygones des communes dépende de la taille de la population :

p4 <- p1 +
  geom_sf(data = gp_com_lim_sf, aes(fill = population)) +
  scale_fill_distiller(palette="OrRd", direction=1) +
  labs(title = "Population par commune")
p4

Ou encore, en fonction de la densité de population :

p5 <- p1 +
  geom_sf(data = gp_com_lim_sf, aes(fill = population_density)) +
  scale_fill_distiller(palette="OrRd", direction=1) +
  labs(title = "Densité de population par commune")
p5

2.3 Cartes de type raster

Précédemment, nous avons traité des objets spatiaux de classe sf, qui peuvent être des polygones, des lignes ou points, dont on connait les coordonnées des contituants dans l’espace. Il existe d’autres types de données spatiales qui correspondent schématiquement à une matrice de points positionnés sur une grille régulière, et de coordonnées (x, y). A chaque point peut être associée une ou plusieurs valeurs. Ce type d’objets correspond à la grande classe des rasters.

Dans le domaine de l’information géographique, typiquement ce type d’objets sont obtenus par des relevés satellites, qui mesurent certaines variables à des points régulièrement espacés sur une grille : données météorologiques, données radiométriques, données d’élévation (altitude) du terrain…

Ici, nous allons utiliser ce type données pour montrer comment les intégrer à une carte sou R. Avec le paquet elevatr nous pouvons télécharger les données d’élévation sur toutes les régions du du globe, à différents niveaux de précision. La commande get_elev_raster télécharge les données à partir de la source de notre choix, et peut utiliser un objet spatial obtenu au préalable pour restreindre les données à la même zone.

Par exemple, utilisons les données OpenStreetMap au niveau admin_level=7 pour extraire les multipolygones correspondant aux 2 arrondissements de la Guadeloupe, Pointe-àPitre et Basse-Terre, ce qui donne à peu près le découpage de la Grande-Terre et de la Basse-Terre, avec les dépendances :

gp_osmdata_l7 <- opq("Guadeloupe") %>%
  add_osm_feature(key="admin_level", value="7") %>% # 
  osmdata_sf()

write_rds(gp_osmdata_l7, "gp_osmdata_l7.rds", compress = "xz")

gp_osmdata_l7 <- read_rds("gp_osmdata_l7.rds")

gp_gtbt_lim_sf <- gp_osmdata_l7$osm_multipolygons %>%
  dplyr::select(name) %>%
  arrange(name) %>%
  st_make_valid()

Verifions à quoi cela ressemble, avec la fonction plot de base :

plot(gp_gtbt_lim_sf, key.width=lcm(4))

Utilisons cet objet gp_gtbt_lim_sf comme “masque” pour le téléchargement et l’enregistrement des données d’élévation :

gp_elev_raster <- get_elev_raster(locations = gp_gtbt_lim_sf,
                                  z = 12,
                                  clip = "locations",
                                  src="aws"
                                  ) %>%
  rast()

# ramener les valeurs négatives à zéro
gp_elev_raster <- app(gp_elev_raster, fun=function(x){ x[x < 0] <- 0; return(x)} )

A noter qu’ici nous utilisons la fonction terra::rast() au lieu de la fonction plus classique raster::raster car la première est plus récente, rapide et efficiente (le paquet terra est voué à remplacer raster).

Une visualisation rapide :

plot(gp_elev_raster, col=terrain.colors(n=40))

Pour le fun, on peut aussi créer une visualisation de l’ombrage du relief à partir des données d’élévation :

gp_ombr_raster <- shade(terrain(gp_elev_raster, v='slope')/180*pi,
                        terrain(gp_elev_raster, v='aspect')/180*pi,
                        angle=40, direction=330
                        )

plot(gp_ombr_raster, col=grey(0:100/100))

Comme les données raster sont assez lourdes, pour les visualisations plus complexes (surtout avec leaflet) il vaut mieux réduire leur résolution. On utilise ici la fonction terra::aggregate avec fact=2 pour diviser la résolution linéaire par 2 :

gp_elev_raster_2 <- aggregate(gp_elev_raster, fact=2)
gp_ombr_raster_2 <- aggregate(gp_ombr_raster, fact=2)

Avec ggplot2 nous poivons superposer la carte des communes sur le fond du relief :

ggplot() +
  layer_spatial(gp_elev_raster_2) +
  scale_fill_gradientn(colors=terrain.colors(n=40), na.value = 0) +
  geom_sf(data=gp_com_lim_sf, fill=NA) +
  theme_bw()

3.1 Avec plotly

L’intérêt de plotly est qu’il procure la fonction ggplotly, qui convertit simplement un objet ggplot2::ggplot() en un object plotly interactif. Or, dans la première partie nous avons créé de tels objets, que nous pouvons donc convertir directement :

dynp1 <- ggplotly(p1)
dynp1

dynp4 <- ggplotly(p4)
dynp4

On peut interagir avec les graphiques grâce aux boutons situés en haut à droite.

Remarquez qu’en passant le curseur de la souris au dessus d’une zone, une étiquette apparaît automatiquement avec la valeur correspondant à la zone la plus proche du curseur. Cependant, si on pointe au centre d’un polygone aucune valeur ne s’affiche, il faut placer le curseur près du contour, et dans certains cas on a du mal à savoir à quelle zone une valeur correpond quand 2 zones sont contiguës.

3.2 Avec leaflet

Avec leaflet l’approche est différente. Il recréer de nouveaux graphique, mais la syntaxe est assez simple et elle utilise les objets de classe sf déjà créés :

leaflet() %>%
  addPolygons(data=gp_com_lim_sf,
              color="black", opacity=1, weight=2,
              fillOpacity = 0,
              )

De prime abord, cette carte est tout simple, elle n’affiche pas les valeurs par défaut, par contre on constate que la navigation est très simple (déplacement, et zoom avec la molette de la souris).

Mais leaflet permet de configurer assez finement la visualisation en ajoutant les commandes adéquates. Et autre avantage, leaflet permet d’afficher plusieurs couches et de sélectionner les couches à visualiser. La syntaxe devient plus compliquée, mais confère à leaflet toute sa puissance.

Ci-dessous, voici une carte beaucoup plus complexe, mais qui permet de visualiser beaucoup d’informations, avec en fond de carte les tuiles OpenStreetMap :

leaflet() %>%
  addTiles(group="Fond OSM (default)") %>%
  addRasterImage(raster::raster(gp_ombr_raster_2), group="Fond relief",
                 colors=grey.colors(100)) %>%
  addRasterImage(raster::raster(gp_elev_raster_2), group="Fond élévation",
                 colors=terrain.colors(n=20, alpha=0.5, rev = F)) %>%
  addPolygons(data=gp_com_lim_sf, group="Limites communales",
             color="black", opacity=0.5, weight=2,
                fillOpacity = 0,
            label = ~population
              ) %>%
  addPolygons(data=gp_com_lim_sf, group="Population density",
              color="black", opacity=0, weight=2,
             fillColor=~colorNumeric(palette = "YlOrRd",
                                      domain = population_density)(population_density),
             fillOpacity = 1,
             label = ~ population_density
              ) %>%
  addCircles(data=gp_com_cen_sf, group="Centres villes",
            radius=500, opacity=0.8, weight=2, color="black",
            fill=T, fillColor="white", fillOpacity = 0.5,
            label = ~as.character(name)
            ) %>%
  addLegend(pal=colorNumeric(palette = "YlOrRd", domain = gp_com_cen_sf$population_density),
            values=gp_com_lim_sf$population_density,
            opacity=0.9, title = "Densité de population (au km2)",
            position = "bottomleft", group="Population density"
           ) %>%
  addScaleBar() %>%
  addLayersControl(
    baseGroups = c("Fond OSM (default)"),
    overlayGroups = c("Fond relief",
                      "Fond élévation",
                      "Population density",
                      "Limites communales",
                      "Centres villes"
                      ),
  options = layersControlOptions(collapsed = T)
)