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

[Zia Ahmed]

  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.  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.
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)
>
-------------- 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/1938ecdb/attachment.vcf>