Utilisation de Python (Fiona, Numpy et Mayavi) et de R pour l'analyse variographique et le krigeage

In [1]:
%matplotlib inline
In [2]:
import numpy as np
import fiona

Ouverture du shapefile avec Fiona

In [3]:
c = fiona.open('/Users/martinlaloux/rs3D/rs3D.shp')
x,y = zip(*[i['geometry']['coordinates']  for i in c])
z = [i['properties'][u'elev'] for i in c]

Activation de rmagic

In [4]:
load_ext rpy2.ipython

Transfert des données à R

In [5]:
x = np.array(x)
y = np.array(y)
z = np.array(z)
%Rpush x y z

création d'un dataframe spatial

In [6]:
%%R
df = data.frame(x, y, z) 
names(df)
[1] "x" "y" "z"
In [7]:
%%R
library(sp)
coordinates(df) = ~x+y

analyse variographique automatique avec automap

In [8]:
%%R
library(automap)
variogramme = autofitVariogram(z~1, df)
plot(variogram)

création d'une grille d'interpolation

In [9]:
%%R
x.range = as.integer(range(df@coords[,1]))
y.range = as.integer(range(df@coords[,2]))
grd = expand.grid(x=seq(from=x.range[1], to=x.range[2], by=500), y=seq(from=y.range[1], to=y.range[2], by=500) )
In [10]:
%%R
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
In [11]:
%%R
plot(grd, cex=0.5)
points(df, pch=1, col='red', cex=0.7)
title("grille d'interpolation et  points")

Krigeage avec automap

In [12]:
%%R
kriging_result = autoKrige(z~1, df, grd)
kriging_result['var_model']
[using ordinary kriging]
$var_model
  model    psill   range kappa
1   Nug   0.0000   0.000   0.0
2   Ste 899.7016 534.094   1.5

In [13]:
%%R
plot(kriging_result)
In [14]:
%%R
automapPlot(kriging_result$krige_output, "var1.pred", sp.layout = list("sp.points", df))

nouvelle grille au krigeage (espacement des points différent)

In [15]:
%%R
grd = expand.grid(x=seq(from=x.range[1], to=x.range[2], by=50), y=seq(from=y.range[1], to=y.range[2], by=50)) 
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
In [16]:
%%R
plot(grd, cex=0.5)
points(df, pch=1, col='red', cex=0.7)
title("grille d'interpolation et  points")
In [17]:
%%R
kriging_result = autoKrige(z~1, df, grd)
[using ordinary kriging]
In [18]:
%%R
plot(kriging_result)

Exportation des résultats

In [19]:
%%R
library(rgdal)
kriging.pred <- kriging_result$krige_output
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar.asc")
rgdal: version: 0.8-16, (SVN revision 498)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/24
Path to GDAL shared files: /Library/Frameworks/GDAL.framework/Versions/1.10/Resources/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: (autodetected)

avec QGIS

Autres tentatives en changeant l'espacement de la grille (pour visualiser son effet)

In [20]:
%%R
grd = expand.grid(x=seq(from=x.range[1], to=x.range[2], by=10), y=seq(from=y.range[1], to=y.range[2], by=10) )
In [21]:
%%R
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
In [22]:
%%R
kriging_result = autoKrige(z~1, df, grd)
[using ordinary kriging]
In [23]:
%%R
plot(kriging_result)
In [24]:
%%R
plot(kriging_result, sp.layout = list(pts = list("sp.points", df)))
In [25]:
%%R  
library(rgdal)
kriging.pred <- kriging_result$krige_output
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred2.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar2.asc")

avec QGIS

In [26]:
%%R
grd = expand.grid(x=seq(from=x.range[1], to=x.range[2], by=5), y=seq(from=y.range[1], to=y.range[2], by=5))
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
kriging_result = autoKrige(z~1, df, grd)
plot(kriging_result)
[using ordinary kriging]
In [27]:
%%R  
library(rgdal)
kriging.pred <- kriging_result$krige_output
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred3.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar3.asc")

avec QGIS

Retour à Python et visualisation 3D de la surface obtenue avec le module Mayavi

In [28]:
from osgeo import gdal
ds = gdal.Open('/Users/martinlaloux/Krigpt3Pred3.asc')
donnees = ds.ReadAsArray() 
In [29]:
from mayavi import mlab  
In [30]:
mlab.figure(size=(400, 480), bgcolor=(0.16, 0.28, 0.46))
Out[30]:
<mayavi.core.scene.Scene at 0x140ce25f0>
In [31]:
fig = mlab.surf(donnees, warp_scale=0.5) 
mlab.savefig(filename='test.png')

avec Mayavi