[R-sig-Geo] Unable to perform intersection operations with vectors

[Roger Bivand]

On Mon, 13 Jul 2020, Coyo Tepe wrote:

> Platform and versions
> Windows 64 bit
> R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Irrelevant > Rstudio Version 1.3.959
> rgdal: 1.5-12 / gdal: 3.0.4
> rgeos version: 0.5-3, (SVN revision 634)
> GEOS runtime version: 3.8.0-CAPI-1.13.1
> Linking to sp version: 1.4-2
> Polygon checking: TRUE

> Problem
>
> We are trying to perform an intersection operation between two vector
> layers in shapefile format and we get an error:
>
> Warning message:
>
> In proj4string(x) : CRS object has comment, which is lost in output
>

No, this is a warning not an error (unless you set options(warn=2)). Read 
for example https://www.r-spatial.org/r/2020/03/17/wkt.html and 
https://rgdal.r-forge.r-project.org/articles/CRS_projections_transformations.html 
to begin to grasp underlying changes that affect everyone's workflows. You 
can no longer assume that using Proj4 strings gives a satisfactory 
representation of CRS. With recent versions of raster, you see:

cat(comment(crs(capa_a)), "\n")

giving the WKT2 (2019) representation.

The warnings are given because the released raster S4 "intersect" method 
for signature (x='SpatialPolygons', y='SpatialPolygons') uses 
sp::proj4string() rather than the actual slot in a test. The code is (from 
getMethod("intersect", c('SpatialPolygons', 'SpatialPolygons'))):

     if (!identical(proj4string(x), proj4string(y))) {
         warning("non identical CRS")
         y using proj4string <- x using proj4string
     }

where each use of sp::proj4string() generates a warning. Here, they are 
identical, and the CRS slot of the output is handled correctly by rgeos 
internally. In the development version of raster, there is no use of 
sp::proj4string(), so no warning is made; however rgeos sees objects with 
no assigned projection, stopping it complaining when a user uses strictly 
planar methods with geographical coordinates. This seems convenient but 
is wrong.

The warnings are in place because workflows may be affected by Proj4 
strings no longer working as before, and being in practice abandoned. In 
sp-based workflows, they are still present, but are not used for 
transformation. Warnings are much better than unannounced changes in 
output.

Roger


>
>
> Below you can have a REPEX:
>
> ´´´
>
> library(rgdal)
>
> library(raster)
>
> library(rgeos)
>
> #layer a
>
> bbox_a<-c(497555,500000,137000,138000)
>
> extent_a<-extent(bbox_a)
>
> extent_a
>
>
>
> capa_a<-as(extent_a,'SpatialPolygons') #creamos un objeto tipo SpatialPolygons
>
> crs(capa_a)<- "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"
>
> crs(capa_a)
>
> shapefile(capa_a,paste("capa_a_re",sep="/"),overwrite=TRUE)
>
> plot(capa_a,col="green") #plot(capa_a,col="green")
>
> #layer b
>
> bbox_b<-c(497995,500050,137100,138100)
>
> extent_b<-extent(bbox_b)
>
> extent_b
>
> capa_b<-as(extent_b,'SpatialPolygons') #creamos un objeto tipo SpatialPolygons
>
> crs(capa_b)<- "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"
>
> crs(capa_b)
>
> shapefile(capa_b,paste("capa_b_re",sep="/"),overwrite=TRUE)
>
> plot(capa_b,add=T)
>
> #buffer
>
> b_buf<-buffer(capa_b,width=50)
>
> plot(b_buf,add=T) #plot(b_buf,add=T,col="blue")
>
> #intersect
>
> b_intersect<-intersect(capa_a,b_buf)
>
>
>
> #Warning messages:
>
> # 1: In proj4string(x) : CRS object has comment, which is lost in output
>
> # 2: In proj4string(y) : CRS object has comment, which is lost in output
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en