[R-sig-Geo] Counting points within polygons

[Agustin Lobo]

You are right the vignette is pretty clear.
I think this is a more elegant solution (BorneoLC is the Sp Pol DF and
fireLC2 the Sp Points DF):

> delme6 <- over(BorneoLC, fireLC2, returnList = TRUE)
> nbp <- sapply(delme6,nrow)
> newdata <- data.frame(BorneoLC at data,nbp=nbp)
> BorneoLC2 <- BorneoLC
> BorneoLC2 at data <- newdata
> BorneoLC2 at data[c(1,1000),]
    LC   ID LC_2         LCname LCabrv     AREA PERIMETER nbp
0    4 4505    4 Lowland Forest    LwF 0.000077  0.036385   0
999  8 6745    8 Lowland Mosaic    LwM 0.007168  0.864000   2

Perhaps having an sp function for the specific task of counting points
within polygons would make sense, I think that
would be quite popular.

Agus


2011/10/18 Edzer Pebesma <edzer.pebesma at uni-muenster.de>:
>
>
> On 10/18/2011 06:32 PM, Agustin Lobo wrote:
>>
>> I'm afraid I do not understand well the help page.
>> It seems to me that what I need is
>> "x = "SpatialPointsDataFrame", y = "SpatialPolygons"
>>  equal to the previous method, except that an argument fn=xxx is
>> allowed, e.g. fn = mean which will then report a data.frame with the
>> mean values of the x points falling in each polygon (set) of y"
>
> No, you need the other way around -- this is what I got wrong with overlay,
> long ago, where argument order did not matter. over(a,b) is the opposite of
> over(b,a)
>
> The text of the "over" documentation is hard to read, and I find the concept
> hard to explain in words, without giving examples. That's why I wrote the
> vignette.
>
>> Thus:
>>>
>>> delme4<- over(fireLC2[,5], BorneoLC,fn=sum)
>>
>> where
>>>
>>> BorneoLC at data[1,]
>>
>>   LC   ID LC_2         LCname LCabrv    AREA PERIMETER
>> 0  4 4505    4 Lowland Forest    LwF 7.7e-05  0.036385
>> and
>>>
>>> fireLC2 at data[1,]
>>
>>   LC ID                      LCname LCabrv nbp
>> 1  7  1 Plantation_Secondary Forest  PltSF   1
>>
>> but I get
>> Error in FUN(X[[1L]], ...) :
>>   only defined on a data frame with all numeric variables
>> Calls: over ... do.call ->  lapply ->  Summary.data.frame ->  lapply ->
>>  FUN
>>
> It would be helpful if you could provide a reproducible script that leads to
> this error, or even the command that gave you this error.
>
>> which I do not understand because fireLC2[,5] is numeric
>> (suggestion: could the user define the field(s) for which he wants fn
>> being applied? )
>
> How about passing an object with numeric attributes (selected)?
>
>> Beyond that, where the help page says:
>> "...which will then report a data.frame with the mean values of the x
>> points falling in each polygon (set) of y",
>> what happens if there is no points falling for a given polygon? Would
>> the record be eliminated or kept as NA ?
>>>
>>> From my point of view, having the option of getting a data.frame with
>>
>> the same nb. of rows and same ordering as
>> the input Sp Polygons object would be an advantage.
>
> That is why over was written to replace overlay.
>
> I modified the example in "over" to illustrate your two cases:
>
> ############ start script
> 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)
> ######## end script
>
> Get some grip on "over". Forget overlay, it will get deprecated. "over"
> doesn't stop with sp; you'll find further "over" methods in packages rgeos
> and spacetime.
>
> Hth,
> --
> 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
>