Réaliser une carte avec R (2)

La première partie du tutoriel nous a permis de réaliser une carte descriptive simple à partir d’images de fonds de cartes téléchargées depuis un serveur en ligne. Dans cette deuxième partie, nous allons apprendre à utiliser des données de systèmes d’informations géographiques (SIG) afin de réaliser des cartes vectorielles affichant des divisions administratives.

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.

Mettre en valeur une région sur une carte

Prennons pour premier exemple le cas où l’on souhaite réaliser une carte simple, sans décorations topographiques, en noir et blanc du Japon dans son contexte est-asiatique. On peut utiliser les fonctions vues dans la première partie en récupérant une carte Stamen de type toner-lite.

require(ggmap)
xlims <- c(110,150)
ylims <- c(20,50)
jpnMap <- get_stamenmap(bbox = c(left = min(xlims), bottom = min(ylims), right = max(xlims), top = max(ylims)), maptype = "toner-lite", zoom = 4)
ggmap(jpnMap)

Imaginons maintenant que l’on veuille mettre en valeur le Japon ou une région du Japon. C’est impossible car le fond de carte est une image sans structure, il n’y a aucune information géographique à part les coordonnées des limites de la carte. Pour remédier à ce problème il faut des données SIG, plus précisément des fichiers de polygones correspondant aux divisions géographiques que l’on veut représenter.

Heureusement le module mapdata contient des données de résolution relativement haute pour le monde entier. Le fond de carte se trace avec la fonction borders à laquelle on passe l’argument "worldHires" (pour avoir des données de bonne résolution), les couleurs désirées, et l’épaisseur du trait avec size. Il faut ensuite faire appel à coord_quickmap() qui prend à nouveau les spécifications de limites. Mettons la carte en mémoire dans une variable asiaMap.

require(mapdata)
asiaMap <- ggplot() +
  borders("worldHires", fill = "gray", colour = "black", xlim = xlims, ylim = ylims, size = .25) +
  coord_quickmap(xlim = xlims, ylim = ylims)
asiaMap

On peut supprimer le fond gris et les lignes en utilisant un thème comme theme_classic().

asiaMap <- asiaMap + theme_classic()
asiaMap

Mettons maintenant le Japon en rouge. Pour cela il suffit de récupérer à part uniquement les polygones correspondant au Japon. Ensuite on les ajoute à la carte avec la fonction geom_polygon().

jpn <- map_data('worldHires', region = 'Japan')
asiaMap + geom_polygon(data = jpn, aes(x = long, y = lat, group = group), color = NA, fill = "red")

On peut améliorer le résultat en jouant sur la transparence avec l’option alpha.

asiaMap + geom_polygon(data = jpn, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .45)

On peut également spécifier des sous-régions.

unique(jpn$subregion)
##  [1] "Haha Jima"         "Sado"              "Shikoku"          
##  [4] "Kakeroma Jima"     "Okino Daito"       "Okushiri To"      
##  [7] "Amakusa Shoto"     "Honshu"            "Kyushu"           
## [10] "Yome Jima"         "Muko Jima"         "Kume Shima"       
## [13] "Nishino Shima"     "Ishigaki Shima"    "Okinawa"          
## [16] "Miyake Jima"       "Tokuno Shima"      "Hokkaido"         
## [19] "Iwo Jima"          "Nakadori Jima"     "Hirado Jima"      
## [22] "Hachijo Jima"      "Kita Iwo Jima"     "Rebun To"         
## [25] "Suwanose Jima"     "Tanega Shima"      "Miyako Jima"      
## [28] "Kikaiga Shima"     "Rishiri To"        "Dogo"             
## [31] "Minami Daito Jima" "Chichi Jima"       "Iki"              
## [34] "Kamino Shima"      "Nakano Shima"      "Shimono Shima"    
## [37] "Kuchino Shima"     "Yaku Shima"        "Yonaguni Jima"    
## [40] "Iriomote Jima"     "Fukue Jima"        "Okino Erabu Shima"
## [43] "Amami O Shima"     "Tori Shima"        "Awaji Shima"      
## [46] "O Shima"

Choisissons par exemple de mettre en couleur l’île de Kyūshū.

ks <- subset(jpn, subregion == "Kyushu")
asiaMap + geom_polygon(data = ks, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .45)

Malheureusement, mapdata ne permet pas pour le Japon d’avoir accès à des subdivisions plus détaillées. Par ailleurs la résolution des données est très bonne pour des cartes à l’échelle des grands pays et des continents, mais elle est insuffisante pour les échelles plus petites. On peut le constater en traçant par exemple une carte d’Okinawa.

xlims <- c(126.65,128.4)
ylims <- c(26,27.2)
ggplot() +
  borders("worldHires", fill = "gray", colour = "black", xlim = xlims, ylim = ylims, size = .25) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On constate encore mieux les problèmes en superposant cette carte en rouge à une carte Stamen.

okiMap <- get_stamenmap(bbox = c(left = min(xlims), bottom = min(ylims), right = max(xlims), top = max(ylims)), maptype = "terrain-background", zoom = 10)
ggmap(okiMap) +
  borders("worldHires", fill = "red", colour = NA, alpha = .45, xlim = xlims, ylim = ylims, size = .25) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
   theme_classic()

Carte haute résolution

Il faut donc des données encore plus précises pour pouvoir travailler à des échelles plus petites. Heureusement de telles données sont disponibles gratuitement sur le site Global Administrative Areas. Il est possible de télécharger les données pour le monde entier (fichiers très volumineux) ou bien pays par pays. Choisissons la deuxième option. Plusieurs formats sont disponibles. Ceux qui nous intéressent sont Shapefile ou bien R (SpatialPolygonsDataFrame). Choisissons le premier, plus universel.

On obtient alors un fichier JPN_adm_shp.zip qu’il suffit de décompresser. Il y a trois jeux de données, 0, 1, et 2, qui sont de précision croissante dans le niveau de subdivision administrative:

  1. pays
  2. préfectures
  3. municipalités

Choisissons le niveau 2. Il faut faire appel au module rgdalet à la fonction readOGR() qui prend comme arguments l’adresse du fichier .shp que l’on souhaite utiliser et le nom de fichier sans extension. On peut tout d’abord inspecter son contenu.

require(maptools)
require(rgdal)
jpn2 <- readOGR(dsn = path.expand("~/repos/cartes/jpn/adm/jpn"), layer = "JPN_adm2") # à modifier en fonction de l'endroit où se trouve le fichier sur l'ordinateur
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/tpellard/repos/cartes/jpn/adm/jpn", layer: "JPN_adm2"
## with 1811 features
## It has 11 fields
head(jpn2@data)
##   ID_0 ISO NAME_0 ID_1 NAME_1 ID_2        NAME_2 NL_NAME_2 VARNAME_2
## 0  114 JPN  Japan   28   Nara 1199         Sango       ???      <NA>
## 1  114 JPN  Japan   28   Nara 1200     Shimoichi       ???      <NA>
## 2  114 JPN  Japan   28   Nara 1201 Shimokitayama      ????      <NA>
## 3  114 JPN  Japan   28   Nara 1202          Soni       ???      <NA>
## 4  114 JPN  Japan   28   Nara 1203      Takatori       ???      <NA>
## 5  114 JPN  Japan   28   Nara 1204    Tawaramoto      ????      <NA>
##   TYPE_2 ENGTYPE_2
## 0  Machi      Town
## 1  Machi      Town
## 2   Mura   Village
## 3   Mura   Village
## 4  Machi      Town
## 5  Machi      Town
unique(jpn2$NAME_1)
##  [1] Nara      Niigata   Yamanashi Aichi     Akita     Aomori    Chiba    
##  [8] Ehime     Fukui     Fukuoka   Fukushima Gifu      Gunma     Hiroshima
## [15] Hokkaido  Hyogo     Ibaraki   Ishikawa  Iwate     Kagawa    Kagoshima
## [22] Kanagawa  Kochi     Kumamoto  Kyoto     Mie       Miyagi    Miyazaki 
## [29] Nagano    Naoasaki  Oita      Okayama   Okinawa   Osaka     Saga     
## [36] Saitama   Shiga     Shimane   Shizuoka  Tochigi   Tokushima Tokyo    
## [43] Tottori   Toyama    Wakayama  Yamagata  Yamaguchi
## 47 Levels: Aichi Akita Aomori Chiba Ehime Fukui Fukuoka Fukushima ... Yamanashi
head(unique(jpn2$NAME_2))
## [1] Sango         Shimoichi     Shimokitayama Soni          Takatori     
## [6] Tawaramoto   
## 1668 Levels: Abashiri Abiko Abira Abu Achi Adachi Aga Agano ... Zushi

On peut facilement extraire les données des îles Okinawa proprement dites. La fonction fortify() transforme les données de polygones shp vers un format plus facile à traiter par ggplot2.

oki <- subset(jpn2, NAME_1 == "Okinawa" & !(NAME_2 %in% c("Minami Daito","Ishigaki","Taketomi","Miyakojima","Yonaguni","Tarama")))
unique(oki$NAME_2)
##  [1] Ogimi          Aguni          Chatan         Ginowan       
##  [5] Ginoza         Haebaru        Higashi        Ie            
##  [9] Iheya          Itoman         Izena          Kadena        
## [13] Kin            Kitanakagusuku Kumejima       Kunigami      
## [17] Motobu         Nago           Naha           Nakagusuku    
## [21] Nakijin        Nanjo          Nishihara      Okinawa       
## [25] Onna           Tokashiki      Tomigusuku     Tonaki        
## [29] Unknown5       Unknown8       Urasoe         Uruma         
## [33] Yaese          Yomitan        Yonabaru       Zamami        
## 1668 Levels: Abashiri Abiko Abira Abu Achi Adachi Aga Agano ... Zushi
okiF <- fortify(oki)

Nous voilà maintenant prêts à réaliser notre carte haute résolution d’Okinawa.

ggplot() + 
  geom_polygon(data = okiF, aes(x = long, y = lat, group = group), color = NA, fill = "grey") +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On peut vérifier la qualité des données en superposant à nouveau cette carte en rouge sur la carte Stamen.

ggmap(okiMap) +
  geom_polygon(data = okiF, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .5) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On peut également zoomer sur une petite île pour constater le degré de précision. Par exemple les îles Tonaki.

tonaki <- subset(oki, NAME_2 == "Tonaki")
tonakiF <- fortify(tonaki)
ggplot() +
  geom_polygon(data = okiF, aes(x = long, y = lat, group = group), color = NA, fill = "grey") +
  geom_polygon(data = tonakiF, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .45) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

ggplot() +
  geom_polygon(data = tonakiF, aes(x = long, y = lat, group = group), color = NA, fill = "grey") +
  coord_quickmap() +
  theme_classic()

Lignes de côtes

Attention nous n’avons pas tracé de lignes de côtes (color = NA) intentionnellement. Les données étant au niveau des municipalités, on obtiendrait le résultat suivant.

ggplot() + 
  geom_polygon(data = okiF, aes(x = long, y = lat, group = group), color = "black", fill = "grey", size = .25) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

Pour des lignes de côtes, il faut ici utiliser le niveau 0.

jpn0 <- readOGR(dsn = path.expand("~/repos/cartes/jpn/adm/jpn"), layer = "JPN_adm0")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/tpellard/repos/cartes/jpn/adm/jpn", layer: "JPN_adm0"
## with 1 features
## It has 24 fields
jpn0F <- fortify(jpn0)
ggplot() + 
  geom_polygon(data = jpn0F, aes(x = long, y = lat, group = group), color = "black", fill = "grey", size = .25) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On peut bien sûr rajouter un coloriage.

ggplot() + 
  geom_polygon(data = jpn0F, aes(x = long, y = lat, group = group), color = "black", fill = "grey", size = .25) +
  geom_polygon(data = okiF, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .5) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

Retour à une échelle plus grande

On peut également revenir à une échelle plus grande et colorier une préfecture japonaise, ici Kagoshima.

jpn1 <- readOGR(dsn = path.expand("~/repos/cartes/jpn/adm/jpn"), layer = "JPN_adm1")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/tpellard/repos/cartes/jpn/adm/jpn", layer: "JPN_adm1"
## with 47 features
## It has 9 fields
kag <- subset(jpn1, NAME_1 == "Kagoshima")
kagF <- fortify(kag)
asiaMap + geom_polygon(data = kagF, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .65)

On peut aussi colorier Kagoshima sur une carte Stamen.

xlims <- c(110,150)
ylims <- c(20,50)
jpnMap <- get_stamenmap(bbox = c(left = min(xlims), bottom = min(ylims), right = max(xlims), top = max(ylims)), maptype = "terrain-background", zoom = 5)
ggmap(jpnMap) + 
  geom_polygon(data = kagF, aes(x = long, y = lat, group = group), color = NA, fill = "red", alpha = .65)

À suivre dans la troisième partie.

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

Une réflexion au sujet de « Réaliser une carte avec R (2) »

  1. Ping : Réaliser une carte avec R (3) | Cipanglossia

Laisser un commentaire

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

Ce site utilise Akismet pour réduire les indésirables. En savoir plus sur comment les données de vos commentaires sont utilisées.