[GIS] Problem points into a shape polygon

Question

I want to determine which points from a set of hundreds of geographic coordinates are located into a shape polygon. I was able to convert all the coordinates (latitude and longitude ) into a SpatialPointsDataFrame and I can read and plot the shapefile in R as well. A graph of my shp file together with the points is:

From similar questions I was trying to use the function point.in.polygon in the sp package as well as the the over() function.

Polygon: polygon1,   
Points: spdf

pt.in.poly <- sp::over( spdf, polygon1, fn = NULL) 

or:

new_shape <- point.in.poly(spdf, polygon1)

but in both cases the result I get is:

Error: identicalCRS(x, y) is not TRUE

The SpatialPointsDataFrame was created with the coordinates: CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") and the shapefile have the CRS arguments: +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 , so they seem to be the same, I am not sure, Im not an expert on this.

If its a problem of coordinates system how can I fix it? or what am I doing wrong?

Answer 1

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)