'CGAL::Precondition_exception' related to self-intersection

[ldemaz]

Hello,

We have been using pprepair quite successfully recently to clean up many hand digitized polygons. However, we got this error on one recently.

$ wget(https://dl.dropboxusercontent.com/u/31890709/mt_21FDWIOV1S7STSB8SJS3.zip) $unzip(mt_21FDWIOV1S7STSB8SJS3.zip)

$ pprepair -i "mt_21FDWIOV1S7STSB8SJS3.shp" -o "mt_21FDWIOV1S7STSB8SJS3_fix.shp" -fix

Adding a new set of polygons to the triangulation... Path: mt_21FDWIOV1S7STSB8SJS3.shp Type: ESRI Shapefile Layers: 1 Reading layer #1 (11 polygons)...

  int         gid
    string      name
    string      completion
    string      assignment
Feature #2 (6 vertices): self intersecting outer boundary #0. Split.
Created 1 rings.
Feature #10 (16 vertices): self intersecting outer boundary #0. Split.

terminate called after throwing an instance of 'CGAL::Precondition_exception' what(): CGAL ERROR: precondition violation! Expr: hierarchy.is_constrained_edge(va,vb) File: /usr/include/CGAL/Constrained_triangulation_plus_2.h Line: 233 Abort (core dumped)

Thanks in advance for your thoughts on this!

[hugoledoux]

I am not sure why CGAL throws an error there, we'll investigate.

But the problem is that the last polygon in the SHP is invalid, it's WKT is:

POLYGON((29.59067730619167236 -28.787288587651922,29.58994774533908156 -28.78723217025416403,29.58932547284739201 -28.78683724761136276,29.58878903104476521 -28.78627306981286793,29.58881048871658592 -28.78557724632377557,29.58891777707746584 -28.78414797270364645,29.58951859189644296 -28.78358378035597198,29.59033398343718702 -28.78364019972809018,29.59084896756798599 -28.78373423194659253,29.59091334058433986 -28.78373423194659253,29.59087042524070554 -28.78373423194659253,29.59087042524070554 -28.78379065123732516,29.59087042524070554 -28.78377184481073314,29.59082750989617949 -28.78386587691133514,29.59082750989617949 -28.78384707049752933,29.5907631368799855 -28.78558664937470724,29.59067730619167236 -28.787288587651922))

If you use prepair we have newer/better code to repair single polygons and there it works, the result is:

MULTIPOLYGON (((29.590677306191672 -28.787288587651922,29.590763136879985 -28.785586649374707,29.590827509896179 -28.783847070497529,29.590827509896179 -28.783865876911335,29.590870425240706 -28.783771844810733,29.590870425240706 -28.783734231946593,29.590848967567986 -28.783734231946593,29.590333983437187 -28.78364019972809,29.589518591896443 -28.783583780355972,29.588917777077466 -28.784147972703646,29.588810488716586 -28.785577246323776,29.588789031044765 -28.786273069812868,29.589325472847392 -28.786837247611363,29.589947745339082 -28.787232170254164,29.590677306191672 -28.787288587651922)))

So that's your solution for the time being, cheers for reporting that bug.

[ldemaz]

Hi Hugo,

Many thanks for your response, and thanks for the prepair solution. If I could follow up to get your advice, we are using pprepair to clean anything from individual polygons to several polygons. We call pprepair out of R using two functions:

library(rgdal) library(rgeos)

callPprepair <- function(dirnm, spdfinname, spdfoutname, crs = crs) {

Function to make system call to polygon cleaning program pprepair

Args:

dirnm: directory where the shapefiles will go

spdfinname: Name of the ESRI shapefile written to disk that needs cleaning, without the .shp extension

spdfoutname: Name of shapefile that pprepair should write to, without the .shp extension

Returns:

spdf of cleaned polygons, pointing to a temporary shapefile written to disk

inname <- paste(dirnm, "/", spdfinname, ".shp", sep = "") outname <- paste(dirnm, "/", spdfoutname, ".shp", sep = "") ppcall <- paste("/usr/local/bin/pprepair -i", inname, "-o", outname, "-fix") ctch <- system(ppcall, intern = TRUE) polyfixed <- readOGR(dsn = dirnm, layer = spdfoutname, verbose = FALSE) polyfixed@proj4string <- crs return(polyfixed) }

createCleanTempPolyfromWKT <- function(geom.tab, crs) {

Function for reading in a spatial geometry from PostGIS and creating a temporary shapefile out of it

Args:

geom.tab: Dataframe with geometry and identifiers in it. Identifier must be 1st column, geometries 2nd col

crs: Coordinate reference system

Returns:

A SpatialPolygonsDataFrame

Notes:

Uses callPprepair function to clean up read in polygons and write them to a temporary location

polys <- tst <- sapply(1:nrow(geom.tab), function(x) { poly <- as(readWKT(geom.tab[x, 2], p4s = crs), "SpatialPolygonsDataFrame") poly@data$ID <- geom.tab[x, 1] newid <- paste(x) poly <- spChFIDs(poly, newid) return(poly) }) polyspdf <- do.call("rbind", polys) polys.count <- nrow(polyspdf) td <- tempdir() tmpnmin <- strsplit(tempfile("poly", tmpdir = ""), "/")[[1]][2] tmpnmout <- strsplit(tempfile("poly", tmpdir = ""), "/")[[1]][2] writeOGR(polyspdf, dsn = td, layer = tmpnmin, driver = "ESRI Shapefile") polyfixed <- callPprepair(td, spdfinname = tmpnmin, spdfoutname = tmpnmout, crs = polyspdf@proj4string) valid.string <- as.character(gIsValid(polyfixed)) return(list("polygons" = polyfixed, "polygoncount" = polys.count, "validity" = valid.string)) }

Would you suggest a workflow that uses prepair first to fix individual polygons, followed by pprepair to clean up overlaps, etc, or if the validity issue we found something that pprepair is likely to be able to clean up?

Following from this, is there a way to use pprepair to read WKT directly, as I see prepair does?

Thanks!

[hugoledoux]

At this moment I would suggest that workflow yes. First it's less likely to crash (although I was surprised, it very rarely happens) and second, prepair gives you more control over the repair. You should perhaps check the "--isr" option that performs snap rounding and permits you to eliminate "spikes" (there were a few, probably unwanted, in the dataset you linked to). "--minarea" could also be useful.

We plan to integrate the prepair code in pprepair in the near future, not sure when we'll have the time though.

pprepair doesn't read directly WKT, although a GeoJSON file should work, or anything OGR.

[kenohori]

I agree with Hugo here. At this moment it's probably a good idea to call prepair for individual polygons before using pprepair. We do have the intention of incorporating the code that solves this problem (I think it's related to #12) into pprepair, but we need to put some more thought into the consequences before doing so.

[ldemaz]

Great, thanks for your responses on that. I will look at revising our code to do this. Thanks for the rapid responses...

[kenohori]

This will definitely be solved by #12, so I'll close this issue.