[GIS] Find which grid cells contain a point

overlayrvector-grid

I have a dataframe with x,y coordinates. I've sampled 10 x,y coordinates from my dataframe:

y <-c(24.637032,21.569813,21.687,19.677183,24.72438,21.433136,24.637032,22.21065,21.443762,25.273882)

x <-c(-156.02597,-112.00287,-112.00287,-105.27803,-159.47324,-157.858,-156.026,-109.68626,-112.12681,-112.15393)

I've converted these points to a SpatialPoints objects:

library(sp)
myPoints <- data.frame(x,y)
coordinates(myPoints) <- c('x', 'y')
proj4string(myPoints) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

I've used these points to make a SpatialGridobject:

bb <- bbox(myPoints)
myPointsGrid <- GridTopology(cellcentre.offset =  bb[, 1] + (cs/2), cellsize = c(1, 1), cells.dim = ceiling(diff(t(bb))/cs))
myPointsGrid <- SpatialGrid(myPointsGrid, proj4string = "+proj=longlat +ellps=WGS84 +datum=WGS84")

I need two things 1. find which grid cells contain a point 2. Find the number of points in each grid cell. I tried…

over(myPointsGrid, myPoints)

…but this has thrown an error. How can I complete this task?

Best Answer

Assuming that your code is correct (I have given the value 1 for the cs variable), the error it gives is:

Error in function (classes, fdef, mtable)  :  unable to find an inherited method for function ‘over’ for signature ‘"SpatialGrid", "SpatialPoints"’

According to the help file, this is what the function is waiting as input:

x: geometry (locations) of the queries
y: layer from which the geometries or attributes are queried

It means that you are giving the wrong order of inputs, first is the point location followed by the layer. This is what you want:

> over(myPoints,myPointsGrid)
  [1]  59 268 213 330   1 222  59 215 268  48

with cs=1

Related Solutions

[GIS] Reprojecting Google Maps (latlong) to Irish Grid TM65 using spTransform in R

The Google Map uses EPSG:3857 (aka Google Mercator) for displaying the map, but the coordinates you gathered from the map are given in lat/lon degrees. So these are not in mercator projection, which has (not real) metres as units, and you have to set EPSG:4326 WGS 84 as source fro your collected data.

For clarification:

A projection transforms the surface of the earth in a form that can be measured in linear units like metres or feet. A geodetic CRS like EPSG:4326 uses degrees to describe the position on the globe. Every projection is based on a geodetic CRS. So if you have degrees as units, rip of the projection information, and use +proj=longlat with the +ellps, +towgs84 and/or +datum information. In most cases +datum contains +ellps and +towgs84 information internally.

[GIS] Problem points into a shape polygon

You have to transform your polygon (shapefile) Coordinate Reference System (CRS) with the function spTransform from the sp package. You can simply do this:

polygon.t <- spTransform(polygon, CRSobj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

I also give you a complete reproducible example below:

Get example data:

# Load libraries
library('rgdal')
library('raster')

# Get SpatialPolygonsDataFrame object example
polygon <- getData('GADM', country = 'URY', level = 0)

# Plot data
plot(polygon)

# Make sample points 
sample.spdf <- spsample(polygon, n = 200, type = "random")
sample.spdf <- spsample(sample.spdf, n = 200, type = "random") # get points outside polygon
sample.spdf <- SpatialPointsDataFrame(sample.spdf, data = data.frame("ID" = 1:200, "V2" = sample(x = 1:1000, size = 200))) # add random data

Plot data:

# Plot points and polygon
plot(polygon)
plot(sample.spdf, add = TRUE, pch = 19, cex = 1)
box()

# Interactive map with Mapview
library(mapview)
mapview(polygon) + mapview(sample.spdf)

Spatial Query: points in polygon with different CRS:

# Apply transformation to sample.spdf
sample.spdf.t <- spTransform(sample.spdf, CRSobj = "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Ask layers CRS
projection(sample.spdf.t) # [1] "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(polygon) # [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

pt.in.poly.sp <- sp::over(sample.spdf.t, polygon, fn = NULL) 
# Error: identicalCRS(x, y) is not TRUE
pt.in.poly.rgeos <- rgeos::gIntersection(sample.spdf.t, polygon) 
# spgeom1 and spgeom2 have different proj4 strings

# Plot points transformed + polygon
plot(polygon)
plot(sampleSP.t, add = TRUE, pch = 19, cex = 0.5)

Spatial Query: points in polygon with same CRS:

# Ask layers CRS
projection(sample.spdf) # [1] [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(polygon) # [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

pt.in.poly.sp <- sp::over(sample.spdf, polygon, fn = NULL)
pt.in.poly.rgeos <- rgeos::gIntersection(sample.spdf, polygon)

# Plot points in polygon
plot(polygon)
plot(sample.spdf, pch = 19, cex = 1, add = TRUE)
plot(pt.in.poly.rgeos, pch = 19, cex = 0.5, col = 'green', add = TRUE)
box()

# Interactive map with Mapview
mapview(polygon) + mapview(sample.spdf) + mapview(pt.in.poly.rgeos)