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.
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)
Best Answer
Well, in case anyone has this problem, here's a solution I just realized that works. By simply:
before doing the over() and it works perfectly. I'm letting the question stay here, in the hope that it might help someone in the future. :)