Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

Any comment or help, would be very apreciated!

share|improve this question
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 Apr 2 '13 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 Apr 2 '13 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 Apr 2 '13 at 15:11
    
Have a look at: grass.osgeo.org/grass64/manuals/v.select.html But many thanks for the alternatives! –  Marco Apr 3 '13 at 7:13

1 Answer 1

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

share|improve this answer
    
@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? –  Simon O'Hanlon Apr 14 '13 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 Apr 14 '13 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 –  Simon O'Hanlon Apr 14 '13 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. –  Simon O'Hanlon Apr 14 '13 at 20:37

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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