M dimension dropped from Intersection of LineStringZM and MultiLineStringZM

[jlhumber]

Hi geos team,

In running tests we see that the intersection between a LineStringZM or MultiLineStringZM with one or more LineStringZM or MultiLIneStringZM drops the M dimension. This simple example mimics the what we see:

from shapely import from_wkt, MultiLineString
line_a = from_wkt('LineString (0 0 0 NaN, 10 0 123 NaN)')
# <LINESTRING ZM (0 0 0 NaN, 10 0 123 NaN)>
line_b = MultiLineString([from_wkt('LineString (4 -5 999 NaN, 5 5 456 NaN, 6 -6 678 NaN)')])
# <MULTILINESTRING ZM ((4 -5 999 NaN, 5 5 456 NaN, 6 -6 678 NaN))>
line_a.intersection(line_b)
# <MULTIPOINT Z ((5.455 0 312), (4.5 0 391.425))>

We think this is specific to NaN values in M but cannot confirm this 100% but are almost certain this is the situation.

ref shapely/shapely#2408

[pramsey]

There is logic that will interpret NaNs in coordinates as missing coordinates, so I am not surprised to see stripping when you feed them in. I'd say at this point you are outside the contract of the library and this is a garbage-in/garbage-out situation.

[jlhumber]

@pramsey thanks for sharing this, now that I understand the behavior I can work with this, feel free to close.

[sgillies]

@pramsey somewhat related question: what value of M should one expect at such an intersection point? The average of the two values interpolated to the intersection point? The difference? Something else? Since M isn't a dimension in the same space as X, Y, and Z, what does intersection even mean for it?

[pramsey]

This is very much angels and pins territory, because each implementor is going to have their own business needs and the rules for Z's and M's might differ from person to person and application to application. Trying to enforce a "leftmost argument wins" rule might be a "reasonable" default for points with conflicting ZM, along with an "interpolate between nearest" rule for introduced points. This reminds me of EMPTY, in that it implicates a lot of work for a very narrow wedge of cases, some of which are quite theoretical.

[jlhumber]

Leftmost seems reasonable and is akin to what we see in commercial applications that handle measures.

Using result = line_a.intersection(line_b) as the example in shapely / geos speak but relating to what is observed in a commercial application:

• line_a is LineString M and line_b is LineString M then result will be Point M with M values from line_a
• line_a is LineString and line_b is LineString M then result will be Point M with M that are NaN
• line_a is LineString M and line_b is LineString then result will be Point M with M values from line_a

[dbaston]

For lines, M values are intended to be interpolated/assigned using the same fairly confusing rules as for Z. IIRC the leftmost argument wins if the constructed point is an endpoint of an input geometry, otherwise it is interpolated.

I recently added comments to the tests at https://github.com/libgeos/geos/blob/main/tests/unit/algorithm/RobustLineIntersectorZTest.cpp#L117 to describe this behavior.

[sgillies]

@dbaston those comments are very helpful! Do the same rules apply to other operations? Difference, symmetric difference, union, etc?

[dbaston]

Those rules are applied whenever vertices are generated from the intersection of lines. There are other cases where Z values are added that do not involve the intersection of lines, for example when calculating the intersection of a XY line located within a XYZ polygon:

dan=# SELECT ST_AsText(ST_Intersection('POLYGON Z ((0 0 5, 10 0 10, 10 10 15, 0 10 20, 0 0 5))', 'LINESTRING (3 8, 9 4)'));
           st_astext            
--------------------------------
 LINESTRING Z (3 8 20,9 4 12.5)

This process is not applied to M values:

dan=# SELECT ST_AsText(ST_Intersection('POLYGON M ((0 0 5, 10 0 10, 10 10 15, 0 10 20, 0 0 5))', 'LINESTRING (3 8, 9 4)'));
      st_astext      
---------------------
 LINESTRING(3 8,9 4)