Re: [R-sig-Geo] How to do Zonal Statistics after Kriging in R given Shapefile of Polygon
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 In the script, you set the wrong CRS to the long/lat points. Set the right CRS, then spTransform to the new one. Check that they match after reprojection by plotting with axes, eg. plot(muni.sp,axes=T) points(muni.sp, col = 'red') Second, when you work in projected coordinates having large numbers, using 1 as an initial value for a variogram range fit will not work; choose sth else, e.g. 10 or so. On 09/13/2013 06:12 AM, James Matthew Miraflor wrote: Hi Jesse, Edzer, and R-SIG-GEO! I was able to make Kriging work, following the suggestions of Jesse and Edzer. However, the results don't turn out to be right - the Kriging and IDW done on a grid is not the same (in fact, it seems to be a mirror image) of the block kriging using a separate map. The grid kriging is in fact closer to the small area estimates (which is the actual sample). You may take a look at the images here: https://www.dropbox.com/sh/q5sev7n2pi6tq8f/tiqiFMjJS5 The R Code is here: https://www.dropbox.com/sh/q5sev7n2pi6tq8f/oM87f6axv5/kriging_pov.r The R Data is here: https://www.dropbox.com/sh/q5sev7n2pi6tq8f/8jC9qmlBZB/POVKRIGE.RData Would you know why block Kriging and block IDW have such a weird results? Thank you for your time and patience. James On Fri, Sep 13, 2013 at 12:26 AM, James Matthew Miraflor james.miraf...@gmail.com wrote: Thanks Jesse and Edzer! It finally worked. On Thu, Sep 12, 2013 at 11:44 PM, Jesse Berman berman.je...@gmail.comwrote: 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)) -- * * *James Matthew B. Miraflor* MS Computer Science College of Engineering University of the Philippines http://politicsforbreakfast.blogspot.com No problem can withstand the assault of sustained thinking. - Voltaire - -- Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251 83 33081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.11 (GNU/Linux) Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/ iQEcBAEBAgAGBQJSMszMAAoJEM1OCHCtOnfxFfEIAIOSncLp+X8cPfhxD3+uZm1g QU+FL2EKrDippniw0i8fGv9+QYIpUf0qxQzMXsDtQj3nqitGSU8J2Hqqlytgj6z2 0gwdlP2Wrv2PWzotH7BJEZVJbl24kchIAXocPwa+FFYBEO+mY8t+iDaGR3sZTHK1 /Nt7kFWN2dRxaCg4Iagv38V3kV5sY/G3X6GREv05DUQ8vdzAzyUzL6byi8rr4E4H dsSo8GegsD7SlqwqSofLAp+MmWlgcgC1l5I+npKF4XG8lUaai3ttVKT0Ulw+4YvU 7qAzlC/8xCOFnlI/Mhp7uVrtx72sl+4NHouS6FSSbF45ZoFtZnPQknLQWUUVklQ= =iAjf -END PGP SIGNATURE- ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] How to do Zonal Statistics after Kriging in R given Shapefile of Polygon
James, it seems like your projections for 'town' and 'muni.sp' are different with one being lat/long and the other being a UTM Zonal. Try setting them the same and see if that solves your problem. Jesse - 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-tp7584330p7584595.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] How to do Zonal Statistics after Kriging in R given Shapefile of Polygon
Thanks Jesse and Edzer! It finally worked. On Thu, Sep 12, 2013 at 11:44 PM, Jesse Berman berman.je...@gmail.comwrote: 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)) -- * * *James Matthew B. Miraflor* MS Computer Science College of Engineering University of the Philippines http://politicsforbreakfast.blogspot.com 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
Re: [R-sig-Geo] How to do Zonal Statistics after Kriging in R given Shapefile of Polygon
Hi Jesse, Edzer, and R-SIG-GEO! I was able to make Kriging work, following the suggestions of Jesse and Edzer. However, the results don't turn out to be right - the Kriging and IDW done on a grid is not the same (in fact, it seems to be a mirror image) of the block kriging using a separate map. The grid kriging is in fact closer to the small area estimates (which is the actual sample). You may take a look at the images here: https://www.dropbox.com/sh/q5sev7n2pi6tq8f/tiqiFMjJS5 The R Code is here: https://www.dropbox.com/sh/q5sev7n2pi6tq8f/oM87f6axv5/kriging_pov.r The R Data is here: https://www.dropbox.com/sh/q5sev7n2pi6tq8f/8jC9qmlBZB/POVKRIGE.RData Would you know why block Kriging and block IDW have such a weird results? Thank you for your time and patience. James On Fri, Sep 13, 2013 at 12:26 AM, James Matthew Miraflor james.miraf...@gmail.com wrote: Thanks Jesse and Edzer! It finally worked. On Thu, Sep 12, 2013 at 11:44 PM, Jesse Berman berman.je...@gmail.comwrote: 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)) -- * * *James Matthew B. Miraflor* MS Computer Science College of Engineering University of the Philippines http://politicsforbreakfast.blogspot.com No problem can withstand the assault of sustained thinking. - Voltaire -- * * *James Matthew B. Miraflor* MS Computer Science College of Engineering University of the Philippines http://politicsforbreakfast.blogspot.com 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
[R-sig-Geo] How to do Zonal Statistics after Kriging in R given Shapefile of Polygon
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,.01) sp_point[,2] - jitter(muni$LAT,.01) 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
