3

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:

enter image description here

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?

2

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

  • Thanks @Guzmán, it was a very clear and great answer, I was able to solve it! – Davido Feb 16 '17 at 14:46
  • @OscarMolina, you are welcome! Please, check the answer as an acceptable answer! =) – Guzmán Feb 16 '17 at 15:08

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.