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

[Edzer Pebesma]

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.pebesma at wwu.de