%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)

%%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']

%%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)

%%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")

%%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)

%%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)

%%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))

fig = mlab.surf(donnees, warp_scale=0.5) 
mlab.savefig(filename='test.png')