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

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

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

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

%%R
plot(predR[["var1.pred"]])
points(df, pch = 19, cex = 0.1)

%%R
plot(predR[["var1.stdev"]])

%%R
contour(predR[["var1.pred"]])