<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
<title></title>
</head>
<body bgcolor="#ffffff" text="#000000">
I am trying to calculating 95th percentile within polygons from a
of set realizations - something like zonal statistics. <br>
How do I calculate <font face="Courier New, Courier, monospace"> 95
th percentile for each polygon over all realizations. </font><br>
Thanks<br>
Zia<br>
<br>
For example:<br>
<br>
<font face="Courier New, Courier, monospace">data(meuse)<br>
data(meuse.grid)<br>
coordinates(meuse) <- ~x+y<br>
coordinates(meuse.grid) <- ~x+y<br>
<br>
# Simulation<br>
nsim=10<br>
x <- krige(log(zinc)~1, meuse, meuse.grid, model = vgm(.59,
"Sph", 874, .04), nmax=10, nsim=nsim)<br>
over(sr, x[,1:4], fn = mean)<br>
<br>
> over(sr, x[,1:4], fn = mean)<br>
sim1 sim2 sim3 sim4<br>
r1 5.858169 5.792870 5.855246 5.868499<br>
r2 5.588570 5.452744 5.596648 5.516289<br>
r3 5.798087 5.860750 5.784194 5.848194<br>
r4 NA NA NA NA<br>
<br>
# 95 th percentile at prediction grid:<br>
x<-as.data.frame(x)<br>
y95<-apply(x[3:nsim],1,stats::quantile,probs = 0.95,na.rm=TRUE)
# 95 th percentile at each prediction grid <br>
<br>
</font><br>
<br>
On 10/18/2011 3:30 PM, Edzer Pebesma wrote:
<blockquote cite="mid:4E9DD3C7.2040706@uni-muenster.de" type="cite">require(sp)
<br>
r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007,
180409,
<br>
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
<br>
332618, 332413, 332349))
<br>
r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142,
179437,
<br>
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
<br>
331133, 331623, 332152, 332357, 332373))
<br>
r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329,
179875,
<br>
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
<br>
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
<br>
329783, 329665, 329720, 329933, 330478, 331062, 331086))
<br>
r4 = cbind(c(180304, 180403,179632,179420,180304),
<br>
c(332791, 333204, 333635, 333058, 332791))
<br>
<br>
sr1=Polygons(list(Polygon(r1)),"r1")
<br>
sr2=Polygons(list(Polygon(r2)),"r2")
<br>
sr3=Polygons(list(Polygon(r3)),"r3")
<br>
sr4=Polygons(list(Polygon(r4)),"r4")
<br>
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
<br>
srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2),
row.names=c("r1","r2","r3","r4")))
<br>
<br>
data(meuse)
<br>
coordinates(meuse) = ~x+y
<br>
<br>
plot(meuse)
<br>
polygon(r1)
<br>
polygon(r2)
<br>
polygon(r3)
<br>
polygon(r4)
<br>
# retrieve mean heavy metal concentrations per polygon:
<br>
# attribute means over each polygon, NA for empty
<br>
over(sr, meuse[,1:4], fn = mean)
<br>
<br>
# return the number of points in each polygon:
<br>
sapply(over(sr, geometry(meuse), returnList = TRUE), length)
</blockquote>
</body>
</html>