Joinn shapefile attribute ID to point in R

Question

I have a CSV file (about 1,000,000 georeferenced points) which contains point_id/latitude/longitude. I also have a shapefile which represents a grid cell_id/xmin/xmax/ymin/ymax. I am trying to join the grid attributes (specifically the cell_id) to the points from my CSV file. Basically I want to know in which cell the points fall.

So my output would look like that point_id/latitude/longitude/cell_id.

Here is the code I have been working on based on this post and that post:

library(sp)
points<-read.csv(file = "/Users/input.csv", stringsAsFactors=FALSE)
coordinates(points)<- ~longitude+latitude  
library(sf)
grid <- read_sf("/Users/grid.shp")
a.data <- raster::intersect(points,grid["cell_id",,,])
points$cell_id <- a.data$cell_id
write.csv(points, file = "/Users/points.csv")

However this does not gives me the cell_id but instead returns the point_id twice.

I have also tried using over but it yielded an error message.

Answer 1

You are mixing and matching different object classes (sp and sf). The raster::intersect function may be able to deal with this but, it is a dangerous assumption and is not good practice. I would use the raster::shapefile function to read grid.shp (resulting in an sp object) and not attempt an index to limit to cell_id. Just let the function work and then deal with the resulting @data data.frame slot. You are trying to be too cleaver and it is causing problems. This really is a fairly simple point in polygon analysis.

First, create some example data (meuse point data) and r (polygon grid).

library(sp)
library(raster)

data(meuse)
  coordinates(meuse) <- ~x+y

r <- raster(extent(meuse), resolution = 100)
  r[] <- 1:ncell(r)
r <- as(r, "SpatialPolygonsDataFrame")
  names(r@data) <- "cell_id"

plot(r)
  plot(meuse,col="red",pch=19,add=TRUE)

Now, we can use the point.in.poly function in the spatialEco package. This will intersect the point and polygon feature geometries and return the associated polygon attributes. By default the returned object will be a SpatialPointsDataFrame however, you can use sp = FALSE to result in an sf class object.

meuse_grid <- spatialEco::point.in.poly(meuse, r)
  names(meuse_grid)
  meuse_grid$cell_id