I have a phenomenon here and am wondering whether this is a bug or work-as-designed.
We receive geometries from a routing service as polylines (saved as JSON) and we want to unite them geometrically with the union function. (We have no Z-coordinates, so we cant't detect if it is a crossing or a bridge or so on.)
So far, everything is fine in normal line segments. Something strange happens when the line crosses, like here
The polyline, which normally consists of exactly one “part”, is split at the crossing point, as soon as the circle closes (at the exact same crossing point wich has already been added before) Now another part is created. So far everything, so good.
However, the part that was traveled first and was previously part 0 ...
... suddenly becomes Part 1 and the newly inserted part is now Part 0.
As far as I know, new parts are always added at the end. And then to top it all off: If the third part is added after the crossing point, then
the first part is the last part (Part 2),
the second section is the first part (Part 0) and
the last section is the middle part (Part 1)
This means that there is no sequence at all. I don't know any longer, which part was added first and which is the last part of my line.
Here's the drop in replacement function for the stuff below:
def ordered_union(*polylines: arcpy.Polyline) -> arcpy.Polyline:
unioned = polylines[0]
for pl in polylines[1:]:
unioned = unioned.union(pl)
pl_joined = arcpy.Polyline(
arcpy.Array(
[point
for polyline in polylines
for part in polyline
for point in part]), spatial_reference=polyline.spatialReference)
return pl_joined.difference(pl_joined.difference(pl_union))
The docstring for polyline.union is as follows:
def union(self, other):
"""Geometry.union(other)
Constructs the geometry that is the set-theoretic union of the input
geometries.
The two geometries being unioned must be the same shape type. Union
operator
other(Object):
A second geometry."""
...
I think the hint here is set theoretic which betrays that they're likely using a set in the background, which in python is essentially randomly ordered (ordered by hash value, so numbers are always in order, but ordering depends on how the has function is implemented for a specific class). To get around this, it seems like you can create a sequential line (add the arrays together, meaning it is continuous), then diff that with the unioned line, then diff the sequential line with that diff (a bit roundabout I know). This sequence seems to maintain ordering at least for the simple test case:
import arcpy
from pprint import pprint
import json
def flatten(polyline: arcpy.Polyline) -> list[arcpy.Point]:
return [point for part in polyline for point in part]
pl1: arcpy.Polyline = arcpy.Polyline(
arcpy.Array([arcpy.Point(0, 0), arcpy.Point(-1, -1)])
)
pl2: arcpy.Polyline = arcpy.Polyline(
arcpy.Array([arcpy.Point(0, -1), arcpy.Point(-1, 0)])
)
pl_union: arcpy.Polyline = pl1.union(pl2)
print("Raw Union")
pprint(json.loads(pl_union.JSON)['paths'], indent=2)
pl_joined = arcpy.Polyline(arcpy.Array(flatten(pl1) + flatten(pl2)), spatial_reference=pl1.spatialReference)
print(f"\nOrdered Join (with connection segment)")
pprint(json.loads(pl_joined.JSON)['paths'], indent=2)
pl_diff = pl_joined.difference(pl_union)
pl_ord = pl_joined.difference(pl_diff)
print("\nProperly Ordered")
for idx, part in enumerate(pl_ord):
print(f"Part {idx}")
pprint([f"[{point.X}, {point.Y}]" for point in part])
Output:
Raw Union
[ [[0, -1], [-0.5, -0.5]],
[[-0.5, -0.5], [-1, -1]],
[[-0.5, -0.5], [-1, 0]],
[[0, 0], [-0.5, -0.5]]]
Ordered Join (with connection segment)
[[[0, 0], [-1, -1], [0, -1], [-1, 0]]]
Properly Ordered
Part 0
['[0.0, 0.0]', '[-0.5, -0.5]']
Part 1
['[-0.5, -0.5]', '[-1.0, 0.0]']
Part 2
['[0.0, -1.0]', '[-0.5, -0.5]']
Part 3
['[-0.5, -0.5]', '[-1.0, -1.0]']
I'm not 100% sure what's happening under the hood here, but it seems like difference doesn't use the same set language as union:
def difference(self, other):
"""Geometry.difference(other)
Constructs the geometry that is composed only of the region unique to
the
base geometry but not part of the other geometry. The following
illustration shows the results when the red polygon is the source
geometry. Difference operator
other(Object):
A second geometry."""
The same intersection identification happens in the difference operation, but it seems like it keeps the ordering, maybe it's iterating the points instead of returning a point/part set? Either way for now it seems that the union operator has some documented but odd behavior (Python sets confuse people too because after 3.8 they're the only builtin sequence that doesn't preserve order)
