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  

Solution 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    

Solution 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')