Vector file with Polygon Z with 2 coords causing errors

[latot]

Hi, I was working with some files today, and I got something really weird, the vector file can be in gpkg or shp and we can get the same result, in particular I was getting errors from R/sf with ParseException: Unexpected EOF parsing WKB when doing some operations, so I start checking the geometries, then I notice that there is geometries that has Polygons with Z vals, but when we check the WKT version...:

292401.9 5888301

There is only two coords, I try loading the WKT of 3° geometry in postgis getting:

SQL Error [XX000]: ERROR: can not mix dimensionality in a geometry
  Hint: "...�" <-- parse error at position 29490 within geometry
  Position: 8

.... No idea how why this is happening.

Thx!

test.zip

[jratike80]

That geopackage seems to have mostly invalid geometries with unclosed rings. How have you managed to create them?
You can test yourself with this command:

ogrinfo test.gpkg -sql "select st_isvalidreason(geom) from test"
INFO: Open of `test.gpkg'
      using driver `GPKG' successful.

Layer name: SELECT
Geometry: None
Feature Count: 34
Layer SRS WKT:
(unknown)
st_isvalidreason(geom): String (0.0)
OGRFeature(SELECT):0
  st_isvalidreason(geom) (String) = Valid Geometry

OGRFeature(SELECT):1
  st_isvalidreason(geom) (String) = Valid Geometry

OGRFeature(SELECT):2
  st_isvalidreason(geom) (String) = Invalid: Unclosed Rings were detected

OGRFeature(SELECT):3
  st_isvalidreason(geom) (String) = Invalid: Unclosed Rings were detected

The faulty rings seem to have all zero coordinates at the end of the geometry. Here is a part of such geometry as a sample:

`...5830947.435195 148608.780061 5830947.059777,
0 0 0, 0 0 0, 0 0 0, 0 0 0, 0 0 0,
0 0 0, 0 0 0, 0 0 0, 0 0 0, 0 0 0,
0 0 0, 0 0 0, 0 0 0, 0 0 0, 0 0 0,
...

[latot]

which geom did you pick? the 3° one only has two coords as poly z, and now that one has the 3 coords, but that coords should be right, append a lot of 0, 0, 0 is still a valid geometry, which could be simplified.

I'm asking here where are the original files, and how they was merged.

[jratike80]

The last vertex of a polygon must be the same as the first vertex. Here is one small ring that ends with 0,0,0 but does not start with 0,0,0.

((292401.858075 5888301.414703 292405.568241,
 5888161.333336 292398.867831 5888042.181091,
 292394.398099 5888035.760441 292357.80551,
 5887983.195835 292287.841197 5887995.572598,
 292264.884107 5888004.165759 292193.795034,
 5888030.774766 292174.549904 5888034.785231,
 292274.745 5888215.8907 292401.621262,
 5888343.527775 292401.858075 5888301.414703,
 0 0 0, 0 0 0, 0 0 0, 0 0 0))

[latot]

how that can happens.... gdal does not allow to read that things, so we should not be able to have this cases.... make_valid also does not fix them.

[rouault]

. gdal does not allow to read that things

what makes you think so ? GDAL is very lenient on reading. It can possibly read non-closed polygons if the input datasource expose them.

[latot]

mmm, you are right, I didn't specify enough good, we can read from a file, but we can't read from WKT.

[jratike80]

Most geometries in the geopackage are bad. Do you have backup data? If you have it would be waste of time to try to use that geopackage.

[latot]

Testing, there is some weird things:

Make Valid does not works directly on geometries.
Drop Z and then Make Valid the geometries fixes it.

Maybe this is more a issue in GEOS?

[elpaso]

I am looking at this issue but I do not understand what is the issue or what is the expected behavior, I wrote this test (which passes) and this seems the expected behavior.

def test_ogr_gpkg_2D_geometry_in_non_z_layer(tmp_vsimem, tmp_path):

    ds = gdal.GetDriverByName("GPKG").Create(
        tmp_vsimem / "tmp.gpkg", 0, 0, 0, gdal.GDT_Unknown
    )

    lyr = ds.CreateLayer("foo", geom_type=ogr.wkbPolygon25D)

    # polygon ok
    feat = ogr.Feature(lyr.GetLayerDefn())
    assert 0 == feat.SetGeometryDirectly(ogr.CreateGeometryFromWkt("POLYGON Z ((0 0 0, 1 0 1, 1 1 2, 0 1 1, 0 0 0))"))
    with gdal.quiet_errors():
        assert 0 == lyr.CreateFeature(feat)


    # missing Z
    feat = ogr.Feature(lyr.GetLayerDefn())
    feat.SetGeometryDirectly(ogr.CreateGeometryFromWkt("POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))"))
    with gdal.quiet_errors():
        lyr.CreateFeature(feat)


    # open polygon
    feat = ogr.Feature(lyr.GetLayerDefn())
    feat.SetGeometryDirectly(ogr.CreateGeometryFromWkt("POLYGON ((0 0, 1 0, 1 1, 0 1))"))
    with gdal.quiet_errors():
        lyr.CreateFeature(feat)


    ds = None

    ds = ogr.Open(tmp_vsimem / "tmp.gpkg")
    lyr = ds.GetLayer(0)
    assert lyr.GetGeomType() == ogr.wkbPolygon25D

    feat = lyr.GetNextFeature()

    feat = lyr.GetFeature(1)
    g = feat.GetGeometryRef()
    gvalid = g.MakeValid()
    assert gvalid is not None
    assert gvalid.GetGeometryType() == ogr.wkbPolygon25D

    # geometry is 2D
    feat = lyr.GetFeature(2)
    g = feat.GetGeometryRef()
    gvalid = g.MakeValid()
    assert gvalid is not None
    assert gvalid.GetGeometryType() == ogr.wkbPolygon

    # open polygon: is not valid!
    feat = lyr.GetFeature(3)
    g = feat.GetGeometryRef()
    gvalid = g.MakeValid()
    assert gvalid is  None

[latot]

Hi, well the expected behavior is getting back be able to "read" the data, here is the weird part, we can read invalid polygons, but when you try to work with them we will get ParseException: Unexpected EOF parsing WKB.

So detect and fix it from GDAL could be a good way to do it.

About https://github.com/OSGeo/gdal/pull/9355/files I think would be better to instead promote to 3d, if all the points only have two coords change the polygon to 2D, because promote to 3d would also take the assumptions about the value of Z. Use Nan as Z value can also cause some confusion in some calcs (?

[jratike80]

Have you already found the reason why such corrupted WKB exists? Did you use some buggy software?

[latot]

it comes from 3 sources that I don't know, also they were merged to one file probably with Arcgis..........