libgeos / geos

Geometry Engine, Open Source
https://libgeos.org
GNU Lesser General Public License v2.1
1.17k stars 354 forks source link

Dimension handling with setPrecision / reducePrecision #1152

Open mwtoews opened 3 weeks ago

mwtoews commented 3 weeks ago

I see a few inconsistencies with GEOSGeom_setPrecision. For example:

  1. $ ./bin/geosop -a "LINEARRING (0 0, 0.1 0, 0.1 0.1, 0 0.1, 0 0)" reducePrecision 1
    LINEARRING Z EMPTY
  2. $ ./bin/geosop -a "LINEARRING (0 0, 0.1 0, 0.1 0.1, 0 0.1, 0 0)" reducePrecisionKeepCollapsed 1
    LINESTRING (0 0, 0 0, 0 0)
  3. $ ./bin/geosop -a "LINEARRING (0 0, 0.1 0, 0.1 0.1, 0 0.1, 0 0)" reducePrecisionPointwise 1
    LINEARRING (0 0, 0 0, 0 0, 0 0, 0 0)

    Here are my guesses for better answers:

  4. Should be LINEARRING EMPTY (without Z dimension)
  5. Should not be LINESTRING, but perhaps be LINEARRING (0 0, 0 0, 0 0)
  6. Probably fine?

More on case 1, GEOSGeom_setPrecision test<10> doesn't fail because ensure_geometry_equals ignores Z/M dims, thus LINEARRING Z EMPTY == LINEARRING EMPTY. There is a similar issue with GEOSGeom_transformXY test<10> where POINT Z EMPTY is returned and compared equal to POINT EMPTY. I'd like ensure_geometry_equals to check if (GEOSHasZ(g1) != GEOSHasZ(g2)) || (GEOSHasM(g1) != GEOSHasM(g2)).

The issue is probably related to GeometryTransformer / GeometryFactory::createLinearRing().

See also trac #1135

mwtoews commented 3 weeks ago

With case 2, I see GeometryTransformer::transformLinearRing switch to a LineString type to ensure a valid LinearRing, since it has three coordinates.

mwtoews commented 3 weeks ago

More generally, the issue is that setPrecision / reducePrecision doesn't properly handle dimensions. A few examples:

for geom_type in POINT LINESTRING LINEARRING POLYGON; do
    geom="$geom_type M EMPTY";
    echo "$geom";
    echo -n "  reducePrecision: ";
    ./bin/geosop -a "$geom" reducePrecision 1;
    echo -n "  reducePrecisionKeepCollapsed: ";
    ./bin/geosop -a "$geom" reducePrecisionKeepCollapsed 1;
    echo -n "  reducePrecisionPointwise: ";
    ./bin/geosop -a "$geom" reducePrecisionPointwise 1;
done

returns all wrong results

POINT M EMPTY
  reducePrecision: POINT EMPTY
  reducePrecisionKeepCollapsed: POINT EMPTY
  reducePrecisionPointwise: POINT Z EMPTY
LINESTRING M EMPTY
  reducePrecision: LINESTRING EMPTY
  reducePrecisionKeepCollapsed: LINESTRING EMPTY
  reducePrecisionPointwise: LINESTRING Z EMPTY
LINEARRING M EMPTY
  reducePrecision: LINEARRING Z EMPTY
  reducePrecisionKeepCollapsed: LINEARRING Z EMPTY
  reducePrecisionPointwise: LINEARRING Z EMPTY
POLYGON M EMPTY
  reducePrecision: POLYGON EMPTY
  reducePrecisionKeepCollapsed: POLYGON EMPTY
  reducePrecisionPointwise: LINEARRING Z EMPTY

Or with non-empty inputs:

for geom in "POINT M (1 2 3)" "LINESTRING M (0 0 1, 0.1 0 1, 0.1 0.1 1, 0 0.1 1, 0 0 1)" "LINEARRING M (0 0 1, 0.1 0 1, 0.1 0.1 1, 0 0.1 1, 0 0 1)" "POLYGON M ((0 0 1, 0.1 0 1, 0.1 0.1 1, 0 0.1 1, 0 0 1))"; do
    echo "$geom";
    echo -n "  reducePrecision: ";
    ./bin/geosop -a "$geom" reducePrecision 1;
    echo -n "  reducePrecisionKeepCollapsed: ";
    ./bin/geosop -a "$geom" reducePrecisionKeepCollapsed 1;
    echo -n "  reducePrecisionPointwise: ";
    ./bin/geosop -a "$geom" reducePrecisionPointwise 1;
done

the results have some OK and some with NaN and others with lost dims:

POINT M (1 2 3)
  reducePrecision: POINT M (1 2 3)
  reducePrecisionKeepCollapsed: POINT M (1 2 3)
  reducePrecisionPointwise: POINT (1 2)
LINESTRING M (0 0 1, 0.1 0 1, 0.1 0.1 1, 0 0.1 1, 0 0 1)
  reducePrecision: LINESTRING EMPTY
  reducePrecisionKeepCollapsed: LINESTRING M (0 0 1, 0 0 NaN)
  reducePrecisionPointwise: LINESTRING (0 0, 0 0, 0 0, 0 0, 0 0)
LINEARRING M (0 0 1, 0.1 0 1, 0.1 0.1 1, 0 0.1 1, 0 0 1)
  reducePrecision: LINEARRING Z EMPTY
  reducePrecisionKeepCollapsed: LINESTRING M (0 0 1, 0 0 NaN, 0 0 NaN)
  reducePrecisionPointwise: LINEARRING (0 0, 0 0, 0 0, 0 0, 0 0)
POLYGON M ((0 0 1, 0.1 0 1, 0.1 0.1 1, 0 0.1 1, 0 0 1))
  reducePrecision: POLYGON EMPTY
  reducePrecisionKeepCollapsed: POLYGON EMPTY
  reducePrecisionPointwise: POLYGON ((0 0, 0 0, 0 0, 0 0, 0 0))