Aujourd’hui je vais vous montrer comment, en utilisant R, on peut faire de belles cartes géographiques.
Et là, vous allez me demander, mais pourquoi faire des cartes géographique ? Et pourquoi avec R ?
Et bien imaginons que, vous, bioinformaticien de terrain, soyez allé échantillonner des animaux à l’autre bout du monde sur plusieurs sites, par exemple des Marsupilami (totalement au hasard !). Vous voulez faire une carte de ces différents sites d’échantillonnage.
Facile ! Il suffit d’utiliser Google Earth, et d’y ajouter les points me direz-vous.
Oui, mais voilà, vous avez 500 points d’échantillonnage. Et 500, ça commence à faire beaucoup à faire à la main… Et puis votre chef étant un ayatollah du libre, vous venez d’être viré par le simple fait d’y avoir pensé !
Et le gros avantage d’utiliser R sera de pouvoir utiliser toutes ses fonctions graphiques sur votre carte (et puis c'est libre !). Nous verrons, à la fin de cet article, comment ajouter des graphiques sur une carte.
Pour faire une simple carte, on va utiliser les packages maps[1] et mapdata[2]. Une fois installés, nous utiliserons la fonction
1 |
map() |
pour créer notre carte. Nous utiliserons les données de la base de données worldHires fournie dans le package mapdata.
Le code suivant donne une carte du monde.
1 2 3 4 |
library(maps) library(mapdata) map('worldHires') |
Les options graphiques de R peuvent s’appliquer. On peut choisir de ne dessiner qu'une partie du monde en utilisant les options xlim et ylim respectivement pour régler la longitude et la latitude. Il faut donc connaître les coordonnées géographique des quatre coins de la carte qui nous intéresse (elles peuvent être trouvées grâce à la fonction
1 |
locator() |
de R). L’option color permettra de modifier la couleur des frontières entre pays, et utilisée avec l’option fill, on pourra colorier les pays.
1 2 3 4 5 6 7 |
map( 'worldHires', col=rainbow(18), fill=T, xlim=c(-19,60), ylim=c(-40,40) ) |
Ouais, mais c’est chiant ton truc. Il faut connaître les coordonnées géographiques dans le système décimal pour dessiner la carte. On peut pas demander le pays que l’on veut ?
Et si ! On peut spécifier le pays que l’on veut dessiner. Par exemple, la commande suivante permet de dessiner une carte du Japon avec un fond grisé.
1 2 3 4 5 6 |
map( 'worldHires', "japan", col='gray80', fill=T ) |
Puisqu’une carte sans échelle ne veut rien dire (comment j’ai bien retenu mes cours de géographie du collège !), on utilisera la fonction
1 |
map.scale() |
pour l’ajouter.
On peut ensuite ajouter les villes sur cette carte, grâce à la fonction
1 |
map.cities() |
.
On pourra ensuite ajouter des points sur la carte à partir des coordonnées géographiques. Il faudra alors disposer de ces coordonnées géographiques dans le système décimal. Si vous avez les coordonnées dans le système sexagésimal, vous pouvez les convertir, par exemple sur ce site : https://tools.wmflabs.org/geohack/
Les coordonnées de Kyoto dans le système sexagésimal sont 35°40’14.6“N, 139°46’18.86“E (merci Wikipedia : https://fr.wikipedia.org/wiki/kyoto). Ces coordonnées donnent dans le système décimal : 35.670724°N pour la latitude et 139.771907°E pour la longitude.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
map('worldHires', "japan", col='gray80', fill=T) map.scale( 134, 26, metric=T, relwidth=0.3 ) map.cities( country='Japan', capitals=1, pch=15, col='red' ) points( 135.7, 35, pch=16 ) text( 135.7, 35.3, label="Kyoto" ) points( 140.5, 37.8, pch=16 ) text( 140.5, 38.2, label="Fukushima" ) |
Je vois d'ici les petits malins qui ont voulu, pour tester, faire une belle carte de la France métropolitaine, et qui n'ont pas été super contents du résultat…
En effet, la commande
1 2 3 4 |
map( "worldHires", "france" ) |
représente la France… en entier ! Donc avec les DOM-TOM (enfin les DROM-COM pour ceux qui ont suivi les changements d'acronymes).
Il faudra donc, dans le cas de la France métropolitane, passer obligatoirement par les coordonnées géographiques :
1 2 3 4 5 6 |
map( "worldHires", "france", xlim=c(-5,10), ylim=c(35,55) ) |
Maintenant revenons à nos marsupilamis. On va réaliser une carte de leur aire de répartition sur laquelle on ajoutera les points d’échantillonnage. Pour les aires de répartition d’une espèce, il est "assez facile" de trouver des données sur internet sous forme de shapefile (.shp) qui sont ceux que l’on utilise pour faire une carte.
On représente les aires de répartition de deux types de population de la même espèce de marsupilami (Marsupilami fantasii) : les marsupilamis jaunes à tâches (dont l'aire de répartition est représentée en vert) et les marsupilamis jaunes uniformes (en rouge). Sur la carte, on voit que les deux populations cohabitent sur une aire géographique limitée.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
library(maptools) repartition = readShapePoly('./repartition_marsu.shp’) map( 'worldHires', xlim=c(-90,-35), ylim=c(-60,15) ) plot( repartition[55,], add=TRUE, col=rgb(20,84,30,150,maxColorValue=255), border=FALSE ) plot( repartition[84,], add=TRUE, col=rgb(201,21,34,150,maxColorValue=255), border=FALSE ) |
Pour faire cette carte, j’ai importé une nouvelle libraire (maptools) dans laquelle se trouve la fonction
1 |
readShapePoly() |
qui permet d’ouvrir des fichiers « de forme ». J’ai récupéré le fichier concernant les "caecilian amphibians"[3] (mais ! On bossait pas sur le marsupilami ?!) sur le site de l’UICN Red List (Liste rouge des espèces en danger)[4]. Une fois ouvert avec R, vous devriez obtenir une data frame contenant une aire de répartition par ligne (et j'ai bien galéré à en trouver deux qui se recoupent).
Une légende pourrait être ajoutée en utilisant la fonction
1 |
legend() |
de R.
Maintenant intéressons-nous en particulier à la zone où les deux populations cohabitent. Sur cette carte, j’ai ajouté les treize points d’échantillonnage dans cette région (dont sept dans la zone de cohabitation). On peut ajouter ces points grâce à la fonction
1 |
points() |
en ayant au préalable importé le fichier csv qui les contient.
Pour zoomer, j’ai cherché les coordonnées des quatre points les plus, respectivement, au nord, à l’ouest, au sud et à l’est de la zone qui m’intéresse grâce à la fonction
1 |
locator() |
de R.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
map( 'worldHires', xlim=c(-65.5,-50), ylim=c(-14,-8) ) plot( repartition[84,], add=TRUE, col=rgb(201,21,34,150,maxColorValue=255), border=FALSE ) plot( repartition[55,], add=TRUE, col=rgb(20,84,30,150,maxColorValue=255), border=FALSE ) points( coord$x, coord$y, pch=16 ) |
C’est bien beau d’échantillonner des bestioles sur le terrain, mais normalement, on n'échantillonne pas pour le plaisir d'échantillonner, mais pour faire quelque chose avec cet échantillonnage… On a donc séquencé le cytochrome b des différents individus échantillonnés dans ces différentes populations, ce qui nous a permis d’identifier huit haplotypes différents.
Hey m4rsu ! On pourrait pas faire une représentation des différents haplotypes identifiés sur chaque site et leur fréquence, genre avec un camembert ?
J'allais justement y venir !
On peut donc ajouter des graphiques en camembert sur une carte grâce à la fonction
1 |
draw.pie() |
qui est fournie avec la librarie mapplots. Il faut disposer d’un fichier qui comporte quatre colonnes : longitude, latitude, haplotype, nombre d’observations. On pourra alors utiliser la fonction
1 |
make.xyz() |
qui créée automatiquement un objet utilisable dans la fonction
1 |
draw.pie() |
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 |
library(mapplots) xyz = make.xyz( haplotypes$x, haplotypes$y, haplotypes$freq, haplotypes$haplo ) map( 'worldHires', xlim=c(-65.5,-50), ylim=c(-14,-8) ) plot( repartition[84,], add=TRUE, col=rgb(201,21,34,150,maxColorValue=255), border=FALSE ) plot( repartition[55,], add=TRUE, col=rgb(20,84,30,150,maxColorValue=255), border=FALSE ) draw.pie( xyz$x, xyz$y, xyz$z, radius=0.3, col=rainbow(8) ) legend( 'topright', legend=c(1 :8), col=rainbow(8), pch=15, ncol=2 ) |
Tadaa ! Bon, je vous laisse faire l’interprétation (bon courage, j'ai inventé les données !), car ce n’est pas l’objet de cet article.
Nous avons donc vu comment faire des cartes avec R. N’étant pas un spécialiste, je vous ai montré quelques possibilités, mais il y en a plein d’autres, et je suis sûr que l’on peut faire des cartes bien plus belles que celles que j’ai montrées ici. Il existe plein d'autres librairies pour faire des cartes, vous devriez y trouver votre bonheur !
Pour finir, je voudrais préciser qu’aucun marsupilami n’a été maltraité pour écrire cet article. En plus, c'est un gros fake : tout le monde sait que les marsupilamis vivent dans la forêt palombienne et non amazonienne.
Merci à mathurin, hedjour, Olivier Dameron et Sylvain P. pour leurs précieux conseils et la relecture de cet article. Merci également à Kim Gilbert pour son article "Making maps with R"[5] qui m'a bien aidé à mes débuts avec les cartes.
Références
[1] Package maps : http://cran.r‑project.org/web/packages/maps/ [2] Package mapdata : http://cran.r‑project.org/web/packages/mapdata/index.html [3] : shapefile des Caecilian Amphibians (ou Gymnophonia) : http://goo.gl/OFGYl1 [4] : IUCN Red List http://www.iucnredlist.org/technical-documents/spatial-data [5] : Making map with R : http://www.molecularecologist.com/2012/09/making-maps-with‑r/
Laisser un commentaire