Réaliser une carte avec R (6)

Suite aux parties 1, 2, 3, 4, et 5, nous allons à présent voir quelques transformations et opérations sur les données géospatiales.

Vous pouvez télécharger le code source complet du notebook R de ce tutoriel ici et le fichier contenant uniquement le code R ici.

Projection

Tout d’abord réalisons une carte du globe à partir des données vectorielles de Natural Earth. À cette échelle nous n’avons pas besoin de données très précises. Nous allons utiliser le thème theme_map() de ggthemes pour plus de clarté.

require(rgdal)
require(ggplot2)
require(ggthemes)
# http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_land.zip
wmap <- readOGR(dsn = path.expand("~/Téléchargements/ne_110m_land"), layer = "ne_110m_land")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/thomas/Téléchargements/ne_110m_land", layer: "ne_110m_land"
## with 127 features
## It has 2 fields
wmapF <- fortify(wmap)
ggplot(wmapF, aes(x = long, y = lat, group=group)) + 
  geom_polygon() + 
  coord_fixed() +
  theme_map()

Rajoutons le graticule et la bordure de la carte.

# http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_graticules_all.zip
grat <- readOGR(dsn = path.expand("~/Téléchargements/ne_110m_graticules_all"), layer = "ne_110m_graticules_15")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/thomas/Téléchargements/ne_110m_graticules_all", layer: "ne_110m_graticules_15"
## with 35 features
## It has 5 fields
bord <- readOGR(dsn = path.expand("~/Téléchargements/ne_110m_graticules_all"), layer = "ne_110m_wgs84_bounding_box")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/thomas/Téléchargements/ne_110m_graticules_all", layer: "ne_110m_wgs84_bounding_box"
## with 1 features
## It has 2 fields
gratF <- fortify(grat)
bordF <- fortify(bord)
ggplot(data = wmapF, aes(x = long, y = lat, group=group)) + 
  geom_polygon(data = bordF, fill = "white", color = "black") + 
  geom_polygon() + 
  geom_path(data = gratF, color = "grey50", linetype = "dashed", size = .25) +
  coord_fixed() +
  theme_map()

Un problème de base de la cartographie est celui de la projection du globe terrestre sur la surface plane d’une carte en deux dimensions. Les fichiers shp sont en général accompagnés d’un fichier prj qui précise la projection utilisée. On peut vérifier également cette projection dans R avec la fonction proj4string(). Il convient de s’assurer que tous les éléments cartographiés utilisent bien la même projection.

proj4string(wmap)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(grat)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(bord)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(wmap) == proj4string(grat)
## [1] TRUE

Sans surprise, il est facile dans R de changer de projection. Il y a deux méthodes. La première utilise la fonction coord map() de ggplot2. Il y toutefois quelques problèmes avec la clôture de polygones dans certaines projections (notamment pour l’Antarctique), et il peut être utile de travailler avec des données transformées sans forcément les cartographier.

map <- ggplot(data = wmapF, aes(x = long, y = lat, group=group)) + 
  geom_polygon(data = bordF, fill = "white", color = "black") + 
  geom_polygon() + 
  geom_path(data = gratF, color = "grey50", linetype = "dashed", size = .25) +
  theme_map()
map + coord_map("mercator", xlim = c(-180,180), ylim = c(-90,90))

map + coord_map("sinusoidal", xlim = c(-180,180), ylim = c(-90,90))

map + coord_map("cylindrical", xlim = c(-180,180), ylim = c(-90,90))

map + coord_map("gilbert", xlim = c(-180,180), ylim = c(-90,90))

map + coord_map("mollweide", xlim = c(-180,180), ylim = c(-90,90))

map + coord_map("cylindrical", xlim = c(-180,180), ylim = c(-90,90))

map + coord_map("polyconic", xlim = c(-180,180), ylim = c(-90,90))

L’autre méthode consiste à transformer les données au préalable avec la fonction spTransform() qui prend comme arguments l’élément géospatial à transformer et la spécification de projection. Cette dernière peut être spécifiée de manière semi-automatique avec la fonction CR() et les codes EPSG des projections et systèmes de réferences officiels et courants (le système GPS a le code 4326, la plupart des cartes web utilisent l’EPSG 3857). Une liste avec exemples de projections et de leur codes est disponible ici .Dans ce cas il faut utiliser dans ggplot la fonction coord_fixed() au lieu de coord_map().

wmapT <- spTransform(wmap, CRS("+proj=robin"))
proj4string(wmapT)
## [1] "+proj=robin +ellps=WGS84"
wmapTF <- fortify(wmapT)
gratT <- spTransform(grat, CRS(proj4string(wmapT)))
gratTF <- fortify(gratT)
bordT <- spTransform(bord, CRS(proj4string(wmapT)))
bordTF <- fortify(bordT)
ggplot(data = wmapTF, aes(x = long, y = lat, group=group)) + 
  geom_polygon(data = bordTF, fill = "white", color = "black") +
  geom_polygon() +
  geom_path(data = gratTF, color = "grey50", linetype = "dashed", size = .25) +
  coord_fixed() +
  theme_map()

wmapT <- spTransform(wmap, CRS("+proj=bonne +lat_1=10"))
wmapTF <- fortify(wmapT)
gratT <- spTransform(grat, CRS(proj4string(wmapT)))
gratTF <- fortify(gratT)
bordT <- spTransform(bord, CRS(proj4string(wmapT)))
bordTF <- fortify(bordT)
ggplot(data = wmapTF, aes(x = long, y = lat, group=group)) + 
  geom_polygon(data = bordTF, fill = "white", color = "black") +
  geom_polygon() +
  geom_path(data = gratTF, color = "grey50", linetype = "dashed", size = .25) +
  coord_fixed() +
  theme_map()

wmapT <- spTransform(wmap, CRS("+proj=aeqd"))
wmapTF <- fortify(wmapT)
gratT <- spTransform(grat, CRS(proj4string(wmapT)))
gratTF <- fortify(gratT)
bordT <- spTransform(bord, CRS(proj4string(wmapT)))
bordTF <- fortify(bordT)
ggplot(data = wmapTF, aes(x = long, y = lat, group=group)) + 
  geom_polygon(data = bordTF, fill = "white", color = "black") +
  geom_polygon() +
  geom_path(data = gratTF, color = "grey50", linetype = "dashed", size = .25) +
  coord_fixed() +
  theme_map()

wmapT <- spTransform(wmap, CRS("+proj=cass"))
wmapTF <- fortify(wmapT)
gratT <- spTransform(grat, CRS(proj4string(wmapT)))
gratTF <- fortify(gratT)
bordT <- spTransform(bord, CRS(proj4string(wmapT)))
bordTF <- fortify(bordT)
ggplot(data = wmapTF, aes(x = long, y = lat, group=group)) + 
  geom_polygon(data = bordTF, fill = "white", color = "black") +
  geom_polygon() +
  geom_path(data = gratTF, color = "grey50", linetype = "dashed", size = .25) +
  coord_fixed() +
  theme_map()

Opérations sur les objets géospatiaux

Imaginons que nous voulions réaliser une carte de l’Asie en couleur mettant en lumière l’archipel des Ryūkyū. Nous avons déjà vu ([1], [2]) comment afficher une carte Stamen et y superposer des polygones en couleur. Il suffit donc de récupérer tout d’abord le fond de carte.

require(ggmap)
xlims <- c(110,150)
ylims <- c(20,50)
asia <- get_stamenmap(bbox = c(left = min(xlims), bottom = min(ylims), right = max(xlims), top = max(ylims)), maptype = "terrain-background", zoom = 5)

Ensuite nous allons trier nos polygones pour ne garder que ceux des îles Ryūkyū. Premièrement, nous allons créer le sous-ensemble des polygones étiquetés comme appartenant à la préfecture d’Okinawa, sauf quelques uns clairement identifiés et qui ne font pas partie des îles Ryūkyū. Ensuite nous allons étiqueter les polygones selon si oui ou non leur centroïde est situé au-dessus ou non d’une certaine latitude, et ne conserver que ceux appartenant à la préfecture de Kagoshima et étiquetés comme en dessous de la latitude choisie.

require(ggmap)
xlims <- c(110,150)
ylims <- c(20,50)
jpn2 <- readOGR(dsn = path.expand("~/repos/cartes/jpn/adm/jpn"), layer = "JPN_adm2")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/thomas/repos/cartes/jpn/adm/jpn", layer: "JPN_adm2"
## with 1811 features
## It has 11 fields
oki <- subset(jpn2, NAME_1 == "Okinawa" & !(NAME_2 %in% c("Minami Daito","Unknown5","Unknown8")))
jpn2@data$south <- coordinates(jpn2)[,2] < 28.5
ama <- subset(jpn2, NAME_1 == "Kagoshima" & south == TRUE)
okiF <- fortify(oki)
amaF <- fortify(ama)

Enfin il suffit de rajouter les polygones sur la carte

ggmap(asia) + 
  geom_polygon(data = okiF, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .5) +
  geom_polygon(data = amaF, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .5) +
  theme_map()

Le résultat n’est pas des plus clairs, et on aimerait tracer une zone autour des Ryūkyū. Commençons par faire l’union de nos polygones d’Okinawa et d’Amami avec la fonction gUnion() du module rgeos. Ensuite extrayons les coordonnées de ses extrémités sous forme de polygone avec la fonction gEnvelope(). Nous pouvons alors simplement cartographier ce polygone (vu la petite taille des données, pas besoin d’utiliser fortify, on peut laisser ggplot2 le faire à la volée).

require(rgeos)
rk <- gUnion(ama, oki)
rkLim <- gEnvelope(rk)
ggmap(asia) + 
  geom_polygon(data = rkLim, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .25) +
  theme_map()

La forme rectangulaire ne rend cependant pas bien compte de l’archipel. Nous allons à la place créer le polygone correspondant à l’enveloppe convexe de notre ensemble de polygones des Ryūkyū.

rkConv <- gConvexHull(rk)
ggmap(asia) + 
  geom_polygon(data = rkConv, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .25) +
  theme_map()

On peut également agrandir cette enveloppe avec un tampon (gBuffer()).

rkBuff <- gBuffer(rkConv, width = .5)
ggmap(asia) + 
  geom_polygon(data = rkBuff, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .25) +
  theme_map()

Annotons notre zone des Ryūkyū avec du texte ancré au centroïde de ce polygone, et affiché légèrement décalé vers la droite.

rkC <- as.data.frame(gCentroid(rkBuff))
ggmap(asia) + 
  geom_polygon(data = rkBuff, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .25) +
  geom_text(data = rkC, aes(x = x, y = y), label = "Iles Ryūkyū", hjust = -.4) +
  theme_map()

Pavage en polygones de Voronoï

La partition d’une surface en polygones de Voronoï est particulièrement utile pour les cartes géolinguistiques. Cette méthode permet d’obtenir à partir d’un ensemble de points un ensemble de polygones recouvrant l’ensemble des points les plus proches de chaque point de référence.

Prennons un exemple simple en considérant que nous avons les coordonnées de 7 lieux à Kyūshū.

kyushuP <- data.frame(id = c("F","O","N","M","S","Q","K"),
                   long = c(130.41,131.61,129.87,131.42,130.3,130.74,130.55),
                   lat = c(33.6,33.24,32.74,31.91,33.25,32.79,31.56))

Affichons-les sur une carte. Pour la suite nous avons besoin de fusionner les polygones des préfectures de Kyūshū. Pour cela nous créons l’union des polygones.

jpn2@data$south <- coordinates(jpn2)[,2] < 30.98
kyushu <- subset(jpn2, NAME_1 %in% c("Fukuoka", "Oita", "Kumamoto", "Saga", "Miyazaki", "Kagoshima", "Naoasaki") & south == FALSE)
kyushu <- gUnaryUnion(kyushu)
kyushuF <- fortify(kyushu)
ggplot() +
  geom_polygon(data = kyushu, aes(x = long, y = lat, group = group), fill = "grey50", color = NA) +
  geom_point(data = kyushuP, aes(x = long, y =lat), color = "red") +
  coord_fixed() +
  theme_map()

Calculons maintenant la partition de Voronoï à partir de ces points. Nous allons pour cela utiliser la fonction voronoi() du module dismo, qui prend comme argument un objet de classe SpatialPoints ou une matrice de coordonnées à deux colonnes.

require(dismo)
vor <- voronoi(kyushuP[,2:3], ext = extent(kyushu))
vorF <- fortify(vor)
ggplot() +
  geom_polygon(data = kyushu, aes(x = long, y = lat, group = group), fill = "grey50", color = NA) +
  geom_polygon(data = vorF, aes(x = long, y = lat, group = group), fill = NA, color = "orange") +
  geom_point(data = kyushuP, aes(x = long, y =lat), color = "red") +
  coord_fixed() +
  theme_map()

Il suffit ensuite de créer l’intersection de l’ensemble de poygones de Voronoï avec les polygones de Kyūshū.

vorK <- gIntersection(vor, kyushu, byid = TRUE)
vorKF <- fortify(vorK)
ggplot() +
  geom_polygon(data = kyushu, aes(x = long, y = lat, group = group), fill = "grey50", color = NA) +
  geom_polygon(data = vorKF, aes(x = long, y = lat, group = group), fill = NA, color = "orange") +
  geom_point(data = kyushuP, aes(x = long, y =lat), color = "red") +
  coord_fixed(ylim = c(30.99, max(kyushuF$lat))) +
  theme_map()

On peut maintenant colorier chaque polygone de la partition.

ggplot() +
  geom_polygon(data = kyushu, aes(x = long, y = lat, group = group), fill = "grey50", color = NA) +
  geom_polygon(data = vorKF, aes(x = long, y = lat, group = group, fill = id), color = NA) +
  geom_point(data = kyushuP, aes(x = long, y =lat), color = "red") +
  coord_fixed(ylim = c(30.99, max(kyushuF$lat))) +
  theme_map() + guides(fill = FALSE)

À suivre dans la 7e partie.

Ce billet peut être cité de la façon suivante: Thomas Pellard, « Réaliser une carte avec R (6) », Cipanglossia, 28 juin 2017, https://cipanglo.hypotheses.org/475.

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *