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:

points<-read.csv(file = "/Users/input.csv", stringsAsFactors=FALSE)
coordinates(points)<- ~longitude+latitude  
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.


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


  coordinates(meuse) <- ~x+y

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


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)
  • Surely if the grid is regular you can work out the raster I,J for a point by arithmetic - subtract xmin, divide by cell width or similar...
    – Spacedman
    Commented Jan 29, 2019 at 21:02
  • 1
    @Spacedman, one would think but, I would not take it for granted that the cell_ids are uniform (they are in my example for simplicity sake). Also, one could have a grid with irregular "cell" sizes across the region. This would make it a non-systematic grid but, I have seen this type of data referred to as a grid. Commented Jan 29, 2019 at 21:08
  • 1
    The data function was just missing the closing. I do not see an error regarding plot margins, is this on your data on the example? Commented Jan 29, 2019 at 21:16
  • @JeffreyEvans I figures out the mistake I was doing on the code that I talked about in my previous comment. So I am able to replicate the first part of your code. However when I try to use the second part witch my data I get this message after running "pts_grid <- spatialEco::point.in.poly(pts, grid)" Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared) : st_crs(x) == st_crs(y) is not TRUE Commented Jan 29, 2019 at 21:54
  • 2
    @MarcelCampion this error indicates that your data is missing a projection string or are in different projections. Check both data using: proj4string(x) Commented Jan 29, 2019 at 21:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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