%matplotlib inline
import numpy as np
import fiona
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]
load_ext rpy2.ipython
x = np.array(x)
y = np.array(y)
z = np.array(z)
%Rpush x y z
%%R
df = data.frame(x, y, z)
names(df)
[1] "x" "y" "z"
%%R
library(sp)
coordinates(df) = ~x+y
%%R
library(automap)
variogramme = autofitVariogram(z~1, df)
plot(variogram)
%%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) )
%%R
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
%%R
plot(grd, cex=0.5)
points(df, pch=1, col='red', cex=0.7)
title("grille d'interpolation et points")
%%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
%%R
plot(kriging_result)
%%R
automapPlot(kriging_result$krige_output, "var1.pred", sp.layout = list("sp.points", df))
%%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
%%R
plot(grd, cex=0.5)
points(df, pch=1, col='red', cex=0.7)
title("grille d'interpolation et points")
%%R
kriging_result = autoKrige(z~1, df, grd)
[using ordinary kriging]
%%R
plot(kriging_result)
%%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)
%%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) )
%%R
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
%%R
kriging_result = autoKrige(z~1, df, grd)
[using ordinary kriging]
%%R
plot(kriging_result)
%%R
plot(kriging_result, sp.layout = list(pts = list("sp.points", df)))
%%R
library(rgdal)
kriging.pred <- kriging_result$krige_output
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred2.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar2.asc")
%%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]
%%R
library(rgdal)
kriging.pred <- kriging_result$krige_output
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred3.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar3.asc")
from osgeo import gdal
ds = gdal.Open('/Users/martinlaloux/Krigpt3Pred3.asc')
donnees = ds.ReadAsArray()
from mayavi import mlab
mlab.figure(size=(400, 480), bgcolor=(0.16, 0.28, 0.46))
<mayavi.core.scene.Scene at 0x140ce25f0>
fig = mlab.surf(donnees, warp_scale=0.5)
mlab.savefig(filename='test.png')