topic Re: Hole lies outside shell and arcpy in Python Questions

Hole lies outside shell and arcpy

irfanlatif — Thu, 21 Nov 2019 10:55:57 GMT

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?

Re: Hole lies outside shell and arcpy

DanPatterson_Retired — Sat, 11 Dec 2021 16:59:47 GMT

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

Re: Hole lies outside shell and arcpy

BruceHarold — Thu, 21 Nov 2019 14:36:39 GMT

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.

Re: Hole lies outside shell and arcpy

JoshuaBixby — Thu, 21 Nov 2019 14:42:49 GMT

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.

Re: Hole lies outside shell and arcpy

JoshuaBixby — Thu, 21 Nov 2019 14:50:06 GMT

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.

Re: Hole lies outside shell and arcpy

KoryKramer — Thu, 21 Nov 2019 14:51:48 GMT

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

Check Geometry—Data Management toolbox | ArcGIS Desktop

Would that work in this case?

Re: Hole lies outside shell and arcpy

BruceHarold — Thu, 21 Nov 2019 15:05:05 GMT

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.

Re: Hole lies outside shell and arcpy

BruceHarold — Thu, 21 Nov 2019 15:14:32 GMT

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'

🙂

Re: Hole lies outside shell and arcpy

JohannesBierer — Thu, 21 Nov 2019 15:31:19 GMT

I ofter can't find self intersections with the Check Geometry Tool, it is a bit unreliable especially with very small mistakes. I can offer several shapes with self intersections not found by Check Geometry Tool (10.4.1) - unfortunately I don't have pro.

Re: Hole lies outside shell and arcpy

JoshuaBixby — Thu, 21 Nov 2019 15:50:41 GMT

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.

Re: Hole lies outside shell and arcpy

irfanlatif — Thu, 21 Nov 2019 16:13:06 GMT

Which software you are using for Check Geometry?

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

Re: Hole lies outside shell and arcpy

JoshuaBixby — Thu, 21 Nov 2019 17:39:41 GMT

As I mentioned above, I put the coordinates into a projected spatial reference with lower tolerances, which appears to have done rounding and allowed Check Geometry to find the self-intersection.

The example geometry you provide raises several issues with Esri tools and underlying geometry model, so in that regard it is a good example, but that doesn't help you easily here since there are no easy workaround for the issues your example polygon raises.

Re: Hole lies outside shell and arcpy

irfanlatif — Sat, 11 Dec 2021 16:59:50 GMT

Thanks, everyone for your input.

As Joshua mention in ESRI environment there is no easy solution. So I have to go outside for help.

I build a small function using ogr to validate geometry. I have used buffer0 which is not a good option. If you have incorrect ring odder, you may end up with duplicate geometry. In my case it works fine as my base data is validated first with ESRI tools and then with ogr.

def ogr_Validate(fc):    
    """
    Used to validate input shapefile geometries with ogr
    INPUTS:
            fc(req) = path to fc            
            example. ogr_Validate(fc)
    IMP: Used buffer0 which in case of incorrect ring oder can produce duplicate geometries.    
    """

    inFc = ogr.Open(fc, update =1)
    inFc_lyr = inFc.GetLayer()
    
    for feature in inFc_lyr:
        geom = feature.GetGeometryRef()            
        if not geom.IsValid():
            print "Cleaning feature:      "+str(feature.GetField(lst_Fields(fc)[0]))        
            feature.SetGeometry(geom.Buffer(0))            
            inFc_lyr.SetFeature(feature) 
            assert feature.GetGeometryRef().IsValid() 
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Reference:

python - Invalid geometries made valid don't remain valid? - Geographic Information Systems Stack Exchange

Re: Hole lies outside shell and arcpy

KoryKramer — Wed, 27 Nov 2019 15:22:10 GMT

The geometry team has reviewed this issue and wanted to convey the following:

The geometry in question is OGC Multipolygon type. If you look carefully it has three rings (with the two rings being in the second part).

While Esri polygon geometry type can represent both OGC Multipolygon and OGC Polygon types, it represents both as a collection of rings with exterior rings being clockwise and holes being counterclockwise. When the Multipolygon in question is converted to Esri format it becomes a three ring polygon, and it is valid as far as Esri simple geometry is concerned.

The OGC specification of simple features is more strict than Esri’s. To help with interoperability we’ve added the OGC option in the Check geometry and Repair Geometry tools. When the OGC option is used, we check additional rules of the OGC spec and ensure that the polygon rings are sorted such that holes follow exterior rings.

What looks like the bowtie highlighted with blue color is actually two rings touching at a vertex.

Here is how to fix the original polygon to be valid for OGC spec:

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) move it here),((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)))

It only requires moving the last ring int the first POLYGON part.

Re: Hole lies outside shell and arcpy

JoshuaBixby — Mon, 02 Dec 2019 03:56:08 GMT

Kory Kramer‌, see my https://community.esri.com/message/894376-check-geometry-fails-in-shared-origin-edge-case thread that highlights a defect in Check Geometry that fails to identify invalid geometries in certain edge cases.