Utilisation de GRASS GIS avec Python, R et Octave

In [1]:
from IPython.display import Image
Image(url='http://www.portailsig.org/sites/default/files/zeropoint_logo.png')
Out[1]:

Définition des variable systèmes pour accéder à GRASS GIS 6.4.x à partir de Python

In [2]:
import os
os.environ['GISBASE'] = '/Applications/GRASS-6.4.app/Contents/MacOS'
os.environ['PYTHONPATH']="/Applications/GRASS-6.4.app/Contents/MacOS/etc/python/"
os.environ['LD_LIBRARY_PATH']="/Applications/GRASS-6.4.app/Contents/MacOS/lib"
os.environ['GIS_LOCK'] = '$$'
os.environ['GISRC'] = '/Users/martinlaloux/.grassrc6'
import grass.script as grass

Accès à un secteur (LOCATION) et à un jeu de données (MAPSET)

simple exemple tiré de GRASS GIS v. 6.4.2: la nouvelle console Python (Shell Python), exemple d'utilisation simple ou avec le module Shapely

In [3]:
gisbase = os.environ['GISBASE']
gisdb="/Users/martinlaloux/grassdata"
location="geol"
mapset="test_python"
import grass.script.setup as gsetup
gsetup.init(gisbase, gisdb, location, mapset)
Out[3]:
'/var/folders/vu/vu099NovH-KZpE5TEa83BU+++TI/-Tmp-/tmpoCupAJ'

Extraction des coordonnées des extrémités d'une ligne

In [4]:
debligne =grass.read_command("v.to.db", flags="p", map="testgrass", type="line", option="start", units="meters" , quiet=True)
In [5]:
debx=float(debligne.split("|")[1])
deby=float(debligne.split("|")[2])
debz=float(debligne.split("|")[3])
In [6]:
finligne=grass.read_command("v.to.db", flags="p", map="testgrass",type="line", option="end", units="meters", quiet=True)
In [7]:
finx=float(finligne.split("|")[1])
finy=float(finligne.split("|")[2])
finz=float(finligne.split("|")[3])
 

Transformation en géométrie Shapely

In [8]:
from shapely.geometry import Point
In [9]:
Point1 = Point(debx,deby)
In [10]:
Point2 = Point(finx,finy)
In [11]:
from shapely.geometry import LineString
In [12]:
line = LineString([Point1,Point2])
In [13]:
%matplotlib inline
In [14]:
print line
LINESTRING (186139.123704173 53082.6541939769, 188199.122798505 53467.7585587325)

Affichage de la ligne avec matplotlib

In [16]:
import matplotlib.pyplot as plt 
fig = plt.figure() 
ax = fig.gca() 
x, y = line.xy
ax.plot(x, y, color='b')
Out[16]:
[<matplotlib.lines.Line2D at 0x1068c4c90>]

Affichage de la ligne avec R

In [17]:
load_ext rpy2.ipython
In [18]:
%%R -w 480 -h 480 -u px
library(ggplot2)
dat <- data.frame(x = rnorm(2), y = rnorm(2), 
                  lab = sample(c('deb', 'fin'), 10, replace = TRUE))
x <- ggplot(dat, aes(x = x, y = y, color = lab)) + geom_point() + geom_line()
print(x)

Affichage de la ligne avec Octave (clone de Matlab)

In [19]:
%load_ext oct2py.ipython
In [20]:
%octave -i x x
[[ 186139.12370417  188199.12279851]]
Out[20]:
array([[ 186139.12370417,  188199.12279851]])
In [21]:
%octave -i y y
[[ 53082.65419398  53467.75855873]]
Out[21]:
array([[ 53082.65419398,  53467.75855873]])
In [22]:
%%octave -s 480,480 -f svg
plot(x,y);
Produced by GNUPLOT 4.2 patchlevel 6 53000 53100 53200 53300 53400 53500 186000 186500 187000 187500 188000 188500

Traitements divers avec R

In [23]:
%Rpush x y
In [24]:
%R lm(y~x)$coef
Out[24]:
<Matrix - Python:0x121de5c68 / R:0x122224808>
[53082.654194, NA_real_, NA_real_, 53467.758559, NA_real_, NA_real_]
In [25]:
Xr = x - x.mean(); Yr = y - y.mean()
pente = (Yr*Yr).sum() / (Xr**2).sum()
interception = x.mean() - x.mean() * pente
(interception, pente)
Out[25]:
(180627.92926224536, 0.034948039908857427)
In [26]:
%R resid(lm(y~x)); coef(lm(x~y))
Out[26]:
<Matrix - Python:0x107606b48 / R:0x10f824e60>
[186139.123704, NA_real_, NA_real_, 188199.122799, NA_real_, NA_real_]
In [27]:
b = %R a=resid(lm(y~x))
In [28]:
print b
  [,1] [,2]
1    0    0

In [29]:
v1 = %R plot(x,y); print(summary(lm(y~x))); vv=mean(x)*mean(y)
print 'v1 est égal à :', v1
v2 = %R mean(x)*mean(y)
print 'v2 est égal à:', v2
Response Y1 :

Call:
lm(formula = Y1 ~ x)

Residuals:
ALL 1 residuals are 0: no residual degrees of freedom!

Coefficients: (2 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    53083         NA      NA       NA
x1                NA         NA      NA       NA
x2                NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom


Response Y2 :

Call:
lm(formula = Y2 ~ x)

Residuals:
ALL 1 residuals are 0: no residual degrees of freedom!

Coefficients: (2 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    53468         NA      NA       NA
x1                NA         NA      NA       NA
x2                NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom


v1 est égal à: [1] 9971473668

v2 est égal à: [1] 9971473668