(C Gumiaux,D Gapais,J.P Brun, 2003) pour des mesures de schistosité.
import fiona
import numpy as np
import matplotlib.pyplot as plt
c = fiona.open('/Users/Shared/telechargement/30_06_14/schisto_select.shp')
x,y = zip(*[i['geometry']['coordinates'] for i in c])
dir = [i['properties'][u'double_dir'] for i in c]
load_ext rpy2.ipython
x = np.array(x)
y = np.array(y)
dir = np.array(dir)
%Rpush x y dir
%%R
library(sp)
df = data.frame(x, y, dir)
coordinates(df) = ~x+y
%%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) )
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
plot(grd, cex=0.5)
points(df, pch=1, col='red', cex=0.7)
title("grille d'interpolation et points")
%%R
library(automap)
variogram = autofitVariogram(dir~1,df)
plot(variogram)
%%R
kriging_result = autoKrige(dir~1, df, grd)
plot(kriging_result, sp.layout = list(pts = list("sp.points",df)))
[using ordinary kriging]
%%R
require(rgdal)
# SpatialPixelsDataFrame
kriging.pred = kriging_result$krige_output
# Write Kriging estimate and variance
writeGDAL(kriging.pred["var1.pred"], "schisto_double_pred_500.tif")
writeGDAL(kriging.pred["var1.var"], "schisto_double_var_500.tif")
# points de la grille interpolée
writeOGR(kriging.pred, ".", "grille_500M", driver="ESRI Shapefile")
Le chargement a nécessité le package : rgdal 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)
%%R
grd2 = expand.grid(x=seq(from=x.range[1], to=x.range[2], by=1000), y=seq(from=y.range[1], to=y.range[2], by=1000) )
coordinates(grd2) = ~ x+y
gridded(grd2) = TRUE
plot(grd2, cex=0.5)
points(df, pch=1, col='red', cex=0.7)
title("grille d'interpolation et points")
%%R
kriging_result2 = autoKrige(dir~1, df, grd2)
plot(kriging_result2, sp.layout = list(pts = list("sp.points",df)))
[using ordinary kriging]
%%R
# SpatialPixelsDataFrame
kriging.pred2 = kriging_result2$krige_output
# Write Kriging estimate and variance
writeGDAL(kriging.pred2["var1.pred"], "schisto_double_pred_1000.tif")
writeGDAL(kriging.pred2["var1.var"], "schisto_double_var_1000.tif")
%%R
require(raster)
predR = stack(kriging_result$krige_output)
Le chargement a nécessité le package : raster
%%R
plot(predR[["var1.pred"]])
points(df, pch = 19, cex = 0.1)
%%R
plot(predR[["var1.stdev"]])
%%R
contour(predR[["var1.pred"]])
Cette démarche ne marche pas à tous les coups