Nous adressons toutes nos pensées à la famille de notre ami Jérôme !
http://www.forumsig.org/showthread.php/43488-Disparition-de-Phoenix
Nous adressons toutes nos pensées à la famille de notre ami Jérôme !
http://www.forumsig.org/showthread.php/43488-Disparition-de-Phoenix
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).
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.:
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 ?
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.
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):
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).
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:
Cette dernière démarche n'utilise pas r.profile. Son principe est le suivant:
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."
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
Sauf mention contraire dans les contenus, l'ensemble de ce site relève de la législation française et internationale sur le droit d'auteur et la propriété intellectuelle.
Le portailSIG est édité par l'association loi 1901 Forum Systèmes d'Information Géographique
dont le siège social se situe à AMIENS
ISSN 2274-4150
Commentaires
Poster un nouveau commentaire