Unable to perform intersection operations with vectors

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Unable to perform intersection operations with vectors

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Unable to perform intersection operations with vectors

Roger Bivand
Administrator
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@proj4string <- x@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
> [hidden email]
> 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: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway