[R-sig-Geo] Calculating 95th percentile within polygons

[Zia Ahmed]

  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) 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20111018/95cd4b4e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: zua3.vcf
Type: text/x-vcard
Size: 281 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20111018/95cd4b4e/attachment.vcf>