Union and multi-part polygons in ArcGIS API for Python

Question

I have a requirement to process some features from a feature layer in ArcGIS Online to merge them together based on unique combinations of attributes and write the output to another feature layer. We want to do this in a scheduled notebook, but have been having some issues with the union operation when it comes to creating multi-part polygons.

When I union two polygons that are disjoint (not intersecting) I expect the union operation to create a single multi-part polygon, as it does in desktop and arcpy. However, when we do this in a standard ArcGIS Online notebook, the result is not a multi-part polygon, but just the first of the two polygons used as input.

The below is a test using an ArcGIS Online notebook. When the input polygons are intersecting, the union works as expected. However, when they are disjoint, it doesn't work as expected:

I have tested the exact same notebook code within ArcGIS Pro and get a multi-part feature when the input polygons are disjoint (non intersecting), which is expected:

I have attached the test python code.

I know there are some differences in the way the API works in ArcGIS Online vs ArcGIS Pro, as arcpy is not used in the standard Online notebooks. However, the documentation for the API seems to suggest union is supported and will use the shapely library when arcpy isn't available.

It seems like a fairly basic thing to want to do, so does anyone know how to get this union working, so that unioning two disjoint polygons results in a multipart polygon, in an ArcGIS Online notebook without arcpy?

Is there something specific I need to do to use shapely?

test-union-polys.zip

SaffiaHossainzadeh1 · Answer

I have found another instance of what might be a similar issue. When using the intersect function with two geometry objects, the ArcGIS Python API in ArcGIS Online only returns part of the full intersection. It seems to be a case of what should be a multi-part polygon.  Here are the two intersecting geometries:

Here is the intersection that is returned:

The full area of intersection should be about 2,600 acres, according to what I'd get if I performed the operation in ArcGIS Pro. Also, I'm not sure why the result I'm getting from the intersect() function is negative. The geometry doesn't have holes in it to create an internal/external ring. It does seem to be multi-part, but it's not the total area that should be returned as the intersection. So I don't know what the intersect() function returned.

I have also submitted a support ticket for my issue.

JohnFannon · Answer

I have now raised a support request for this issue.

After further investigations and experimentation, it seems that union will use the shapely union function when arcpy is not available, which also involves a conversion behind the scenes between Esri and shapely geometry. It seems that after unioning, the resulting polygon has both incorrect ring orientation (exterior rings should be clockwise orientation) for Esri and the srid is also set to 4326 regardless of what the input srids were.

I've had to come up with a fix to re-orient the polygons following union, which seems to work in most cases and the geometry is then valid and can be added to a feature layer successfully.

However, there also seem to be multiple issues with the shapely union function, especially when using multi-part polygons as input, and this still results in invalid geometries in some cases.

Overall I'm not that impressed with what should really be a straightforward spatial operation in the ArcGIS API for python. I should not be having to write code to fix the resulting geometry issues.