# Computing number of points in a raster grid cell in R

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)

'SDE_DATA_IN_F7DSTRBND_2011', stringsAsFactors = F)

# 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
``````

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)
`````` 