[R-sig-Geo] How to do Zonal Statistics after Kriging in R given Shapefile of Polygon

Jesse Berman berman.jesse at gmail.com
Thu Sep 12 17:44:15 CEST 2013

Edzer is correct about this.  I was just looking at spTransform and adjusted
the data into a Robinson global projection.  If you set muni.sp and
brgy_poly to this projection the block kriging now works.  

You'll have to adjust your variogram, as it becomes a singular model, but
for demonstration purposes it works.  As an aside, you may want to consider
numbers of point pairs when setting your variogram.  

robin<-"+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
+units=m +no_defs"
muni.sp2<-spTransform(muni.sp, CRS(robin))
brgy_poly2<-spTransform(brgy_poly, CRS(robin))

## Variogram plot
v.obj<-variogram(POVERTY_INCIDENCE~1, locations=coordinates(sp_point),
data=muni.sp2, cloud=F)

## Assuming exponential model
v.exp <- fit.variogram(v.obj,vgm(psill=1, model='Exp', range=1))

## Krige and Plot
kr <- krige(POVERTY_INCIDENCE~1, muni.sp2, newdata=brgy_poly2, model=v.exp)
spplot(kr, "var1.pred", col.regions = rev(topo.colors(20)))

Jesse D Berman, PhD
Yale University
School of Forestry and Environmental Studies
Post-Doc Fellow
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/How-to-do-Zonal-Statistics-after-Kriging-in-R-given-Shapefile-of-Polygon-tp7584330p7584597.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

More information about the R-sig-Geo mailing list