topic Re: Issues with shapely and polygons having holes in ArcGIS API for Python Questions

Issues with shapely and polygons having holes

J_R_Matchett — Mon, 05 Mar 2018 16:38:21 GMT

UPDATE: After installing v1.4.2 , I'm no longer experiencing the issue with polygons containing holes; however, multi-part polygons still are not handled properly (only one part is extracted during conversion to a Shapely polygon).

I 'm having some issues with the ArcGIS Python API v1.4 and polygons with holes (interior rings). Attached is a jupyter notebook with an example where I retrieve a US Census block containing a hole. I am unable to intersect the polygon with another polygon (a TopologyException is raised). Additionally, the area property of the polygon appears incorrect, returning the sum of the exterior and interior rings rather than the difference. It appears that the as_shapely method for polygon geometries with holes isn't working properly.

Re: Issues with shapely and polygons having holes

JoshuaBixby — Mon, 05 Mar 2018 21:47:16 GMT

I will have to dive deeper into it, but several online GeoJSON checkers/validators do not like the GeoJSON being put out by that service, which could indicate something is wrong with that specific geometry or the service as a whole. Beyond a specific geometry issue, this could be an issue caused by ring ordering being different in Esri's and other spatial libraries.

Re: Issues with shapely and polygons having holes

JoshuaBixby — Sun, 12 Dec 2021 09:12:27 GMT

Before getting too far into your notebook, the following code to create poly_int by intersecting works just fine for me:

from arcgis.gis import GIS
from arcgis.features import FeatureLayer
from arcgis.geometry import Polygon as Polygon
from shapely.geometry import Polygon as ShapelyPolygon
from shapely.validation import explain_validity


census_layer = FeatureLayer(url=r'https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Tracts_Blocks/MapServer/12')
census_block = census_layer.query("GEOID = '060470010051020'", out_sr=3857)
block_poly = census_block.df.SHAPE[0]
type(block_poly)

clip_poly = Polygon({'rings':
    [[[-13414668, 4483240],
      [-13413951, 4483228],
      [-13413965, 4482677],
      [-13414675, 4482686],
      [-13414668, 4483240]]],
    'spatialReference': {'wkid': 3857}})

gis = GIS()
map = gis.map('Riveria Holiday Mobile Est, Merced, California', zoomlevel=15)

poly_int = clip_poly.intersect(block_poly)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Your notebook error indicates shapely is causing the error, but that is odd since shapely isn't with the code leading up to that point.

Re: Issues with shapely and polygons having holes

J_R_Matchett — Tue, 06 Mar 2018 18:46:15 GMT

Thanks for taking a look at my issue, Joshua. After some further investigating, it appears that my issue is related to running in a python environment without the desktop ArcGIS python packages (e.g. ArcPy) available. I primarily run my code on macOS X, but I tried in a python environment on Windows 10 without ArcPy available and had the same issues, but things worked fine when running in the default ArcGIS Pro python environment. So, it appears proper shapely functionality requires desktop ArcGIS.

Re: Issues with shapely and polygons having holes

JoshuaBixby — Sun, 12 Dec 2021 09:12:30 GMT

The following code is executed when running a spatial method locally (not from a web service) on an ArcGIS Geometry object:

        if HASARCPY:
            if isinstance(second_geometry, Geometry):
                second_geometry = second_geometry.as_arcpy
            return self.as_arcpy.intersect(other=second_geometry,
                                           dimension=dimension)
        elif HASSHAPELY:
            if isinstance(second_geometry, Geometry):
                second_geometry = second_geometry.as_shapely
            return Geometry(self.as_shapely.intersection(
                other=second_geometry).__geo_interface__)
        return None‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

If ArcPy is installed, it is used to run the spatial method while Shapely is used if ArcPy is not installed. I don't have time to dig deeper right now, but I suspect the issue is line #8 because an invalid geometry is being returned to second_geometry.

UPDATE: There is a bug in the __geo_interface__ methods of the API. When ArcPy is installed, the API lets ArcPy handle __geo_interface__ calls, but the API uses native code (the code with a bug) when ArcPy is not installed. The following is from a Python interpreter without access to ArcPy:

Python 3.5.3 |Continuum Analytics, Inc.| (default, May 15 2017, 10:43:23) [MSC v.1900 64 bit (AMD64)] on win32
>>> from arcgis.geometry import Geometry
>>> from pprint import pprint
>>> 
>>> # simple square with a hole in the middle
... json = {
...     "rings": [
...         [
...             [0, 0],
...             [10, 10],
...             [0, 10],
...             [10, 10],
...             [10, 0],
...             [0, 0]
...         ],
...         [
...             [3, 3],
...             [7, 3],
...             [7, 7],
...             [3, 7],
...             [3, 3]
...         ]
...     ],
...     "spatialReference": {
...         "wkid": 3857
...     }
... }
>>> pprint(json)
{'rings': [[[0, 0], [10, 10], [0, 10], [10, 10], [10, 0], [0, 0]],
           [[3, 3], [7, 3], [7, 7], [3, 7], [3, 3]]],
 'spatialReference': {'wkid': 3857}}
>>> 
>>> pprint(Geometry(json).__geo_interface__)
{'coordinates': [[[(0, 0), (10, 10), (0, 10), (10, 10), (10, 0), (0, 0)]],
                 [[(3, 3), (7, 3), (7, 7), (3, 7), (3, 3)]]],
 'type': 'MultiPolygon'}
>>> ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

If you treat round and square brackets the same, which they are for __geo_interface__ purposes, you will notice the API has inserted an extra closing bracket on the outer coordinates. The extra closing bracket converts the geometry from a polygon with a hole in it to a multi-part polygon with overlapping parts, hence the errors you are seeing in Shapely.

Re: Issues with shapely and polygons having holes

J_R_Matchett — Wed, 07 Mar 2018 00:58:05 GMT

I was digging through the API code, and it seems that a Geometry's __geo_interface__ property always returns a "MultiPolygon" type (when ArcPy isn't available), even for non-multi polygons (see lines 330–344 in the geometry/_types.py file). When that result is passed on to shapely's shape function (see lines 21–46 in shapely's geometry/geo.py file), a shapely MultiPolygon will be created, which treats all rings as separate polygons, instead of a shapely Polygon, which treats the first ring as the exterior (aka shell) and the subsequent rings as interiors (aka holes).

Re: Issues with shapely and polygons having holes

JoshuaBixby — Wed, 07 Mar 2018 14:20:06 GMT

The semantics of Esri's geometry model don't line up with, well, just about any other geometry model. I wrote about this in a blog post last summer: /blogs/tilting/2017/06/10/a-case-of-missing-prefixes-esris-geometries .

Whether WKB, WKT, __geo_interface__, GeoJSON, etc...; Esri polygons will be described as multi-polygons, even if they are just a polygon. Since an OGC Simple Feature MULTIPOLYGON is perfectly valid if it only has a single geometry, having polygons described as multi-part polygons mostly works.

The current bug in the ArcGIS API for Python isn't describing a polygon as a multi-part polygon, that is expected with Esri software, the bug has to do with how the multi-part polygon is being structured. As I mentioned earlier, an extra bracket is being inserted which turns a multi-part polygon with a single polygon with a hole in it (what you want) into a multi-part polygon with two polygons with no holes. Since the two polygons overlap each other, the multi-part polygon is invalid, and running spatial methods on it/them fails.

Re: Issues with shapely and polygons having holes

J_R_Matchett — Wed, 07 Mar 2018 18:56:52 GMT

Thanks very much for the insights, Joshua! As a practical solution to the problems I'm encountering, I've wrote an alternative method, as_shapely2, to handle the issues I'm running across in my analyses using US Census Blocks. It's based on the arcgis documentation from Esri that states exterior rings are oriented clockwise and interior rings counter-clockwise. The attached notebook includes the method definition plus some various tests.

EDIT: Updated as_shapely2 method to use shapely's LinearRing for checking ring direction plus other refinements.