Skip to content

Unable to perform intersection operations with vectors

4 messages · Coyo Tepe, Roger Bivand

Original
Coyo Tepe
# 
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

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,
Roger Bivand
#
On Mon, 13 Jul 2020, Coyo Tepe wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      Coyo Tepe
    

    

      
      #
    

  


  

      
Dear Doctor Roger Bivand,

First I have to say thank you for your quick response and secondly I
am sorry for the confusion. I did not mean you to think I don?t trust
your comments. Not at all, you are the reference in this matter..


I first sent an incomplete message and when I realized it I sent the
complete one,  that?s why the message was sent twice and I did not
read your reply before sending this second message.


Regarding your message, you are right, this is a warning, not an error
and I agree with you that is much better to have warnings that
problems at the end of the process and not knowing what could be
wrong.


I was not aware thae using the Proj4 strings should not be relied on
and of course you are right about the comment with shapefile() call,
this should not be part of the REPEX. Thanks also for the
recommendation regarding the downgrading or R packages. I have read
the links you shared regarding the GDAL and PROJ development and is
more clear to me now.


Thank you!!  and again  I am really sorry for the misunderstanding,

El mar., 14 jul. 2020 a las 6:36, Roger Bivand (<Roger.Bivand at nhh.no>) escribi?:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 15 Jul 2020, Coyo Tepe wrote:

        

            

      
      
      
    

        
OK, not a problem.

            

      
      
      
    

        
I mis-understood this - you might have added a comment at the head of your 
post. Had I looked at 
https://stat.ethz.ch/pipermail/r-sig-geo/2020-July/thread.html, I would 
have seen that you were replying to yourself, not to my reply.

            

      
      
      
    

        
Excellent! Now you and others following the thread have been alerted. What 
we do not know is which workflows will be affected - if we did, we could 
check and only warn those affected. So awareness of changes in 
sp::CRS/raster::crs/sf::st_crs will be needed to let users check whether 
their current workflows (written before the changes) are affected.

Roger