Grass Gis et R: superposition des couleurs d'un raster quelconque sur un profil topographique élaboré à partir d'un MNT

Sur le ForumSIG, a été posée une question qui peut être résumée de la manière suivante:

Comment superposer sur un profil topographique, élaboré à partir d'un MNT, les couleurs provenant de n'importe quel autre raster ?

J'ai été amené à proposer une solution avec la combinaison des logiciels Grass Gis et R. C'est celle-ci qui va être développée dans la suite.

La résolution de ce problème peut en effet être très intéressante pour de nombreuses applications pratiques avec des cartes d'occupation du sol, en pédologie avec les toposéquences, ou avec des cartes géologiques.

Le cas des vecteurs ne sera pas abordé ici (grass.osgeo.org/wi/How_to_create_an_elevation_profile).

Profil interactif (ou d.profile)

Grass possède un outil interactif pour créer des profils topographiques (et autres). La ligne de profil est directement tracée sur le moniteur (ici ligne bleue, sur un fragment de carte géologique superposé à un MNT). 

ligne de profil dynamique

Une fois le profil tracé, Grass autorise le choix du raster à traiter.:

• Dans le cas d'un raster 3D comme un MNT,  l'altitude - élévation sera tracée en y et la distance cumulée en x: 

• Dans le cas d'un raster simple, non 3D (ici carte géologique), la valeur tracée en y sera généralement la catégorie (cat):

Les profils obtenus peuvent être exportés au format eps. Comment alors les superposer, c'est-à-dire avoir les couleurs de la carte géologique sur le profil topographique ?

Première solution: r.profile

La première solution est l'utilisation de la fonction r.profile qui permet de créer et de traiter un profil de manière programmée et de récupérer les données. Les coordonnées des points peuvent être entrées de manière interactive dans le moniteur ou en fournissant directement les coordonnées x, y:

Par exemple, dans le cas d'une ligne limitée par 2 noeuds (x1, y1, x2, y2), la commande:

r.profile input=MNT  output=profile.pt profile=171200,118000,171200,108000

donnera le résultat sous la forme de points ( on peut rajouter res=x, un point tous les x unités de mesure, ou c'est la résolution de la région qui sera utilisée) sous la forme distance cumulée, altitude-élévation pour un raster 3D (ici avec une résolution de région de 30 m et le fichier profile.pt comme résultat):

    0.000000 206.000000
   30.000000 209.000000
   60.000000 209.000000
   [....]

Pour les autres rasters (non 3D), le résultat sera sous la forme distance cumulée, catégorie (généralement la variable cat):

    0.000000 12
   30.000000 12
   60.000000 12
   [....]

En pratique, ce sont les mêmes résultats que ceux obtenus avec la commande de profil interactif, sous une forme plus facile à traiter, mais avec la possibilité d'ajouter d'autres paramètres comme l'ajout du flag -c. Celui-ci permet d'obtenir les couleurs des points sur les rasters utilisés, sous la forme RRR.GGG.BBB (RGB):

r.profile -c input=MNT output=profile.pts profile=171200,118000,171200,108000

   0.000000 206.000000 228:066:050
  30.000000 209.000000 228:066:050
  60.000000 209.000000 228:066:050
  [....]

Toute ligne droite sans segments (2 noeuds seulement) peut être utilisée. S’il y a plusieurs segments (et plusieurs noeuds), des problèmes seront rencontrés (surtout dans le calcul de la distance cumulée). La solution est alors de transformer l'ensemble en polyligne, comme le montre http://osgeo-org.1803224.n2.nabble.c...td5298928.html.

En conclusion, si l'on travaille dans une même région, les valeurs des distances cumulées seront toujours les mêmes quelque soit le raster. Les seules variables qui vont changer sont les paramètres z ou cat et les couleurs RGB.

Application pratique

La ligne de profil est numérisée (ligne rouge) sur l'ensemble MNT-carte géologique. L'objectif final est d'obtenir une ligne de profil topographique (MNT) coloriée avec les couleurs de la carte géologique.

 ligne de profil numérisée (moniteur et Nviz)

Pour envoyer les paramètres de la ligne dans r.profile, plusieurs solutions sont possibles (voir http://osgeo-org.1803224.n2.nabble.c...td5298928.html.). L'une d'entre elles est de transformer la ligne en points puis de les traiter (à l'aide de « pipes », fr.wikipedia.org/wiki/Tube_Unix):

• transformation de la ligne en points: v.to.points  avec le flag -i (interpolation)
• exportation des points en fichier texte: v.out.ascii
• utilisation de ces points dans r.profile avec la commande:

cat  sortie_de_v.out.ascii | r.profile -c input=raster output=profil.pts

Une autre est de traiter directement la ligne, sans passer par la phase d'exportation des points avec la commande:

v.out.ascii in=ligneprofil format=standard | grep  '^ ' | r.profile -c in=raster output=profil.pts

Ce dernier processus est appliqué au MNT et à la carte géologique superposée.

# création des points avec couleur RGB du MNT
GRASS 6.4 (geol):~ > v.out.ascii in=ligneprofil format=standard | grep  '^ '|\ 
r.profile -c in=MNT res=30 output=profil.pts
 
# création des points avec couleur RGB de la carte géologique
GRASS 6.4 (geol):~ > v.out.ascii in=ligneprofil format=standard | grep  '^ '|\
 r.profile -c in=geol res=30 output=profiligeol.pts

Deux fichiers sont obtenus, l'un avec l'altitude-élévation et les couleurs des points sur le MNT et l'autre avec les couleurs de ces mêmes points sur la carte géologique ( la catégorie ne nous intéresse pas ici).

La suite du traitement consiste donc à combiner les 2 fichiers pour obtenir comme variables à traiter la distance cumulée, l'élévation et la couleur RGB des points sur la carte géologique, puis de dessiner le profil obtenu.

La démarche pourrait être faite avec de nombreux logiciels, comme un tableur ou gnuplot, mais j'ai choisi de le faire avec R qui s'intègre très bien avec Grass, comme cela a déjà été montré sur ce Portail. L'utilisation de R pour créer les profils à partir des points obtenus dans Grass a déjà été illustrée comme nous le verrons par la suite.

La solution sera exposée ici de manière « basique » avec les fonctions standards de R. La première démarche est de transformer le format RRR:GGG:BBB en 3 colonnes r, g et b, pour pouvoir les traiter dans R (c'est-à-dire de 228:066:050 à 228 066 050).

Le processus est exécuté directement avec sed:

GRASS 6.4 (geol):~ >sed -i -e "s/\:/ /g" profil.pts
GRASS 6.4 (geol):~ >sed -i -e "s/\:/ /g" profilgeol.pts

R est lancé dans le shell Grass et le premier fichier est traité:

GRASS 6.4 (geol):~ > R
[....]
 
# lecture de la table "profil MNT"
> m1 <- read.table('profil.pts', header=F, sep=" ",
+ col.names=c('distance','elev','r','g','b'))
 
# plot des points
> plot(m1$distance,m1$elev,pch=21,xlab='Distance (m)', ylab='Elevation (m)',
+  main='Profil brut)

R permet de colorier les points-cercles (paramètre pch=21) avec des couleurs RGB. Les couleurs des points sur le MNT sont d'abord utilisées (colonnes r, g, b):

> plot(m1$distance,m1$elev,pch=21,col = rgb(m1$r/255,m1$g/255,m1$b/255),
+ bg = rgb(m1$r/255,m1$g/255,m1$b/255),xlab='Distance (m)', 
+ ylab='Elevation (m)', main='Profil brut (couleur MNT)')

R permet de combiner 2 objets ayant une même colonne commune (ici distance cumulée) avec la fonction merge, ce qui permet de traiter directement les couleurs des points sur la carte géologique :

# lecture de la table "profil géologie"
> m2 <- read.table('profilgeol.pts', header=F, sep=" ", col.names=c('distance',
+ 'cat','r','g','b'))
 
# fusion des 2 tables sur base de la variable distance qui est la même:
# le data frame résultant aura donc les paramètres distance et élévation de la table 
# "profil MNT" et les couleurs r,g,b de la table "profil géologie"
> res <- merge(m1[,c('distance','elev')],m2[,c('distance','r','g','b')],
+ by="distance", all=T)
 
#plot
> plot(res$distance,res$elev,pch=21,col = rgb(res$r/255,res$g/255,res$b/255),
+ bg = rgb(m2$r/255,m2$g/255,m2$b/255),xlab='Distance (m)', ylab='Elevation (m)', 
+ main='Profil géologie')

Comment obtenir une ligne de profil coloriée et non plus les points coloriés seuls ? L'ajout du paramètre ty="o" permet de dessiner une ligne reliant les points, mais ce n'est pas la solution.

> plot(res$distance,res$elev,pch=21,col = rgb(res$r/255,res$g/255,res$b/255),
+ bg = rgb(m2$r/255,m2$g/255,m2$b/255),xlab='Distance (m)', ylab='Elevation (m)', 
+ main='Profil géologie', ty="o")

La solution est de faire intervenir la fonction graphique segments de R. Celle-ci permet de relier 2 points par un segment de ligne, colorié ou non. La démarche est effectuée ici en assignant la couleur de point(i) au segment point(i)-point(i+1):

> plot.new()
> plot.window(xlim=range(m2$distance),ylim=range(m2$elev))
 
# plot des segments de ligne, lwd = épaisseur des segments de ligne
> i = 0
> while(i <= length(res$distance)-1){
+  i = i + 1
+  segments(res$distance[i],res$elev[i],res$distance[i+1],res$elev[i+1],
+ col=rgb(res$r[i]/255,res$g[i]/255,res$b[i]/255),lwd=3)
+ }
> axis(1); axis(2); box();
> title(main="Profil géologie", xlab="Distance (m)", ylab="Elevation (m)")

Si l''exagération verticale n'est pas souhaitée, il suffit de préciser asp=1 dans plot (tout est paramétrable, voir notamment www.statmethods.net/advgraphs/parameters.html).

Problème et deuxième solution (adaptée de celle de Dylan Beaudette)

Un problème se pose néanmoins dans les cas où les limites des ensembles cartographiés sont cruciales (tel est le cas d'une carte géologique avec les limites des formations cartographiées). En effet, r.profile génère les points exclusivement en fonction de la résolution d'une région ou à partir du paramètre res, sans s'occuper du reste. Diverses solutions sont alors possibles:

• jouer avec la résolution des points avec le paramètre res, sans aucune garantie que les points obtenus se trouvent sur les limites;
• placer sur la ligne de profil les points de contact des ensembles cartographiés, mais le résultat sera une ligne multinoeuds avec les problèmes déjà soulevés;
• modifier la démarche de Dylan Beaudette dans « Raster profile along arbitrary line segments » (http://casoilresource.lawr.ucdavis.edu/drupal/node/375) pour tenir compte des couleurs.

Cette dernière démarche n'utilise pas r.profile. Son principe est le suivant:

• transformation de la ligne en points: v.to.points  avec le flag -i (interpolation) ou choix des points (limites par exemple);
• transformation des points en points 3D avec la commande v.drape sur le MNT
• exportation des points 3D en fichier texte: v.out.ascii
• éventuellement représentation graphique de ces points 3D dans R: 

• et enfin utilisation de son script "Process profile data, and plot in R", que je vous laisse découvrir:

application de la méthode de Dylan Beaudette à notre exemple

Comment obtenir les couleurs de la carte géologique sur ce profil qui est obtenu à partir du MNT ? La solution est r.what qui permet d'obtenir pour un ou des points les valeurs de catégories ou les couleurs RGB d'un raster ou de rasters (ici MNT et carte géologique).

v.out.ascii in=pts_3d fs = ' '  | r.what -r input=MNT,geol > fichier

Le résultat sera un fichier avec les coordonnées x, y, z, la catégorie et les couleurs RGB de chaque point sur le MNT et la carte géologique :

167566.79603283|118063.59206566|211.60614082 1|211|255:185:000|12|228:066:050
167601.7495212|118034.80683995|212.17301181 1|212|255:181:000|12|228:066:050
[...]

Au final, un seul fichier à traiter et plus de limitation sur le nombre éventuel de segments. Une petite adaptation du script de Dylan Beaudette donne le résultat final:

Une petite remarque de l'auteur, néanmoins:

"Note that the elevation profiles are calculated in a rather non-standard way, as elevation is plotted against 3D distance along the original line segment. This means that resulting transect length calculations are sensitive to the resolution of the raster (or distance between sample points), across a rough surface: i.e. smaller grid spacing will result in a longer effective (3D) transect length."

Conclusions

Les graphiques R peuvent directement être tracés/exportés en divers formats, bitmap ou vectoriels qui peuvent être utilisés dans n'importe quel autre logiciel. Les graphiques pourraient certainement être améliorés avec l'utilisation de paquets spécifiques comme plotrix, ggplot2 ou autres (voir cran.r-project.org/web/views/Graphics.html)

L'ensemble du processus pourrait aussi faire l'objet d'une commande Grass « autonome », pour autant que R soit installé et disponible en ligne de commandes. La valeur cat (catégorie) pourrait aussi être utilisée (si le raster est bien « classé » ou usage de la commande r.reclass) en fixant par après des couleurs suivant la catégorie.

Et voilà, j'espère encore une fois avoir montré l'intérêt de la combinaison de deux logiciels tels que Grass et R. Charly, dans les réponses à la question sur le ForumSIG, a proposé une autre solution élégante avec gnuplot et je suis sûr que la démarche pourrait aussi être effectuée en Python avec matplotlib.

Tous les traitements ont été effectués sur Mac OS X avec Grass 6.4 et R version 2.9

Copyright (C) Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled "GNU Free Documentation License".