[GIS] Overlaying re-projected points and polygons in R

coordinate systemcorruptoverlayrshapefile

I am trying to overlay point data with polygons in R in order to find the polygons which intersect with the points provided.

The assciated prj file from the polygons shows their coordinate reference system: "DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]]".

The file is available here http://mapas.ine.pt/download/index2011.phtml.
Then I have sample points. The point coordinates come from the Google API and are transformed to WGS 1984.

The point data set is here https://www.dropbox.com/s/1r1rzmhjhqcsarv/Points.csv?dl=0.

But when I try to overlay both, an dataframe is returned filled with NAs.

When I plotted both files in QGIS for visual checking, I found that all most all points are aggregated in one polygon. However they should spread over Portugal,(even beyond the polygon outcrop which is Lisbon)

There should be a corruption of the coordinate system at one point, but I cannot find it!

I set a common coordinate reference system for polygons and points (WGS 1984) and repaired the polgyon with cleangeo. Below the commented code I used.

(A) Load points and polygons

(A1) Define common CRS

CRS.new <- CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

(A2) Read points from csv and build SpatialPointDataFrame

Points <- read.csv("Points.csv", header = TRUE, sep = ";")                         

Points_coords <- Points[,c(3,4)]
Population <- SpatialPointsDataFrame(coords = Points_coords, data = Points,proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

(A3) Read polygons from BGRI dataset

Portugal_subsections <- readShapeSpatial("BGRI11_LISBOA.shp")  
Portugal_subsections@proj4string <- CRS.new

(A3) Fix polygons with cleangeo

report <- clgeo_CollectionReport(Portugal_subsections) ## Check topology error
clgeo_SummaryReport(report)

Portugal_subsections_clean <- clgeo_Clean(Portugal_subsections) ## Clean topology errors

report2 <- clgeo_CollectionReport(Portugal_subsections_clean)  ## Check again
clgeo_SummaryReport(report2)

(B) Saving files

(B1) Portugals subsections

Portugal_subsections_clean@proj4string  <- CRS.new
writeOGR(obj=Portugal_subsections_clean, dsn = ".", Portugal_subsections_clean", driver="ESRI Shapefile")

(B2) Saving Population as point shape

colnames(Population@data) <- c("Index","Zip", "lon", "lat", "Population")
writeOGR(obj=Population, dsn = ".", "Population", driver="ESRI Shapefile")

(C) Overlay Population with BGRI subsections

Portugal_subsections_clean@proj4string <- CRS.new  ## Assure same CRS (WGS 84)
Population@proj4string <- CRS.new  ## Assure same CRS (WGS 84)

a <- over(Population, Portugal_subsections_clean)

Best Answer

Try the commented and reproducible example below! I use the sf package because is much faster than rgdal to open Shapefiles. Also, I used the great package mapview for interactive visualization in R. I transformed the points object to polygon's default projection before doing the intersection.

Note: I downloaded the Points.csv file from Dropbox link first and saved it to working directory.

# Load libraries ----------------------------------------------------------

library('sf')
library('mapview')

# Download data -----------------------------------------------------------

# Download polygon
download.file("http://mapas.ine.pt/download/files/2011/nuts2/lisboa2011.zip", destfile = "lisboa2011.zip")
unzip("lisboa2011.zip")

# Load data ---------------------------------------------------------------

polygons <- sf::read_sf("BGRI11_LISBOA.shp")
points <- read.csv(file = "Points.csv", sep = ";")

# Convert points df to sf object
points <- st_as_sf(points, coords = c("lon", "lat"))
st_crs(points) <- "+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0" # assign CRS to points

# Plot data
mapview::mapview(polygons, zcol = "DTMN11") + mapview::mapview(points, color = "white", col.regions = "black")

# Overlay Population with BGRI subsections --------------------------------

# Transform points data
pointsT <- st_transform(points, crs = st_crs(polygons))
print(st_crs(pointsT))

# Plot transformed data to check is ok
mapview::mapview(polygons, zcol = "DTMN11") + mapview::mapview(pointsT, color = "white", col.regions = "black")

# Intersection
overlayOutput <- st_intersects(pointsT, polygons)

# Build intersected polygons object
intersectedPolys <- polygons[unlist(overlayOutput), ]

mapview::mapview(polygons, zcol = "DTMN11") + mapview::mapview(pointsT, color = "white", col.regions = "black") + mapview::mapview(intersectedPolys, color = "black", col.regions = "black")

Note: the intersected polygons are colored in "black".

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)