[R-sig-Geo] gUnion: no outgoing dirEdge found

[Roger Bivand]

On Fri, 17 Jun 2011, Brian J. Stults wrote:

> Hello,
>
> I am trying to create a union and intersection between two shapefiles.
> They consist of census tracts in the state of Delaware - one file for
> 1990 and the other for 2000.  They are available here:
>
> http://www2.criminology.fsu.edu/~stults/misc/delaware_tracts.zip
>
> I would like to generate a shapefile that includes all identical
> polygons between the two layers as well as all portions that intersect
> (i.e. "union" in Qgis).  I believe gUnion is the function I need to use.

I can replicate the GEOS error:

> union <- gUnion(de1990, de2000)
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_union") :
   TopologyException: no outgoing dirEdge found at -77.0822 38.9222

but in:

> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
...

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgeos_0.1-8     stringr_0.5     maptools_0.8-9  lattice_0.19-26
[5] sp_0.9-83       foreign_0.8-44

loaded via a namespace (and not attached):
[1] grid_2.13.0  plyr_1.5.2   tools_2.13.0

no error exit - is this a script running on a node, if so, could you catch 
errors with try() or similar means? I am in doubt about whether your 
approach is sensible, and suspect that running gEquals(de1990, de2000, 
byid=TRUE) first to establish the duplicates, then tackle the ones which 
are not seen as equal. Most of the polygons are the same:

eq <- gEquals(de1990, de2000, byid=TRUE)
sum(apply(eq, 2, any))
sum(apply(eq, 1, any))
plot(de1990, col="cyan", border="transparent")
plot(de1990[apply(eq, 2, any),], add=TRUE)
nde1990 <- de1990[!apply(eq, 2, any),]
nde2000 <- de2000[!apply(eq, 1, any),]

then step through int row by row (or column by column) to union the ones 
for which the correct predicate is TRUE for that combination. I'm unsure 
whether you want to dissolve to the 2000 file or the other way round, or 
to create sub-units permiting either to be reconstructed, so I can't 
suggest further steps now. The predicates should let you find the 
geometries on which topology operations are to be performed, which should 
then work provided that the shapefiles are clean. The trick seems to be in 
finding the right predicates, then the right operation, but it needs 
interactive analysis to work through use cases and instrument a script to 
be sure that the correct joins are being made.

Hope this helps,

Roger


>
> Here is my code:
>
> library(maptools)
> library(rgeos)
>
> de1990 <- readShapePoly("shapefiles/subset_1990_de.shp")
> de2000 <- readShapePoly("shapefiles/subset_2000_de.shp")
>
> union <- gUnion(de1990, de2000)
>
>
> I receive the following error (full output attached):
>
>> union <- gUnion(de1990, de2000)
> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_union") :
>  TopologyException: no outgoing dirEdge found at -77.0822 38.9222
> Calls: gUnion -> RGEOSBinTopoFunc -> .Call
> Execution halted
>
> I confirmed that the polygons are valid using gIsValid, but maybe that's
> not what this error is indicating.
>
> I am able to successfully do the union in Qgis.  However, I eventually
> need to do this for the entire U.S. for several decades, and just one
> union of the full U.S. takes several days to complete on my desktop.  I
> would rather run it in R on my university's supercomputer, which would
> significantly reduce the run-time.
>
> I am hoping that this does not require editing the shapefile, since I
> could not possibly do that for all tracts in the U.S.
>
> Thanks for any suggestions,
> Brian
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no