Hole lies outside shell and arcpy

[irfanlatif]

hi there,

In my workflow, I am converting our mapset in ESRI (ArcGIS 10.7) FileGeodb (topological clean and no geometry error) to Postgis.

In PostGIS I am getting a geometry validation error  "Hole lies outside shell" and I am thinking how to find this error using arcpy.

Here is polygon where I am getting this error.

MULTIPOLYGON(((18.580058 54.387989,18.580688 54.387949,18.581455 54.38745,18.58103 54.387236,18.581292 54.387046,18.581613 54.387193,18.581928 54.387025,18.581818 54.386618,18.581261 54.386702,18.580977 54.38657,18.581072 54.386234,18.580021 54.385724,18.579978 54.386042,18.580181 54.386148,18.579153 54.386818,18.579249 54.386861,18.578992 54.387036,18.578708 54.386886,18.578114 54.386992,18.57806 54.387864,18.578553 54.388089,18.580058 54.387989)),((18.580058 54.387989,18.58009 54.388582,18.580315 54.388695,18.581319 54.388631,18.580058 54.387989),(18.580058 54.387989,18.579111 54.387517,18.579764 54.38739,18.580048 54.387528,18.580058 54.387989))) 

I know st_makevalid() will correct this problem. What I want, is to correct my source data which is in FileGeodb.

any suggestions?

[DanPatterson_Retired]

I had to add an enclosing set of brackets to make the two multipart shapes represent 1 shape

z1 = [
    [[[18.580058, 54.387989], [18.580688, 54.387949], [18.581455, 54.38745],
      [18.58103, 54.387236], [18.581292, 54.387046], [18.581613, 54.387193],
      [18.581928, 54.387025], [18.581818, 54.386618], [18.581261, 54.386702],
      [18.580977, 54.38657], [18.581072, 54.386234], [18.580021, 54.385724],
      [18.579978, 54.386042], [18.580181, 54.386148], [18.579153, 54.386818],
      [18.579249, 54.386861], [18.578992, 54.387036], [18.578708, 54.386886],
      [18.578114, 54.386992], [18.57806, 54.387864], [18.578553, 54.388089],
      [18.580058, 54.387989]]], 
    [[[18.580058, 54.387989], [18.58009, 54.388582], [18.580315, 54.388695],
      [18.581319, 54.388631], [18.580058, 54.387989]], 
      [[18.580058, 54.387989], [18.579111, 54.387517], [18.579764, 54.38739],
      [18.580048, 54.387528], [18.580058, 54.387989]]]
      ]‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The two blocks ( lines 2 - 9 and lines 10 - 13) 

weren't closed in  the [ on line 1 and on line 14

Sorry for converting tuples of tuples to lists of lists... just easier to visualize for me

[BruceHarold]

Irfan is correct that ArcGIS does not detect that the geometry is OGC-invalid, here is how it looks:

Therefore ArcGIS's Repair Geometry is not an option.  I'll ask the team to look at that.

In Pro we have the osgeo and ogr modules, somewhere in there is the MakeValid() function.

With Data Interoperability/FME you can de-aggregate the offending bowtie feature which will then need more work.

[JoshuaBixby]

Bruce Harold‌, I am confused by your statement:

ArcGIS does not detect that the geomety is OGC-invalid

If I run the geometry through Check Geometry using ArcGIS 10.7.1, it flags the polygon as self-intersecting, which is OGC-invalid.  Unfortunately Repair Geometry fails to fix it, but detecting an issue and being able to fix it are separate issues.

UPDATE:  As mentioned in another response, I created the polygon in a spatial reference that rounded the original coordinates, thus creating the self-intersection (or maybe one could argue showed the self-intersection).  When using a geographic spatial reference with higher precision, the Check Geometry tool does see the polygon as invalid.

[BruceHarold]

I'm confused too!  I created a feature class from the WKT using Workbench and it made the below

Here it is as WKT

Aside from the precision this looks like what went in, and 10.7.1 Check Geometry reports nothing for me.  Like you say it needs manual work.

[JoshuaBixby]

I figured out the issue.  I wasn't paying attention to the values themselves and squeezed the polygon into a projected coordinate space, which did some rounding and created the self-intersection.  Oddly, the WKT back out from the rounded polygon wasn't rounded itself, but maybe Check Geometry is applying the XY tolerance and not the ArcPy WKT constructor.

[JoshuaBixby]

You have a true multipart geometry (since Esri exports all polygon WKT as MULTIPOLYGON, the 'MULTI' doesn't mean the polygon has multiple parts), and one of the two parts in invalid.

Specifically, the first part is valid and the second part is invalid because it self-intersects.  In the screenshot below, red vertices are from the first part while blue vertices are from the second part.

If I run Check Geometry, it catches the invalid polygon and states self-intersection as the problem.  Once a polygon is invalid, there is no silver bullet for fixing it.  Sure, there are methods that can clean it up, but assumptions are made along the way by the tools, and those assumptions might not work for the user.  Since the shape of the polygon means something specific to the record in the data set, cleaning up an invalid manually is the only way to ensure the updated polygon accurately represents what it is supposed to represent.

[irfanlatif]

Which software you are using for Check Geometry?

I am also runing check geometry (ArcGIS 10.7) but not catching any thing.

[KoryKramer]

ArcGIS Pro's Check Geometry has an OGC validation method.

Check Geometry—Data Management toolbox | ArcGIS Desktop 

Would that work in this case?

[BruceHarold]

Well Pro with the OGC option eports 'non-simple' as an issue, which is better than using 10.7.1, but not as good as the GeometryValidator in Workbench which provides issue details and  location:

GeometryTrait(string) : `_issues{0}.issue_found' has value `Fails OGC Valid'

GeometryTrait(64 bit real) : `_issues{0}.location_sample.x' has value `18.579111'

GeometryTrait(64 bit real) : `_issues{0}.location_sample.y' has value `54.387517'

GeometryTrait(string) : `_issues{0}.repair_state' has value `Not Fully Repaired'

GeometryTrait(string) : `_issues{0}.supplementary_info' has value `Hole Outside Shell'

GeometryTrait(string) : `_issues{0}.supplementary_info{0}.details' has value `Hole Outside Shell'

GeometryTrait(64 bit real) : `_issues{0}.supplementary_info{0}.location.x' has value `18.579111'

GeometryTrait(64 bit real) : `_issues{0}.supplementary_info{0}.location.y' has value `54.387517'

🙂