[R-sig-Geo] Counting points within polygons

[Agustin Lobo]

Many thanks Maria for the very interesting pointer!

I've been able to find this solution within R:

> delme5 <- over(fireLC2[,5], BorneoLC)
> delme5[1,]
  LC    ID LC_2                      LCname LCabrv     AREA PERIMETER
1  7 17582    7 Plantation_Secondary Forest  PltSF 0.379161    35.415
just to remind you:
> fireLC2 at data[1,]
  LC ID                      LCname LCabrv nbp
1  7  1 Plantation_Secondary Forest  PltSF   1
> BorneoLC at data[1,]
  LC   ID LC_2         LCname LCabrv    AREA PERIMETER
0  4 4505    4 Lowland Forest    LwF 7.7e-05  0.036385

> nbp <- table(delme5$ID)
> nbp <- data.frame(ID=as.numeric(names(nbp)), nbp=as.numeric(nbp))
> nbp[1:3,]
    ID nbp
1 4758   1
2 4760   1
3 4836   2
> BorneoLC3 <- mijoin(BorneoLC,nbp,by.x=2, by.y=1)
> BorneoLC3 at data$nbp[is.na(BorneoLC3 at data$nbp)] <- 0
> BorneoLC3 at data[1:3,]
  LC   ID LC_2         LCname LCabrv     AREA PERIMETER ID nbp
0  4 4505    4 Lowland Forest    LwF 0.000077  0.036385 NA   0
1  8 4524    8 Lowland Mosaic    LwM 0.000346  0.094675 NA   0
2  4 4525    4 Lowland Forest    LwF 0.000020  0.018000 NA   0

> writeOGR(BorneoLC3,dsn="/media/Iomega_HDD/JACOB/ExFireBorneo/RFireBorneo/BorneoLC3",layer="BorneoLC3",driver="ESRI Shapefile")

The sahpe file is correctly displayed in qgis


Agus


2011/10/18 Maria Zwart <m.c.zwart at newcastle.ac.uk>:
> Hi,
>
> There is an easy way to do this, but it is outside R. I use GME (geospatial modelling environment). It can also run R commands however.
>
> The function in that program you are looking for is: countpntsinpolys
> More information on: http://www.spatialecology.com/gme/gmecommands.htm
>
> Best,
>
> Mieke
>
>>-----Original Message-----
>>From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-
>>project.org] On Behalf Of Agustin Lobo
>>Sent: 18 October 2011 15:25
>>To: r-sig-geo
>>Subject: [R-sig-Geo] Counting points within polygons
>>
>>Hi!
>>
>>Given a Sp Points DF and an Sp Polygons DF, I want to add a field of
>>number of points for each polygon
>>to the Sp Pol. DF table. This is what I'm doing (using objects from
>>the help page for sp::overlay()
>>(meuse is the Sp Points DF and srdf the Sp Pol. DF)
>>
>>> delme3 <- overlay(meuse,srdf)
>>> delme3
>>  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
>>1 1 1 1
>> [38] 1 1 1 1 1 1 2 1 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3
>>3 3 3 3
>> [75] 3 3 3 3 3 3 3 3 2 1 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2
>>2 2 2 2
>>[112] 2 2 2 2 2 2 3 2 2 2 2 2 1 2 2 2 2 1 1 2 2 1 2 2 2 2 3 3 3 3 3 3 3
>>3 3 3 3
>>[149] 3 3 3 3 3 3 3
>>newdata <- data.frame(srdf at data,nppoints=table(delme3))
>>srdf2 <- srdf
>>srdf2 at data <- newdata
>>
>>An alternative:
>>newdatap <- data.frame(meuse at data, nbp=rep(1,nrow(meuse)))
>>meuse2 <- meuse
>>meuse2 at data <- newdatap
>>delme4 <- overlay(meuse2[,13],srdf,fn=sum)
>>newdata <- data.frame(srdf at data,nppoints=delme4)
>>srdf2 <- srdf
>>srdf2 at data <- newdata
>>
>>I wonder if there is a better way or, best, an standard R spatial
>>function doing that already.
>>
>>Thanks!
>>
>>Agus
>>
>>_______________________________________________
>>R-sig-Geo mailing list
>>R-sig-Geo at r-project.org
>>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>