5

After realising that the GEOS-library in Grass and the R package rgeos is incredibly slow, i found the packe "PBSmapping" to do the same thing, a Union of two Polygons and to cut the parts with no overlay-areas.

The problem is to export the created Polyset within R. I use the following code:

#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)

... and get this error message:

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

As CRS I've used UTM and GK, neither of both worked.

Is there another way to export a Polyset-Class?

4
  • 1
    Has your original file a correct CRS ? I don't see that you specifiy it in your code. What is the output of 'print(proj4string(area))' ?
    – Curlew
    Commented Apr 2, 2013 at 10:13
  • "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel" ... this is the CRS with which i projected the shapefiles and it's also the output of proj4string.
    – Marco
    Commented Apr 2, 2013 at 12:27
  • 1
    GRASS GIS don't use GEOS but a completely different topological algorithm. The union function of GEOS was supplanted by the Cascaded Union function, much faster (ported from the JTS, Java Topological Suite) which is used in PostGIS, or in Python modules as Shapely. Try OpenJump which uses the JTS with cascaded Union
    – gene
    Commented Apr 2, 2013 at 15:11
  • Have a look at: grass.osgeo.org/grass64/manuals/v.select.html But many thanks for the alternatives!
    – Marco
    Commented Apr 3, 2013 at 7:13

1 Answer 1

3

Are you sure that rgeos and GEOS is very slow? I doubt that, it is highly optimised C code. In addition the package PBSmapping uses the GPC (General Polygon Clipper) library, which has some restrictive licensing. You also claim that you want "a Union of two Polygons and to cut the parts with no overlay-areas". This sounds like an intersection, which is specified as "INT"

Before you boldly claim that GEOS is very slow, let's try an experiment, using an GEOS intersection, a GPC intersection and a GPC union in case that is what you really wanted:

# 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 )

enter image description here

# 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

Across 100 runs, the GPC engine and PBSmapping package is around 10 times slower! In short, avoid possible sticky licensing issues depending on your use case, and gain some speed by using rgeos!

Edit

OP wants to keep the parts that do not overlap. Which could be a difference

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

enter image description here

Or if you want to keep both the bits that overlap and the bits that don't overlap (a simply union), then

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

enter image description here

4
  • @Marco did this help and do you get similar results on your system? Were you looking for an intersection, or perhaps I misunderstand your question? Commented Apr 14, 2013 at 6:27
  • I'm rather looking for an union, but with non-overlapping/touching parts left behind. This is your suggested intersection: Intersection ... and this is what i'm looking for: Union The red parts should be maintaned. I have benchmarks myself from both of the libraries, GEOS and GPC, and in my case the Union with GEOS was significant slower (GPC: 5 minutes, GEOS: around 2 days).
    – Marco
    Commented Apr 14, 2013 at 9:45
  • @Marco so what you want is the difference? I don't understand when you say you want NON-overlapping parts left behind. Does that mean you want to remove overlapping parts? Otherwise you just want to join the full extent of both polygons. I have added a GEOS union. It is far quicker than using GPC from the PBS package. For further help or insight I suggest making your shapefiles available online if you are able? Cheers Commented Apr 14, 2013 at 20:30
  • gUnion, still appears to be much faster in my simple test. Please post the code you are using for unioning polygons with rgeos if you want anyone to look at optimising it. It is very odd that rgeos would take 2 days over gpc. Commented Apr 14, 2013 at 20:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.