[GIS] Computing number of points in a raster grid cell in R

Question

I have a country shapefile and a dataframe with coordinates of bird sightings. My goal is to "rasterize" the shapefile and compute the number of points within each grid cell. Specifically, I want to end up with a grid-cell level dataframe containing lat/long of the centroid of each grid-cell, the district census code of that grid cell (which is contained in the shapefile attribute table), and the number of points falling in the grid cell. I would like to keep cells with no points in them and leave number of points as "NA". This is what I have so far. I get stuck after finding the district code of each cell. I guess the next step is to convert the raster to a polygon grid and then overlay the points on it? Is there an easier way?

# Load Packages
require(data.table)
require(rgdal)
require(rgeos)
require(sp)
require(raster)

# Load india 2011 districts
india.districts.2011 <- readOGR(paste(dist_shpf, "maps/india-district", sep=""), 
                            'SDE_DATA_IN_F7DSTRBND_2011', stringsAsFactors = F) 
# Read Bird Data
df.sample <- fread('ebd_IN_200001_201810_relOct-2018.txt')

# Initialize Grid
grid <- raster(extent(india.districts.2011))
res(grid) <- .25
proj4string(grid) <- proj4string(india.districts.2011) 
india.grid <- rasterToPolygons(grid) #SpatialPolygonsDataFrame

# Get centroids
india.grid.centroids = gCentroid(india.grid ,byid=TRUE)
india.grid.centroids$c_code_2011 <- 
over(india.grid.centroids,india.districts.2011)$district.code #SpatialPointsDataFrame

Answer 1

You can use raster::cellFromXY to get which raster cell a point is in, so you don't have to create polygons for each cell. Here's a simple reproducible example:

First I'll make a sample raster over this x and y range:

> xmn = 17;xmx = 23;ymn=42;ymx=49
> r = raster(matrix(1:12,3,4), xmx=xmx, xmn=xmn,ymx=ymx,ymn=ymn)

Next some sample points over the space:

> pts = data.frame(x=runif(10,xmn,xmx), y=runif(10,ymn,ymx))

See what they look like:

> plot(r)
> points(pts)

Next a function to do the work. Writing functions is a good way of making code that you can easily test on small data sets (like the one I just made) before running it on anything else. Always do this.

 pointcount = function(r, pts){
  # make a raster of zeroes like the input
  r2 = r
  r2[] = 0
  # get the cell index for each point and make a table:
  counts = table(cellFromXY(r,pts))
  # fill in the raster with the counts from the cell index:
  r2[as.numeric(names(counts))] = counts
  return(r2)
 }

Results:

 r2 = pointcount(r, pts)
 plot(r2)
 points(pts)