Overlaying re-projected points and polygons in R?

Question

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

1. 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.

2. 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)

Answer 1

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".

Answer 2

I couldn't find the polygons file at the site you shared. I'm using this one for illustration purposes:

http://mapas.ine.pt/download/files/2011/municipios/BGRI2011_1106.zip

library(sp)
library(rgdal)
# fread is fast & I love data.table syntax
library(data.table)

# readOGR is great because it automatically reads the prj file
Lisboa = readOGR('.', 'BGRI2011_1106')
pts = fread('Points.csv')
ptsSPDF = SpatialPointsDataFrame(
  pts[ , cbind(lon, lat)],
  # I always use this CRS for Google-derived lat/lon
  data = pts, proj4string = CRS('+init=epsg:4326'),
  match.ID = FALSE
)

#convert
ptsSPDF = spTransform(ptsSPDF, proj4string(Lisboa))

#plot for illustration
plot(Lisboa)
plot(ptsSPDF, add = TRUE, col = 'red')