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)
aldo_tapia
  • 13,560
  • 5
  • 30
  • 56
fabo
  • 61
  • 1
  • 4

1 Answers1

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
aldo_tapia
  • 13,560
  • 5
  • 30
  • 56
  • 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