Exportación de superposición espacial de paquete de R PBSmapping

Question

Después de darse cuenta de que el SEGURO de la biblioteca en el Césped y el paquete de R rgeos es increíblemente lento, me encontré con el packe "PBSmapping" hacer la misma cosa, una Unión de dos Polígonos y cortar las piezas con la no superposición de áreas.

El problema es exportar el creado Polyset dentro de R. Yo uso el siguiente código:

#require 
require("PBSmapping")
require("maptools")
require("rgdal")

# import .shp
area <- readShapeSpatial("area")
area.ps <- combinePolys(SpatialPolygons2PolySet(area))

area_WAL <- readShapeSpatial("area_WAL")
area_WAL.ps <- combinePolys(SpatialPolygons2PolySet(area_WAL))

# UNION
union.area <- combinePolys(joinPolys(area_WAL.ps, area.ps,"UNION", maxVert = 1.0e+06))

#export polyset as shape
export <- PolySet2SpatialPolygons(union.area)

... y aparece este mensaje de error:

In PolySet2SpatialPolygons(union.area) :
unknown coordinate reference system

Como CRS he usado UTM y GK, ni de tanto trabajó.

Hay otra forma de exportar un Polyset de Clase?

Cualquier comentario o ayuda sería muy apreciada!

Answer 1

¿Estás seguro de que rgeos y GEOS es muy lento? Dudo de eso, es altamente optimizado código C. Además, el paquete PBSmapping utiliza la GPC (General Polígono Clipper) de la biblioteca, que tiene algunas restricciones de la licencia. También afirman que desea "una Unión de dos Polígonos y cortar las piezas con la no superposición de áreas". Esto suena como una intersección, que se especifica como "INT"

Antes de atrevidamente afirmar que los GEOS es muy lento, vamos a hacer un experimento, utilizando un SEGURO de intersección, una GPC intersección y una GPC de la unión en caso de que es lo que realmente quería:

# Simple squares
p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
p2 = readWKT("POLYGON((0.5 0,1.5 0,1.5 1.5,0.5 1.5,0.5 0))")

#GEOS Intersection
pI <- gIntersection( p1 , p2 )

#Which look like
plot(p1 , col = "#ffd9d9" , xlim = c(0,2) , ylim = c(0 , 2 ) )
plot(p2 , col = "#d9d9d9" , add= T )
plot(pU , add = T , lty = 2 , bor = 2 , lwd = 2 )

# PBS equivalent
p1.ps <- SpatialPolygons2PolySet(p1)
p2.ps <- SpatialPolygons2PolySet(p2)

#GPC Intersection
pI.ps <- joinPolys(p1.ps, p2.ps , "INT" )

#GPC Union
pU.ps <- joinPolys(p1.ps, p2.ps , "UNION" )

# And if we benchmark this
library(microbenchmark)
microbenchmark( gIntersection( p1 , p2 ) , joinPolys(p1.ps, p2.ps , "UNION" ) , joinPolys(p1.ps, p2.ps , "INT" ) , times = 100 )

#Unit: microseconds
#                            expr      min       lq   median        uq       max neval
#           gIntersection(p1, p2)  311.117  359.586  424.504  447.6015   580.586   100
#                   gUnion(p1, p2) 394.259 399.9645 404.6595 411.6865 660.992   100
# joinPolys(p1.ps, p2.ps, "UNION") 3661.995 3862.172 4241.355 5260.2670 64874.749   100
#   joinPolys(p1.ps, p2.ps, "INT") 3625.948 3790.256 4042.752 4931.4790 64758.213   100

A través de 100 pistas, la GPC del motor y PBSmapping paquete es de alrededor de 10 veces más lento! En definitiva, evitar la posible pegajosa problemas de licencias dependiendo de su caso de uso, y ganar algo de velocidad mediante el uso de rgeos!

Editar

OP quiere mantener a las partes que no se superponen. Lo que podría ser una diferencia

pD1 <- gDifference( p1 , p2 )
pD2 <- gDifference( p2 , p1 )
pD <- gUnaryUnion( pD1 , pD2 )
plot( pD , lty = 2 , col = "#d93a3a" )

O si desea mantener tanto los bits que se superponen y los bits que no se superponen (una simple unión), entonces

pU <- gUnion( p1 , p2 )
plot(pU , col = "#d9ffd9")