2

I am trying to intersect two spatial objects in R (point and polygon) while preserving the contained information in the dataframes of both.

However, the intersection produces an empty polygon, while both objects do perfectly overlay in QGIS 2.16.2.

After I checked the coordinates of the spatial polygons (test_polys), I find them strangely altered, such as coordinates like -44524.15 , 167543.0.

My efforts to change the coordinate system (Step 3) did not succeed. Some ideas how to fix it and explanations, where the coordinate modifications could originate from?

I use R.Studio Version 1.0.143 and rgdal 1.2-8 as well as raster 2.5-8.

Please find the input files in this dropbox folder

And the code I use:

# 1) Initialization -----------------------------------------------------

library(rgdal)
library(sp)
library(raster)

setwd("SET YOUR DIRECTORY") 
list.files()

# 2) Read data ----------------------------------------------------------

test_polys  <- readOGR(dsn = "C:/YOUR PATH/Desktop/Stackexchange", layer = "Tests_10subsections")
test_points <- readOGR(dsn = "C:/YOUR PATH/Desktop/Stackexchange", layer = "calc_centroids")

test_polys@polygons # Strange coordinates already at this step, such as "-44524.15 167543.0"

# 3) CRS Homogenization -------------------------------------------------

proj4string(test_polys)       <- CRS("+proj=longlat +datum=WGS84 +no_defs  +ellps=WGS84 +towgs84=0,0,0")
# ERROR "Error in `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) : Geographical CRS given to non-conformant data: -45033.0103 167554.5511"

test_polys@proj4string        <- CRS("+proj=longlat +datum=WGS84  +no_defs +ellps=WGS84 +towgs84=0,0,0")
test_points@proj4string       <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

    test_points  <- spTransform(test_points, CRS(test_polys))
# Error in CRS(test_polys) : 
# no method for coercing this S4 class to a vector
# In addition: Warning message:
# In is.na(projargs) : is.na() applied to non-(list or vector) of type   #'S4'

# 4) Intersection while preserving point data adding polygon data -------

intersected  <- raster::intersect(test_points, test_polys)
# EMPTY "data" with: [1] id         Roof_area  DTMN11     FR11       SEC11      SS11       BGRI11     LUG11      LUG11DESIG
#<0 rows> (or 0-length row.names)
2

Both layers don't have the same CRS, you have to reproject one layer prior to intersect. In QGIS you can see both layers overlapping with on-the-fly (OTF) projection of vector and raster layers, but also you can't intersect in QGIS if you don't have the same CRS.

For R, just simply use spTransform():

test_polys2 <- spTransform(x = test_polys, CRSobj = test_points@proj4string)

intersected  <- raster::intersect(test_points, test_polys2)

intersected
## class       : SpatialPointsDataFrame 
## features    : 101 
## extent      : -8.669497, -8.612373, 41.14529, 41.17578  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 9
## names       : id, Roof_area, DTMN11, FR11, SEC11, SS11,      BGRI11,  LUG11, LUG11DESIG 
## min values  :  1,        55,   1312,   01,   001,   03, 13120100103, 007437,      Porto 
## max values  : 99,      1277,   1312,   15,   014,   03, 13121500403, 007437,      Porto 
  • Great thank you, that was working. Any idea why readOGR produces the mutated coordinates as shown under point 3) in my explanation? – fabo Aug 25 '17 at 12:17
  • But they are not mutated. It is a different coordinate system (Transverse Mercator). Check the relevant .prj file or test_polys@proj4string (right after reading the data into R). – Janina Aug 25 '17 at 12:34
  • As @Janina says, they aren't mutated. The CRS are different, one considers degrees as unit and the other one, meters. Take a look into this question – aldo_tapia Aug 25 '17 at 12:46

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.