unable to plot two spatial objects with identical CRS in R

Question

I have two spatial objects, one is a spatial polygon object and the other one is a .csv file that I turned into a spatial points object. The first one is an official shape file from the chilean government for one of its communes, the other one was created by geocoding with the HERE API, street adresses of the same commune.

First I loaded the spatial polygon object with readOGR from the :

quilpue <- readOGR( dsn= getwd() , layer="quilpue-rgdal", 
                encoding = "UTF-8") 

Then I loaded the .csv file into R, and converted it into a spatial points object with the coordinates()function from the sp package.

pointsCoords<- read.csv("../quilpueR/quilpueLayer.csv", header = TRUE)
coordinates(pointsCoords) <- ~Longitude+Latitude

Then I checked the projection of each object.

proj4string(quilpue) 
proj4string(pointsCoords)

"+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" and NArespectively.

The only projection that work for pointsCoords was CRS("+init=epsg:3857"). Therefore I assigned that projection to quilpue

proj4string(pointsCoords) <- CRS("+init=epsg:3857")
quilpue_prj <- spTransform(quilpue, CRSobj = CRS(proj4string(pointsCoords)))

Nonetheless, when I check the extension of both objects with extent() from raster() package, they do not overlap.

extent(quilpue_prj)
class       : Extent 
xmin        : -7957703 
xmax        : -7946463 
ymin        : -3907594 
ymax        : -3898059 

extent(pointsCoords)
class       : Extent 
xmin        : -71498550 
xmax        : -71334950 
ymin        : -33133030 
ymax        : -32769810 

Therefore, when I try to plot them together, they do not overlap. i only get the plot of the first object that I choose to plot.

plot(quilpue_prj)
plot(pointsCoords, add = TRUE) 

To check if there was a problem with the shapefile, or .csv file, I opened both on Maptitude another GIS software, and it manage to automatically overlay them. I would like to be able to do the same in R.

Answer 1

I think you've gone through some unnecessary steps in order to answer your question.

In order to plot them together they simply need to be in the same CRS:

library(rgdal)    

quilpue <- spTransform(quilpue, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) #or any CRS you wish to use

pointsCoords <- spTransform(pointsCoords, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"")) #or any CRS you wish to use

plot(quilpue)
plot(pointsCoords, add = T)

A helpful thing to also do is to check the extent of your features in order to make sure they align. Sometimes it happens for the features in question to be in the same CRS, but due to processing or many conversion, they end up with a screwy extent. Check this with:

extent(quilpue)
extent(pointsCoords)

Answer 2

I manage to solve the problem but I don´t actually understand why it worked. After loading the .csv file and using

coordinates(pointsCoords) <- ~Longitude+Latitude 

to create the spatial points object, I used the projection() function from rasterpackage to assign it a projection:

projection(pointsCoords) = "+init=epsg:4326"

Then, I transformed the projection of the spatial polygon object quilpue, first to "+init=epsg:3857" and then to "+init=epsg:4326":

quilpue <- spTransform(quilpue, 
             CRSobj = CRS("+init=epsg:3857"))

quilpue <- spTransform(quilpue, 
                       CRSobj = CRS("+init=epsg:4326"))

With bbox() I check the range of each spatial object:

bbox(pointsCoords)

            min       max
Longitude -71498550 -71334950
Latitude  -33133030 -32769810

bbox(quilpue)

    min       max
x -71.48526 -71.38429
y -33.09254 -33.02075

And notice they were very similar, and that pointsCoords was contained within quilpue. The only caveat was that quilpue coords had a "." after the first two digits, so I used gsub to add "." to the coords in pointsCoords.

dfcoords <- as.data.frame(pointsCoords@coords)
dfcoords$Longitude <- as.numeric(gsub("([[:digit:]]{6,6})$", ".\\1", 
                      dfcoords$Longitude)) 
dfcoords$Latitude <- as.numeric(gsub("([[:digit:]]{6,6})$", ".\\1", 
                      dfcoords$Latitude)) 
coordinates(dfcoords) <-  ~Longitude+Latitude

And assign the modified coordsto the original ones.

pointsCoords@coords <- dfcoords@coords

Then I was able to use over() and plot the spatial objects.

df_over <- over(quilpue_prj, pointsCoords)
plot(quilpue)
plot(pointsCoords, add = TRUE)