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

As I replied before (and I know with complete certainty, because I wrote 
and maintain the code), this is a warning, not an error. You (no 
affiliation, I have no idea what work/research you are doing) and 
everybody else using R spatial software needs to be aware that Proj4 
strings (the ones starting with "+proj=") are almost defunct and should 
not be relied on. Hence the warning. If your workflow is affected by 
degradation as described in 
https://www.r-spatial.org/r/2020/03/17/wkt.html and 
https://rgdal.r-forge.r-project.org/articles/CRS_projections_transformations.html, 
you will suffer increased imprecision by about 150m.

You should in principle trust information from package authors and 
maintainers much more than random google search hits. I do not contribute 
to stackoverflow, and while the answer you refer to does concern this 
issue, it is basically wrong in saying that regular users do not need to 
pay attention. I have noted why in that thread. Downgrading R/packages 
makes absolutely no sense, and you should only consider installing old 
versions of R and packages to reproduce old results.

Again, you  are not seeing any errors in your example, only two warnings 
in the released version of raster, which are (wrongly) suppressed in the 
development version.

Please demonstrate that anything here is an error, and show the output of 
traceback in your case, which the reprex obviously does not capture. Your 
reprex also contains totally irrelevant calls to raster::shapefile(), 
which just writes the polygons to disk but are otherwise without effect.

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
>
> ´´´
>
> We have found in the following link:
>
> https://gis.stackexchange.com/questions/364667/r-4-0-1-not-sure-i-understand-this-message-warning-message-in-proj4stringx
>
> some comment about this problem: “This is not related to R 4.0.1 but
> to rgdal 1.5-8 and the migration to gdal 3 and proj 6. This is a very
> long and complex process that impact hundreds maybe thouthands of
> packages. All the packages are not yet up-to-date with what is coming.
> As a regular user you do not need to pay too much attention on these
> warnings but should be aware of what is coming in the geospatial
> world.
>
> You can have a look to ?rgdal::set_thin_PROJ6_warnings() to eliminate
> those warnings.”
>
> However, this warning for us is a problem since it does not allow us
> to perform the intersection operation between the 2 layers.
>
> What we have tried is to pass our R and rgdal to a previous version
> (Version rgdal 1.4-8) and R version 3.6.1 (2019-07-05) - "Action of
> the Toes" but we would like to know if this is the way to go or there
> other options.
>
> Currently we have managed to get the 3.6.1 R version and but we have
> trouble passing to the Rgdal version 1.4-8, hoping that this version
> does not have these problems with the warning that we can not perform
> the operation of intersection.
>
> Regarding how to install an older rgdal version we found an interesting link
>
> (https://stackoverflow.com/questions/51749860/installing-rgdal-backdate-rgdal-version-versus-update-gdal-version-and-how)
> that explains how to do it with the following commands:
>
> require(devtools)
> install_version("rgdal", version="1.4-8")
> But we get the following error:
> ERROR: lazy loading failed for package 'rgdal'
> * removing 'C:/Users/pepe/Documents/R/win-library/3.6/rgdal'
> * restoring previous 'C:/Users/pepe/Documents/R/win-library/3.6/rgdal'
> Error: Failed to install 'unknown package' from URL:
>  (converted from warning) installation of
> package ‘C:/Users/pepe /AppData/Local/Temp/RtmpqCxtfw/remotes26c7f616614/rgdal’
> had non-zero exit status
>
> Thank you for your help,
>
> _______________________________________________
> 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