0

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.