R: How can I count how many points are in each cell of my grid?

Question

I have made a reference grid, cells 50x50m, based on GPS locations of a collared animal. I want to do the equivalent to a spatial join in ArcGIS, and count the number of points in each cell.

I have made a reference grid, using a SpatialPointsDataFrame object (the data frame is already projected, using UTM coordinate system)

RESO <- 50 # grid resolution (m)
BUFF <- 500 # grid extent (m) (buffer around location extremes) 
XMIN <- RESO*(round(((min(dat.spdf$Longitude)-BUFF)/RESO),0))
YMIN <- RESO*(round(((min(dat.spdf$Latitude)-BUFF)/RESO),0))
XMAX <- XMIN+RESO*(round(((max(dat.spdf$Longitude)+BUFF-XMIN)/RESO),0))
YMAX <- YMIN+RESO*(round(((max(dat.spdf$Latitude)+BUFF-YMIN)/RESO),0))
NRW <- ((YMAX-YMIN)/RESO)
NCL <- ((XMAX-XMIN)/RESO)
refgrid<-raster(nrows=NRW, ncols=NCL, xmn=XMIN, xmx=XMAX, ymn=YMIN, ymx=YMAX) 
refgrid<-as(refgrid,"SpatialPixels")

To make sure my grid was in the same projection as the SpatialPoints:

proj4string(refgrid)=proj4string(dat.sp.utm) #makes the grid the same CRS as point

count.point function in adehabitatMA seems like the function that will do the trick

cp<- count.points(dat.spdf, refgrid)

But I get this error:

Error in w[[1]] : no [[ method for object without attributes

Is this not the right route to take to achieve my goal? Or how can I resolve this error? Or would the over function (sp package) be more suitable?

output from SpatialPointsDataFrame (dat.spdf)

dput(head(dat.spdf, 20))
structure(list(Latitude = c(5.4118432, 5.4118815, 5.4115713, 
5.4111541, 5.4087853, 5.4083702, 5.4082527, 5.4078161, 5.4075528, 
5.407321, 5.4070598, 5.4064237, 5.4070621, 5.4070555, 5.4065127, 
5.4065134, 5.4064872, 5.4056724, 5.4038751, 5.4024223), Longitude = c(118.0225467, 
118.0222841, 118.0211875, 118.0208637, 118.0205413, 118.0206064, 
118.0204101, 118.0209272, 118.0213827, 118.0214189, 118.0217748, 
118.0223343, 118.0227079, 118.0226916, 118.0220733, 118.02218, 
118.0221843, 118.0223316, 118.0198153, 118.0196021), DayNo = c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L)), .Names = c("Latitude", "Longitude", "DayNo"), row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 15L, 16L, 
17L, 18L, 19L, 20L, 21L), class = "data.frame")

And summary:

summary(dat.spdf)
Object of class SpatialPointsDataFrame
Coordinates:
           min      max
Longitude 610361.0 613575.5
Latitude  596583.5 599385.2
Is projected: TRUE 
proj4string : [+proj=utm +zone=50 +ellps=WGS84]
Number of points: 5078
Data attributes:
Latitude       Longitude       DayNo      
Min.   :5.396   Min.   :118   Min.   :  1.0  
1st Qu.:5.404   1st Qu.:118   1st Qu.: 92.0  
Median :5.406   Median :118   Median :183.0  
Mean   :5.407   Mean   :118   Mean   :182.6  
3rd Qu.:5.408   3rd Qu.:118   3rd Qu.:273.0  
Max.   :5.422   Max.   :118   Max.   :364.0  

Stridulous

Did you verify that refgrid is of class (or inherited class) SpatialPixels ? – Mariandi

Yes, seems to be (sorry about the formatting...) summary(refgrid) Object of class SpatialPixels Coordinates: min max x 609850 614100 y 596100 599900 Is projected: TRUE proj4string : [+proj=utm +zone=50 +ellps=WGS84] Number of points: 6460 Grid attributes: cellcentre.offset cellsize cells.dim s1 609875 50 85 s2 596125 50 76 – Stridulous

If you leave refgrid as a raster, you can build on table(cellFromXY(refgrid, dat.spdf)). – Escent

Thanks @jbaums, I have the table that looks like the number of points within each cell. How would I get to display spatially? – Stridulous

I can't test this now but maybe r <- raster(refgrid); r[] <- 0; tab <- table(cellFromXY(refgrid, dat.spdf)); r[as.numeric(names(tab))] <- tab – Escent

I think this is almost it @jbaums, but how would I link up the cell names from the table that have the count value now attributed to it, back to the cells corresponding GPS coordinates? – Stridulous

Answer 1

Here's one way to do it, first tabulating the frequency of cell numbers represented by points, then assigning these frequencies to the cells' values, and finally extracting the cells' coordinates and values.

library(raster)
r <- raster(xmn=0, ymn=0, xmx=10, ymx=10, res=1)
r[] <- 0
xy <- spsample(as(extent(r), 'SpatialPolygons'), 100, 'random')
tab <- table(cellFromXY(r, xy))
r[as.numeric(names(tab))] <- tab

Now we have something like this:

plot(r)
points(xy, pch=20)

We can extract the cells' coordinates with coordinates() and their values with values(r) or simply r[]:

d <- data.frame(coordinates(r), count=r[])

head(d)
##     x   y count
## 1 0.5 9.5     0
## 2 1.5 9.5     1
## 3 2.5 9.5     1
## 4 3.5 9.5     3
## 5 4.5 9.5     2
## 6 5.5 9.5     3    

Escent

Answer 2

The rasterize function can do that for you:

library(raster)
r <- raster(xmn=0, ymn=0, xmx=10, ymx=10, res=1)
xy <- spsample(as(extent(r), 'SpatialPolygons'), 100, 'random')

x <- rasterize(xy, r, fun='count')

Proliferous