Hi!
I was trying to do some zonal statistics after kriging based on poverty
incidence, but I can't seem to get the hang of it. In the code below, I am
trying to estimate town level poverty incidence after I kriged the actual
observed poverty incidence at muni.csv. I have the shapefiles of the towns
at base//town.shp, and data at base//town.dbf. I am trying to get the
average poverty incidence in the area covered by the town polygons.
Would anyone know how? Thanks!
James
####################################################
muni <- read.delim("muni.csv", sep=",")
sp_point <- matrix(NA, nrow=nrow(muni),ncol=2)
sp_point[,1] <- jitter(muni$LON,.000001)
sp_point[,2] <- jitter(muni$LAT,.000001)
colnames(sp_point) <- c("LON","LAT")
## Create spatial object
muni.sp <-
SpatialPointsDataFrame(coords=sp_point,muni,proj4string=CRS("+proj=utm
+zone=51 +datum=WGS84"))
par(mar=rep(0,4))
plot(muni.sp,pch=1,cex=muni.sp$POVERTY_INCIDENCE)
## Variogram plot
v.obj<-variogram(POVERTY_INCIDENCE~1, locations=coordinates(sp_point),
data=muni.sp, cloud=F)
plot(v.obj,type='b',pch=16)
## Assuming exponential model
v.exp <- fit.variogram(v.obj,vgm(psill=1, model='Exp', range=1))
plot(v.obj, v.exp, pch = 16,cex=.5)
v.fin = v.exp
## Ordinary kriging
grd <- Sobj_SpatialGrid(muni.sp)$SG
plot(grd,axes=T,col="grey")
points(muni.sp)
kr <- krige(POVERTY_INCIDENCE~1, muni.sp, grd, model=v.fin)
idw.out <- idw(POVERTY_INCIDENCE~1,muni.sp,grd,idp=.2)
r = raster(idw.out)
town <- readShapePoly("base\\town.shp",IDvar="TOWN_CODE",
proj4string=CRS("+proj=longlat +ellps=clrk66"))
zonal(r, town, stat='mean')
####################################################
--
*
*
*James Matthew B. Miraflor*
MS Computer Science
College of Engineering
University of the Philippines
http://politicsforbreakfast.net
________
"No problem can withstand the assault of sustained thinking." - Voltaire
________
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
