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

[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 using nhh.no>) escribió:
>
> 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