Grass Gis: région, résolution, résolution 3D, drapage, ou l'apparente quadrature du cercle...

La notion de région, propre à Grass est souvent mal comprise ou mal assimilée par les débutants, surtout lorsqu'ils travaillent dans Grass à partir de Quantum Gis ou lorsqu’ils veulent visualiser un raster en (pseudo) 3D avec NVIZ ou avec la nouvelle interface wxPython.

En pratique, une région n'est qu'une fenêtre de traitement (N,S,E,W) avec une certaine résolution. Une fois définie, tous les traitements postérieurs ne s'appliqueront exclusivement qu'à ce qui est dans la région et rien d'autre (pour les rasters uniquement et c'est très perturbant lorsqu'on débute sur Grass). Sa définition ne modifiera en rien les rasters originaux. Si l'un d'entre eux a, au départ, une résolution différente, Grass va effectuer un rééchantillonage temporaire à la volée de ce qui est visible en utilisant l'algorithme du voisin le plus proche (nearest neighbor resampling), sans que l'on s'en rende compte.

Quel est l'avantage d'une telle conception ? Une région peut être assimilée à une espèce de filtre entre les données brutes et les traitements que l'on veut effectuer dans une zone géographique donnée. Par exemple, dans mon cas, j'ai des secteurs (locations) avec les MNT SRTM ou Aster avec leurs emprises et leurs résolutions originales et je travaille dans des jeux de données (mapsets) avec plusieurs régions limitées en fonction de mes intérêts.

L'erreur la plus commune est de lancer des traitements avec une région mal adaptée ou mal définie,  ce qui au final ne donnera rien comme résultat (y compris pour les visualisations). Il en découle que si l'on veut visualiser dans Quantum Gis des données Grass "hors région", le résultat sera souvent une image désespérément blanche...

C'est la commande g.region [paramètres] (grass.osgeo.org/gdp/html_grass64/g.region.html) qui permet de créer une région. Sa forme la plus simple est:

g.region rast=unrasterparmidautres

Elle créera une région avec les caractéristiques du raster choisi (emprise et résolution) mais il est possible de le faire en modifiant directement les paramètres (voir dans la suite).

Pour expliquer plus clairement ces notions, nous allons détailler une démarche qui pose beaucoup de problèmes aux débutants (que j'étais il y a quelques années...), comment draper un raster sur un MNT avec des résultats corrects. Nous utiliserons un fragment scanné de carte, noir et blanc, que l'on veut superposer à un MNT, ici avec NVIZ.

g.region et emprise N, S, E, W d'un raster (et éventuellement utilisation de r.mask)

Le problème se présente de la manière suivante (la barre noire dans le coin supérieur gauche donne l'échelle, 700m):  l'emprise du raster noir et blanc (image scannée et géoréférencée) est beaucoup plus petite que celle du MNT (Aster de la NASA dans le cas présent).

Les propriétés des couches sont les suivantes (seuls les paramètres de résolution sont importants et seront détaillés ici, l'emprise des couches étant visible sur les figures):

Pour le MNT

Rows:            8369                                                      
Columns:      11854                                                      
Total Cells:   99206126                                                   
[...]                                
N: [...]    S:   [...]   Res: 19.91493275         
E:  [...]    W:  [...]  Res: 19.91493275          
Range of data:    min = 99.670000  max = 695.960022

Pour le raster (noir et blanc)

Rows:            885                                                        
Columns:      1130                                                       
Total Cells:   1000050                                                    
[...]                                 
N: [...]    S: [...]   Res: 1.00093415        
E: [...]    W: [...]  Res: 1.00093415      
Range of data:    min = 0  max = 1

Res indique qu'une cellule a une dimension de 1.00093415 x 1.00093415 unités de mesure (N S - EW ou yx). Pour comprendre comment ces paramètres sont calculés, voir www.portailsig.org/content/python-georeferencement-et-worldfiles-tfw-jgwpngw

Tant les étendues que les résolutions sont différentes et il parait logique de ne travailler que sur la zone de la couche raster et pas sur celle de l'ensemble du MNT. La commande

g.region rast=raster

va limiter une fenêtre de travail (visualisée par le carré rouge):

Après la commande, Grass n'effectuera donc les traitements postérieurs que sur les rasters ou les portions de rasters qui sont dans la région (avec les dimensions et la résolution définies). Tout ce qui se situe à l'extérieur sera ignoré.

Une autre solution plus "radicale " est d'utiliser un masque avec les commandes r.mask (grass.osgeo.org/gdp/html_grass64/r.mask.html) ou g.copy (grass.osgeo.org/gdp/html_grass64/g.copy.html) :

r.mask input=raster  ou g.copy rast=raster, MASK

Ces commandes auront pour effet de créer un masque du raster nommé MASK. Grass traitera uniquement les données situées à l'intérieur de la zone masquée, toutes les autres étant considérées comme si leurs valeurs étaient NULL. La résolution dépendra toujours de la région.

Ces 2 commandes permettent donc de créer des zones de travail très précises, "à la carte", pour effectuer les traitements.  Ces démarches sont valables pour l'emprise géographique d'une zone étudiée et permettent, par exemple, de découper très facilement le MNT par le raster. Mais qu'en est-il des résolutions originelles ? Le rééchantillonage temporaire à la volée paraîtra correct en 2D mais souvent horrible en 3D, surtout si les résolutions sont très différentes comme dans le cas présent.

g.region et résolutions des rasters: 3D (NVIZ) et problèmes

NVIZ permet de régler deux types de résolution, coarse et fine. Les paramètres sont assez faciles à comprendre, coarse s'adresse à la résolution de l'image "en mouvement" (déplacements, rotations, etc.), fine à l'image "fixe", mais quelle sont leurs relations avec la résolution d'une région choisie ?

fine = 1 signifie que la résolution choisie est celle préalablement fixée par g.region et non celle d'un raster particulier

Que va t'il se passer dans notre cas avec ces résolutions très différentes ?

MNT vu avec g.region fixée à partir du MNT originel (fine= 1 ou Res: 19.91493275 ):

MNT vu avec g.region fixée à partir du raster noir et blanc (fine=1 ou Res:1.00093415):

Les problèmes de résolutions différentes sont directement visibles par les effets "d'escalier" dus à la taille des pixels mais aussi à l'algorithme de rééchantillonage choisi par Grass, comme la suite le montrera.

Si l'on essaie de draper le raster noir et blanc sur le MNT, les résultats vont être encore plus catastrophiques:

Raster vu avec g.region fixée à partir du MNT originel (fine=1 ou Res: 19.91493275 ):

Raster vu avec g.region fixée à partir du raster noir et blanc (fine=1 ou Res:1.00093415):

Dans les 2 cas, la représentation est mauvaise, raster noir et blanc trop "grossier" (résolution pratiquement x 19 !) dans le premier cas et toujours MNT avec "escaliers" dans le deuxième, ce qui limite la qualité du rendu final.

g.region et résolutions des rasters: 3D (NVIZ) et solutions

Si le seul intérêt est l'amélioration de la vue 3D (dans NVIZ), la solution est de modifier définitivement (et pas temporairement comme pour la représentation 2D) la résolution d'une des couches pour la ramener à la première. Ce raster rééchantilloné ne sera par contre pas valide pour les traitements du fait des artéfacts éventuels créés par l'interpolation.

Dans le cas présent, ramener celle du raster noir et blanc à celle du MNT ne conduira qu'à sa représentation "grossière" déjà figurée. Il faut donc rééchantilloner le MNT mais il est évident qu'appliquer ce traitement à l'ensemble du MNT serait très couteux en temps et en mémoire machine.

C'est ici que la notion de région acquiert toute son importance puisqu'elle va permettre de ne travailler exclusivement que dans la fenêtre choisie. Si un MASK a été spécifié, la procédure sera plus rapide. Ce sera fait sans modifier le MNT originel en créant une copie locale de la zone étudiée à la résolution voulue. Il faut alors utiliser une des commandes de la famille r.resamp, r.resamp.interp, interpolation par algorithmes du voisin le plus proche, bilinéaire ou bicubique ou r.resamp.rst, interpolation par splines avec tension (redularized spline with tension).

1) solution avec r.resamp.rst:

La procédure est la suivante:

• créer une région avec l'emprise du raster noir et blanc et la résolution originale du MNT (19.91493275)
•  g.region n=[du raster NB] e=[du raster NB] s=[du raster NB] w=[du raster NB] res=[du MNT à rééchantilloner: 19.91493275]
• effectuer le rééchantillonage avec  la commande r.resamp.rst avec la résolution de 1.00093415 du raster noir et blanc, soit directement avec l'original du MNT, soit avec un MASK
  
•  r.resamp.rst input=mnt ew_res=1.00093415 ns_res=1.00093415 elev=mntres (nom du MNT résultant)
• revenir à la résolution du raster noir et blanc (1.00093415)
• g.region rast=raster

et le résultat obtenu pour mntres est:

Rows:         885                                                       
Columns:      1130                                                      
Total Cells:  1000050                                                   
[...]                               
N: [= celle du raster NB]    S: [= celle du raster NB] Res: [= celle du raster NB : 1.00093415]      
E: [= celle du raster NB]    W: [= celle du raster NB] Res: [= celle du raster NB : 1.00093415]
Range of data:    min = 333.307556  max = 421.131683

Maintenant que l'étendue et la résolution sont les mêmes,quels sont les résultats avec mntres ?

MNT ajusté seul avec fine=1 ou Res:1.00093415 :

raster noir et blanc drapé avec fine=1 ou Res:1.00093415:

Il restera à jouer avec les paramètres de luminosité pour enrichir le rendu.

2) solution avec r.resamp.interp

La démarche avec r.resamp.interp est légèrement différente puisqu'il faut cette fois-ci se placer dans une région avec l'emprise et la résolution du raster noir et blanc.

• créer une région avec l'emprise et la résolution du raster noir et blanc (1.00093415)
• g.region rast=raster
• effectuer le rééchantillonage avec la commande r.resamp.interp
• r.resamp.interp input=mnt output=mntinterp method=bilinear (ou bicubic ou nearest)

Les résultats obtenus avec les méthodes bilinear et bicubic sont équivalents à ceux obtenus avec r.resamp.rst, ceux avec nearest non (toujours l'effet escalier du à l'algorithme, ce qui explique les problèmes dus au rééchantillonage à la volée).

                method  = bilinear                                                  method = bicubic

method= nearest

Cette fois les résultats sont enfin conformes à ceux espérés.

retour à la case départ

Mais, que se passe-t-il si l'on modifie la valeur de fine en conservant la même résolution pour la région (1.00093415)  ?

 fine = 5  et donc Res = 5 x 1.00093415 = 5.00467075            fine = 10 et donc Res = 10 x 1.00093415 = 10.0093415

et finalement avec fine = 19 et donc Res = 19 x 1.00093415 = 19.01774885, ce qui correspond plus ou moins à la résolution du MNT originel

Comme on pouvait s'y attendre logiquement, le rendu du raster se dégrade peu à peu jusqu'à arriver à celui de la première représentation figurée, et la boucle est bouclée, car ceci explique les problèmes dus aux différences de résolutions.

Conclusions

J'ai essayé ici d'expliquer les démarches comme j'aurais voulu que l'on me les explique dès le départ, ce qui n'est pas évident au vu des longues recherches nécessaires pour trouver de la documentation claire dans le cas de NVIZ et la 3D. Les tutoriels sont souvent soit trop simples, voire simplistes, soit trop compliqués ou trop spécialisés. Heureusement, il existe les listes de l'OSGeo très réactives comme celles de Grass-User ou Grass-Dev qui permettent de poser des questions et d'obtenir des réponses rapides.

J'espère aussi que cet exemple montre l'importance d'une bonne compréhension des principes de Grass. Ils peuvent paraitre très différents de ceux d'autres SIG et ce n'est qu'une fois bien maitrisés que l'on pourra profiter de toutes ses potentialités. Il illustre aussi les erreurs que l'on peut commettre si on se lance directement à l'aventure (comme je l'ai fait dans le temps avec Grass v.5, avec les catastrophes que je vous laisse deviner...). C'est encore plus critique si on le fait à partir de Quantum Gis. 

Draper n'importe quel raster (pour autant que sa résolution de départ soit bonne !) sur un MNT en représentation 3D et ce, quelques soient leurs résolutions, n'a donc plus de secrets. Il faut néanmoins souligner que la démarche peut demander beaucoup de temps et de place sur le disque pour les fichiers temporaires si l'emprise de la région est importante...

avec orthophotos de la même zone

 

Tous les traitements ont été effectués sur Mac OS X avec Grass 6.4 RC5 de www.kyngchaos.com/ 

Site officiel : GRASS
Site officiel : listes Grass