Grass Gis, ses amis Python et R, ou la puissance d'un shell pas si "ringard" que ça... (traitement des données, géostatistiques and Co.)

On reproche, souvent à juste titre, à l'interface de Grass d'être beaucoup moins ergonomique ou intuitive que d'autres logiciels SIG, tels ArcGIS ou même Quantum Gis. 

Et pourtant, une belle interface graphique n'est souvent que l'arbre qui cache la forêt.  Ce qui fait la puissance de Grass est sa ligne de commande, son shell.  Une fois maitrisé, il permet, entre autres choses, d'enchainer les commandes à l'aide de "pipes" ( fr.wikipedia.org/wiki/Tube_Unix) ou de faire des scripts d'automatisation de traitements (la plupart des commandes de Grass sont des scripts).

Bien vite on se rend compte que l'utilisation de la ligne de commandes est beaucoup plus "souple" que l'interface graphique. Le parallèle avec ArcInfo, lignes de commandes et scripts AML,  est souvent préjudiciable pour Grass, car cette approche est vue maintenant comme dépassée, voire archaïque comparée aux beaux habits actuels des SIG. Et pourtant ESRI indique encore que  "ArcInfo is the most complete desktop GIS" (www.esri.com/software/arcgis/arcinfo/index.html), donc l'approche ne doit pas être si dépassée que ça.

Le shell Grass offre une possibilité unique,  celle d'utiliser d'autres programmes pour effectuer les traitements des données Grass et de les récupérer par la suite, et tout ceci, sans le quitter (geoinformatics.fsv.cvut.cz/gwi/On_scripting_GRASS_GIS:_Building_location independent_command_line_tools).

Parmi ceux-ci, il y a Python et R mais comment cela fonctionne-t-il ? Contrairement au billet précédent sur Grass où je n'avais pas abordé l'aspect codes, ici, shell oblige, il faudra mettre ses mains dans le cambouis.

Voici une commande classique du shell qui permet de connaitre les caractéristiques d'un fichier vecteur:

GRASS 6.4.0RC5 (geol):~ > v.info map=newcondZ@moi layer=1
+----------------------------------------------------------------------------+
| Layer:           newcondZ@moi                                              |
| Mapset:          moi                                                       |
| Location:        geol                                                      |
|[....]                                                                      |
|----------------------------------------------------------------------------|
|   Type of Map:  vector (level: 2)                                          |
|                                                                            |
|   Number of points:       533             Number of areas:      0          |
|   Number of lines:        0               Number of islands:    0          |
|   Number of boundaries:   0               Number of faces:      0          |
|   Number of centroids:    0               Number of kernels:    0          |
|                                                                            |
|   Map is 3D:              Yes                                              |
|   Number of dblinks:      1                                                |
|                                                                            |
|         Projection: x,y                                                    |
|               N:            117618    S:         113774.47                 |
|               E:            176201    W:         165638.91                 |
|               B:            -64.89    T:               268                 |
|[....]                                                                      |
+----------------------------------------------------------------------------+
GRASS 6.4.0RC5 (geol):~ > 

Grass - Python, main dans la main

La suite est alors très simple, dans l'invite de commande, on lance le programme désiré, Python ici:

GRASS 6.4.0RC5 (geol):~ > python
Python 2.5.4 
[....]
>>>

et on se retrouve dans Python, mais à l'intérieur de Grass. Comment alors l'utiliser avec les données de Grass puisque, à priori, il n'y a pas de communication ? En pratique, il est possible d'utiliser Python de deux manières (grass.osgeo.org/wiki/GRASS_and_Python):

• soit comme simple interface (Swig) aux commandes de Grass download.osgeo.org/grass/grass6_progman/swig/ - grass.osgeo.org/wiki/PythonSwigExamples
• soit en pur Python avec un module/librairie ( download.osgeo.org/grass/grass6_progman/pythonlib.html)

En voici quelques exemples:

- avec des vecteurs

>>> from grass.script import vector as grassv
>>> grassv.vector_db('newcondZ')
{1:
{'layer': '1', 'name': '', 'database':
'/Users/toto/grassdata/geol/moi/dbf/', 'driver': 'dbf', 'key':
'cat', 'table': 'newcondZ'}}
>>> grassv.vector_info_topo('newcondZ')
{'kernels':
0, 'lines': 0, 'centroids': 0, 'boundaries': 0, 'points': 533, 'faces':
0, 'primitives': 533, 'islands': 0, 'nodes': 533, 'map3d': 1, 'areas':
0}
>>> grassv.vector_columns('newcondZ')
{'int_3':
'INTEGER', 'int_2': 'INTEGER', 'int_1': 'INTEGER', 'cat': 'INTEGER',
'dbl_2': 'DOUBLE PRECISION', 'dbl_3': 'DOUBLE PRECISION', 'dbl_1':
'DOUBLE PRECISION', 'dbl_4': 'DOUBLE PRECISION'}
>>> grass.vector_db_select('newcondZ')['values'][2]
['2', '169097.19', '113891.51', '189.76', '192', '80', '1', '0.5']
>>> 

- avec des rasters

>>> from grass.script import raster as grassr
>>> grassr.raster_info('s_newfrst3@moi')
{'north':
118869.900569, 'timestamp': '"none"', 'min': -456.08587646484398,
'datatype': 'FCELL', 'max': 265.500732421875, 'ewres':
14.008076920000001, 'vertical_datum': '', 'west': 164160.653425,
'units': '', 'title': ' (s_newfrst3)', 'east': 176641.849961, 'nsres':
13.976472729999999, 'south': 112552.534895}
>>> 
# on quitte Python et on se retrouve dans la shell Grass 
GRASS 6.4.0RC5 (geol):~ >

Tous les traitements sont donc possibles avec des modules comme Numpy, Scipy, Matplotlib ou Shapely. Les résultats peuvent ensuite être renvoyés dans Grass.

Traitement avec Grass puis traitement avec Scipy et Matplotlib avec les données Grass:

traitement en Grass avec
r.contour input=s_newfrst3 output=s_newfrst_contours minlevel=-230 maxlevel=200 step=50 cut=20

même traitement (ici, traitement brut, sans projection et sans définir d'intervalle, contrairement à r.contour) avec Scipy  et Matplotlib (à droite le shell Grass)

Tous ces traitements en Python seront de plus en plus développés dans les nouvelles versions de Grass comme la 7 qui est en chantier (où il sera possible de piloter tout Grass en Python comme le démontrent la plupart des commandes classiques, voir trac.osgeo.org/grass/browser/grass/trunk/scripts). Il est déjà possible de développer des scripts comme ceux de grass.osgeo.org/wiki/GRASS_and_Python#Interfacing_with_NumPy où Numpy est utilisé pour faire des traitements propres à ce module et où les résultats sont renvoyés dans Grass. La gamme des possibilités de Grass se voit donc augmentée considérablement.

Si Python ne doit plus être présenté, R oui. Il s'agit d'un environnement mathématique et d'un langage (Open Source) pour le traitement et l'analyse statistique des données et leur visualisation (3D comprise). Il est très utilisé dans le monde scientifique (et autre) et moins connu dans le monde des Sigistes (il suffit d'une simple recherche sur Google (vu la difficulté de chercher "R" seul, il y a un moteur de recherche personnalisé Google  www.rseek.org/). Récemment un article lui a été consacré dans le New York Times "Data Analysts Captivated by R's Power" (www.nytimes.com/2009/01/07/technology/business-computing/07program.html), et il a suscité de nombreuses réactions d'utilisateurs, célèbres ou non  (bits.blogs.nytimes.com/2009/01/08/r-you-ready-for-r/).

Présentation très sommaire de R

R est un monde en soi, gigantesque, et dont l'application dans le domaine géospatial ne représente qu'une petite partie (fr.wikipedia.org/wiki/R_%28logiciel%29). Dans sa version de base, il dispose de tout ce qui est nécessaire pour les traitements statistiques classiques. Tout comme Python, ses possibilités se voient considérablement accrues par l'adjonction de très nombreuses librairies ("packages" ou extensions, quasi 2000). Ces librairies permettent à R de pratiquement tout traiter depuis la statistique multivariée jusqu'à l'analyse complète des données écologiques en passant par la biométrie, l'analyse des séries chronologiques, l'analyse bayésienne ou l'économétrie. Tous les domaines d'applications sont résumés sous forme de thèmes dans cran.r-project.org/web/views/.

   exemple de traitement de données sismiques, refait d'après learnr.wordpress.com/2009/07/20/ggplot2-version-of-figures-in-lattice-multivariate-data-visualization-with-r-part-6/

R possède bien évidemment des librairies pour traiter les données géospatiales (cran.r-project.org/web/views/Spatial.html,r-spatial.sourceforge.net/, geodacenter.asu.edu/r-spatial-projects/). Elles offrent, entre autres, la possibilité d'effectuer les traitements géostatistiques où R est très utilisé dans le monde de la recherche. Par leur entremise, R sait traiter les données SIG mais aussi interagir directement avec les programmes. En sus de ceux Open Source comme Grass, il peut aussi être couplé à ArcGIS (package RPyGeo, voir aussi www.forumsig.org/showthread.php ) ou à ArcInfo (RArcInfo). R sait aussi extraire des images de Google Maps pour les utiliser dans les traitements (RgoogleMaps). Tout cela indique que R est fort utilisé dans le domaine. Par exemple, Kamran Safi explique sur son site sites.google.com/site/spatialr/  comment il a été conduit à utiliser R en lieu et place d'ArcGis suite à son abandon par l'institution où il travaillait. Il y a cependant un problème pour l'utilisateur débutant. Chaque librairie s'est en effet développée parallèlement, avec ses formats et structures propres, et il est nécessaire d'apprendre à jongler entre ces éléments. Ce manque d'homogénéité est bien explicité par S.Dray sous le titre "Le grand bazar" dans biom3.univ-lyon1.fr/R/pdf/ter5.pdf. Tout ce qu'il est possible de faire est synthétisé dans le cours geography.uoregon.edu/GeogR/topics/index.html

R sait aussi se connecter à toutes les bases de données depuis Oracle jusqu'à Microsoft Access en passant par SQLite ou PostgreSQL/Postgis.

"Foin d'une interface graphique élaborée !",  tout se fait dans le langage de R. Plusieurs livres ont été écrits à son sujet (www.r-project.org/doc/bib/R-books.html , avec, notamment pour l'analyse géospatiale, "Applied Spatial Data Analysis with R" www.springer.com/public+health/epidemiology/book/978-0-387-78170-9). Le chapitre 10 du livre "Open Source GIS: A GRASS GIS Approach" (www.grassbook.org/) est consacré à ce sujet. Des sites comme bm2.genes.nig.ac.jp/RGM2/index.php ou addictedtor.free.fr/graphiques/ se proposent de montrer toutes ses facettes sous forme visuelle (plus de 21000 graphiques pour le premier cité), accompagnés des commandes permettant de les réaliser.  Ses domaines d'application sont donc extrêmement variés et ce n'est pas la maigre petite application qui suit qui lui rendra justice...

Grass - R, main dans la main

C'est le même principe qu'avec Python, il faut d'abord charger le package spgrass6 pour communiquer avec Grass. Le premier traitement présenté ici est une adaptation de celui de grass.osgeo.org/wiki/R_statistics  avec mes données.

GRASS 6.4.0RC5 (geol):~ > R
R version 2.9.2 (2009-08-24)
Copyright (C) 2009 The R Foundation for Statistical Computing
[....]
> library(spgrass6) # chargement du package spgrass6 qui permet d'interagir "en direct" avec Grass.
[....]
# lecture d'un raster et d'un MNT provenant de Grass dans la variable montrait
> G <- gmeta6() # paramètres de Grass, région etc.
> montrait <- readRAST6(c("geologie", "mnt"),cat=c(TRUE, FALSE), ignore.stderr=TRUE,plugin=NULL)
 
# puisque Grass utilise des rasters "intelligents" (catégories), il y a moyen d'obtenir le nombre de cellules dans chaque catégories (cat) qui sont ici les formations géologiques
> table(montrait$geologie)
 ACO  AIFA   ANS   BAU  ....
20116  2067  4485 12733 ....
 
# on peut visualiser dans R les rasters de Grass
> image(montrait, attr = 2, col = terrain.colors(20))
 
# l'auteur s'amuse ensuite à créer un diagramme "en boîte" (box plot) de la répartition des formations géologiques en fonction de l'altitude
> boxplot(montrait$mnt ~ montrait$geologie, medlwd = 1)
 
#création d'un nouveau MNT 
> montrait$nvmnt <- montrait$mnt * 2
 
# renvoi du nouveau MNT dans Grass
> writeRAST6(montrait, "mnt2", zcol="nvmnt", ignore.stderr=TRUE)
 
# on quitte R
> q()
GRASS 6.4.0RC5 (geol):~ > r.info map=mnt2@moi 
 +----------------------------------------------------------------------------+
 | Layer:    mnt2@moi                       Date: Mon Mar 22 22:12:59 2010    |
 | Mapset:   moi                            Login of Creator: yo              |
 | Location: geol                                                             |
 |[....]                                                                      |
 |   Type of Map:  raster               Number of Categories: 255             |
 |   Data Type:    FCELL                                                      |
 |   Rows:         191                                                        |
 |   Columns:      270                                                        |
 |   Total Cells:  51570                                                      |
 |        Projection: Lambert Conformal Conic                                 |
 |            N: 118892.8197047    S: 107481.49881847   Res: 59.74513553      |
 |            E: 178100.24901493    W: 161969.06242183   Res: 59.74513553     |
 [....]
GRASS 6.4.0RC5 (geol):~ >

             MNT vu dans Grass                                                                                 MNT vu dans R

diagramme en boîte ("boxplot", "boîte à moustache") de la répartition des formations géologiques en fonction de l'altitude. La taille de la boîte représente la surface sur la carte d'une formation

Nous allons utiliser le principe précédent pour faire une interpolation 3D par krigeage ordinaire à l'aide d'un des packages, gstat, en appliquant les traitements proposés par casoilresource.lawr.ucdavis.edu/drupal/node/438 ou blog.andrewl.net/articles/2008/05/02/digging-kriging.

C'est un exemple très simple qui ne se veut nullement une introduction à la géostatistique ou autre chose et qui se limite à un aspect purement exemplatif de la méthode sans aller plus loin (le krigeage suppose que certaines hypothèses soient validées). Néanmoins, en simplifiant, le traitement se fait en deux étapes, d'abord l'estimation d'un variogramme (on parle aussi de semi-variogramme) adapté aux données étudiées, puis l'interpolation par krigeage en fonction du variogramme choisi. Ce sera fait ici avec les paramètres par défaut pour un variogramme non directionnel. Pour corser le tout, et bien que les données, le traitement et le modèle choisis soient loin d'être les plus appropriés, les points traités seront ceux utilisés dans www.portailsig.org/content/grass-gis-et-paraview-la-modelisation-geologique-3d-pour-les-demunis-ou-les-enthousiastes. Encore une fois, seule la démarche compte ici. Tous ces éléments sont amplement développés dans le livre "A Practical Guide to Geostatistical Mapping" qui utilise R en combinaison avec Saga Gis Ilwis, Grass voire Google Earth (spatial-analyst.net/book/, www.nceas.ucsb.edu/files/scicomp/GISSeminar/UseCases/CreateKMLFiles/ConvertShapefileToKML.html). Il peut être librement téléchargé à www.lulu.com/product/download/a-practical-guide-to-geostatistical-mapping/6379057.

Points vus dans Grass et dans R

Première étape: variogramme

# chargement de la librairie gstat
> library(gstat)
# lecture des points à traiter:
> mespoints = readVECT6('newcondZ@moi')
> coor <-coordinates(mespoints)
> coor
       coords.x   coords.y   coords.z 
  [1,]  169099.0  113900.0    239.00
  [2,]  169098.1  113895.8    214.38
  [3,]  169097.2  113891.5    189.76
  [...]
> z <- coor[,3]
> z 
 [1] 239.00 214.38 189.76 [...]
 
# l'interpolation nécessite une grille qui sera remplie par les valeurs interpolées
> grd <- gmeta2grd()
> grille <- SpatialGridDataFrame(grd, data=data.frame(k=rep(1,G$cols*G$rows)), proj4string=CRS(mespoints@proj4string@projargs))
 
# initialisation d'un objet gstat avec la valeur z à traiter 
> g <- gstat(id="elev", formula=z ~ 1, data=mespoints)
 
# vue des divers variogrammes (en fonctions de divers paramètres)
> plot(variogram(g,[...])

grille générée et points analysés

variogrammes directionnels obtenus en fonction des paramètres choisis.

Deuxième étape: interpolation par krigeage ordinaire (en fonction d'un variogramme linéaire simple retenu ici bien qu'il n'est pas le plus adapté)

# ajustement à un modèle de variogramme linéaire simple en utilisant les paramètres par défaut
> v.fit <- fit.variogram(variogram(g) ,vgm(model="Lin") )
 
# adaptation de l'objet gstat en fonction du modèle choisi  
> g <- gstat(g, id="elev", model=v.fit )
# interpolation par krigeage ordinaire
> krig_OK <- predict(g, id='elev', newdata=grille)
 
# renvoi des résultats de l'interpolation dans Grass
> writeRAST6(krig_OK, 'elev.ok', zcol='elev.pred')
# carte de variance d'estimation minimale du krigeage ou variance de krigeage
> writeRAST6(krig_OK, 'elev.ok_var', zcol='elev.var')
> q()
GRASS 6.4.0RC5 (geol):~ >

résultats graphiques dans R de l'interpolation avec courbes de niveau et points originels et variance de krigeage

interpolation et variance de krigeage en 3D avec le package rgl

Contrôlons le résultat dans Grass

GRASS 6.4.0RC5 (geol):~ > r.info map=elev.ok@moi 
 +----------------------------------------------------------------------------+
 | Layer:    elev.ok@moi                    Date: Sat Mar 20 13:11:01 2010    |
 | Mapset:   moi                            Login of Creator: yo              |
 | Location: geol                                                             |
 |[....]                                                                      |
 |   Type of Map:  raster               Number of Categories: 255             |
 |   Data Type:    FCELL                                                      |
 |   Rows:         168                                                        |
 |   Columns:      269                                                        |
 |   Total Cells:  45192                                                      |
 |        Projection: Lambert Conformal Conic                                 |
 |            N: 118202.17276904    S:  108164.99   Res: 59.74513553          |
 |            E: 178054.96845757    W: 161983.527   Res: 59.74513553          |
 |   Range of data:    min = -148.065659  max = 262.488892                    |
 [....]                                                                       
 
GRASS 6.4.0RC5 (geol):~ > 

 Comparaison du résultat obtenu avec l'interpolation faite avec v.surf.rst:

      surface obtenue avec v.surf.rst                                                    surface obtenue avec R et le krigeage simple

Visualisation avec Nviz (exagération verticale x3). Il faut rappeler, dans le cas présent, que le traitement avec v.surf.rst a été optimisé pour que la surface obtenue "colle" véritablement aux points traités alors que le krigeage est brut. La surface doit encore être "nettoyée" avec r.mapcalc pour obtenir celle limitée aux points.

Conclusions sur les bons amis et perspectives

On le voit les possibilités sont immenses, il suffit de consulter la table des matières de l'ouvrage "Applied Spatial Data Analysis with R" (340 pages)  déjà cité, et ce, tant avec R qu'avec Python.

Et la boucle est bouclée, car la souplesse du shell de Grass et de ses scripts permet d'intégrer tous ces traitements, s'ils sont intéressants, comme nouveaux modules standards. C'est le cas avec le processus Grass-R que nous venons de faire et qui a été intégré dans un futur module standard v.krige qui sera disponible dès la version 6.5. Il peut déjà être utilisé dès la version 6.4 RC5 (grass.osgeo.org/wiki/V.krige_GSoC_2009l), mais il faut l'installer soi-même. Il est nécessaire, bien entendu, que R soit aussi installé.

v.krige, image tirée de grass.osgeo.org/wiki/V.krige_GSoC_2009

Pour de plus amples renseignements et exemples sur l'emploi de Grass et de R, voir le site casoilresource.lawr.ucdavis.edu/drupal/node/100 (menu de gauche pour les détails), mais il ya beaucoup, beaucoup de sites qui montrent l'utilisation des packages spatiaux, car tous les objets SIG peuvent être traités, tant avec R qu'avec Python de manière indépendante grâce aux modules (GDAL/OGR, rgdal, ...), mais le couplage et l'interactivité seront perdus.

R peut aussi être utilisé à partir de Quantum Gis à l'aide d'extensions en Python. Il a déjà été utilisé avec Python sur le Portail dans www.portailsig.org/content/installation-d-un-serveur-wps-sous-windows-avec-pywps-et-r, (je confirme, pour mon plaisir inutile, que R est accessible à partir de Python dans Grass grâce au module Rpy2 , R dans Python dans Grass...).

En espérant que ce tout petit aperçu des possibilités de R dans le domaine spatial vous donnera des envies d'aller y faire un petit tour... 

Pourtant la majorité des utilisateurs passeront sans doute à côté, préférant les beaux habits...

"However, for me personally and for most users, both individual and organizational, the much greater cost of software is the time it takes to install it, maintain it, learn it and document it. On that, R is an epic fail. It does NOT fit with the way the vast majority of people in the world use computers. The vast majority of people are NOT programmers. They are used to looking at things and clicking on things" (www.thejuliagroup.com/blog/)

Opinion qui a évidemment valu de nombreux retours de flammes...

Image tirée de www.drewconway.com/zia/ et  [R] "is for making new things. Point and click is for redoing old things. I often need to make new things to analyze my data."(twitter.com/fdaapproved/status/12200318208)

Tous les traitements ont été effectués sur Mac OS X avec Grass 6.4 RC5 de www.kyngchaos.com/ , Python 2.5.4, et R version 2.9, paquet gstat version 0.9-63