Difference with LineString Z and Polygon erroneously affects Z coordinates

[abey79]

While computing the difference between a single-segment LineString Z with a Polygon which happens to share a vertex (but otherwise isn't overlapping), the resulting LineString has one of the Z coordinate altered.

Here is the geometry:

Here is how to reproduce:

import numpy as np
import matplotlib.pyplot as plt
from shapely.wkt import loads, dumps


def main():
    ls = loads(
        "LINESTRING Z (-0.08365735370739519 0.08221315533414053 -1.010732607765682, -0.1533588930799144 0.0653358494958651 -1.010390339846296)"
    )
    p = loads(
        "POLYGON ((-0.2239294216623028 -0.01825216019427908, -0.1515331919718202 -0.03513918728914136, -0.1533588930799144 0.0653358494958651, -0.2239294216623028 -0.01825216019427908))"
    )
    result = ls.difference(p)


    plt.plot(*ls.coords.xy, 'r-')
    plt.plot(*p.boundary.xy, 'g-')
    plt.axis("equal")
    plt.show()

    print(dumps(result))
    print(np.array(result) == np.array(ls))


if __name__ == "__main__":
    main()

...and the result:

LINESTRING Z (-0.0836573537073952 0.0822131553341405 -1.0107326077656820, -0.1533588930799144 0.0653358494958651 -1.0107326077656820)
[[ True  True  True]
 [ True  True False]]

Am I missing something or is this a bug? (macOS 10.14.6, MacPort's Python 3.7.4, Shapely 1.6.4.post2)

[abey79]

I did a quick test with shapely 1.7a2 and MacPort's GEOS 3.7.2 with the same result.

[sgillies]

GEOS doesn't intend to be "Z correct" (for lack of a better phrase) and so Shapely can't help here.

[mwtoews]

I just checked this with GEOS 3.9.1, and this appears to be "fixed" upstream. Perhaps they have a test to check this behaviour, I don't know for sure, but generally don't always assume Z dimension are correctly handled.