%matplotlib inline
import numpy as np
from pandas import DataFrame, Series
from geostatsmodels import utilities, kriging, variograms, model, visualization
import matplotlib.pyplot as plt
from scipy.stats import norm
import urllib2
import os.path
import zipfile
import StringIO
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(gstat)
library(sp)
%%R
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) )
%%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
require(automap)
Le chargement a nécessité le package : automap
%%R
kriging_result = autoKrige(z~1, df, grd)
[using ordinary kriging]
%%R
plot(kriging_result)
%%R
x.range[1]
[1] 252338
%%R
x.range[2]
[1] 252338
%%R
x1=seq(from=x.range[1], to=x.range[2], by=500)
%%R
x1
[1] 251215 251715 252215
%%R
df@coords[,1]
[1] 252290.2 251238.0 251845.8 251581.9 251918.4 251994.6 251647.7 252338.6 [9] 252149.2 251477.9 251328.6 251533.5 252166.1 251215.5 251717.0 251931.4 [17] 251761.0 252035.2 252194.9 251588.1 251622.6 251304.6 251447.5 251413.6 [25] 251235.0 251871.1 251879.5 251236.9 251484.2 251924.0 251503.5 251333.5 [33] 252294.2 251533.5 251556.0 251278.0 251826.5 251822.9 252279.1 252213.3 [41] 251405.1 251521.9 251290.3 252152.3 251675.4 251633.5 252197.8 252177.5 [49] 251296.6 251600.4 251265.1 252199.1 252198.7 252210.1 252326.3 251845.5 [57] 251713.0 252144.7 252188.7 251268.9 251357.0 251381.6 251308.1 251571.6 [65] 251293.1 251513.0 251377.8 251389.6 251359.1 251816.6 251404.1 251917.7 [73] 252263.8 251964.4 251650.8 251857.4 252228.6 251798.2 251909.3 251248.9 [81] 251570.5 251396.2 251963.5 251739.3 251266.3 251757.8 251993.8 251628.2 [89] 251819.3 252212.7 251487.1 252100.7 251672.3 252076.6 252301.3 251346.3 [97] 251809.0 251806.8 251956.3 252129.7
%%R
x.range = as.integer(range(df@coords[,1]))
x.range
[1] 251215 252338
%%R
summary(x.range)
Min. 1st Qu. Median Mean 3rd Qu. Max. 251200 251500 251800 251800 252100 252300
%%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) )
%%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)
[using ordinary kriging]
%%R
plot(kriging_result)
%%R
require(rgdal)
# Pull sp SpatialPixelsDataFrame
kriging.pred <- kriging_result$krige_output
# Write Kriging estimate and variance
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar.asc")
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
plot(kriging.pred["var1.pred"])
%%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
# Pull sp SpatialPixelsDataFrame
kriging.pred <- kriging_result$krige_output
# Write Kriging estimate and variance
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred10.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar10.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) )
%%R
coordinates(grd) = ~ x+y
gridded(grd) = TRUE
%%R
kriging_result = autoKrige(z~1, df, grd)
[using ordinary kriging]
%%R
plot(kriging_result)
%%R
# Pull sp SpatialPixelsDataFrame
kriging.pred <- kriging_result$krige_output
# Write Kriging estimate and variance
writeGDAL(kriging.pred["var1.pred"], "Krigpt3Pred5.asc")
writeGDAL(kriging.pred["var1.var"], "Krigpt3DVar5.asc")