On 10/18/2011 11:45 PM, Zia Ahmed wrote:
Your are correct. Now, what I am going to do: first calculate 95th percentile of all realizations at prediction grid and then calculate the mean for each polygon.
this implies you end up with one single value per polygon. Call it A.
Is it possible to calculated weighted mean of the this 95th percentile? weight may be area of polygon or any value, for example population density for each polygon.
if this population density is in the polygon attributes, then this is trivial, as the number of records equals that of output A.
Thanks Zia On 10/18/2011 5:12 PM, Edzer Pebesma wrote:Zia, You may want to rethink the question. Each realization has a 95 percentile within a particular polygon. Over the set of realizations of some aggregated value for a polygon, you can take a 95 percentile. These are two different things. The first is a spatial aggregation, the second an aggregation over the (sampled) probability distribution. On 10/18/2011 11:07 PM, Zia Ahmed wrote:I am trying to calculating 95th percentile within polygons from a of set realizations - something like zonal statistics. How do I calculate 95 th percentile for each polygon over all realizations. Thanks Zia For example: data(meuse) data(meuse.grid) coordinates(meuse) <- ~x+y coordinates(meuse.grid) <- ~x+y # Simulation nsim=10 x <- krige(log(zinc)~1, meuse, meuse.grid, model = vgm(.59, "Sph", 874, .04), nmax=10, nsim=nsim) over(sr, x[,1:4], fn = mean) > over(sr, x[,1:4], fn = mean) sim1 sim2 sim3 sim4 r1 5.858169 5.792870 5.855246 5.868499 r2 5.588570 5.452744 5.596648 5.516289 r3 5.798087 5.860750 5.784194 5.848194 r4 NA NA NA NA # 95 th percentile at prediction grid: x<-as.data.frame(x) y95<-apply(x[3:nsim],1,stats::quantile,probs = 0.95,na.rm=TRUE) # 95 th percentile at each prediction grid On 10/18/2011 3:30 PM, Edzer Pebesma wrote:require(sp) r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 332618, 332413, 332349)) r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 331133, 331623, 332152, 332357, 332373)) r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110), c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 329783, 329665, 329720, 329933, 330478, 331062, 331086)) r4 = cbind(c(180304, 180403,179632,179420,180304), c(332791, 333204, 333635, 333058, 332791)) sr1=Polygons(list(Polygon(r1)),"r1") sr2=Polygons(list(Polygon(r2)),"r2") sr3=Polygons(list(Polygon(r3)),"r3") sr4=Polygons(list(Polygon(r4)),"r4") sr=SpatialPolygons(list(sr1,sr2,sr3,sr4)) srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2), row.names=c("r1","r2","r3","r4"))) data(meuse) coordinates(meuse) = ~x+y plot(meuse) polygon(r1) polygon(r2) polygon(r3) polygon(r4) # retrieve mean heavy metal concentrations per polygon: # attribute means over each polygon, NA for empty over(sr, meuse[,1:4], fn = mean) # return the number of points in each polygon: sapply(over(sr, geometry(meuse), returnList = TRUE), length)
-- Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/geostatistics e.pebe...@wwu.de _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
