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:
- pays
- préfectures
- municipalités
Choisissons le niveau 2. Il faut faire appel au module rgdal
et à 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.
1 réponse
[…] la première et la deuxième partie du tutoriel nous avons vu comment réaliser des cartes de qualité et comment mettre en valeur des régions. […]