yeesian / ArchGDAL.jl

A high level API for GDAL - Geospatial Data Abstraction Library
https://yeesian.github.io/ArchGDAL.jl/stable/
Other
142 stars 26 forks source link

Can't create valid linearring #325

Open henrik-wolf opened 2 years ago

henrik-wolf commented 2 years ago

While working through the Documentation I was not able to create a linearring which would pass the validity test. MWE:

julia> using ArchGDAL

julia> ring1 = ArchGDAL.createlinearring([(0.,0.), (0.,1.), (1.,1.)])
Geometry: LINEARRING (0 0,0 1,1 1)

julia> ring2 = ArchGDAL.createlinearring([(0.,0.), (0.,1.), (1.,1.), (0.,0.)])
Geometry: LINEARRING (0 0,0 1,1 1,0 0)

julia> ArchGDAL.isvalid(ring1)
false

julia> ArchGDAL.isvalid(ring2)
false

julia> ArchGDAL.closerings!(ring1)
Geometry: LINEARRING (0 0,0 1,1 1,0 0)

julia> ArchGDAL.isvalid(ring1)
false

Interestingly, Polygons can be valid, if the contained linearrings have the same first and last point, even if the rings themselves are not valid. (and if the other conditions for polygons are met)

visr commented 2 years ago

Interesting. Using LibGEOS.jl this works as expected:

julia> using LibGEOS

julia> ring1 = readgeom("LINEARRING (0 0,0 1,1 1)")
ERROR: GEOSError
        IllegalArgumentException: Points of LinearRing do not form a closed linestring

julia> ring2 = readgeom("LINEARRING (0 0,0 1,1 1,0 0)")
LinearRing(Ptr{Nothing} @0x0000000075ee4360)

julia> isValid(ring2)
true

In GDAL.jl I guess this would reproduce it, but it gives an error code back for unsupported geometry.

using GDAL
ring2 = Ref(C_NULL)
wkt = "LINEARRING (0 0,0 1,1 1,0 0)"
GDAL.ogr_g_createfromwkt([wkt], C_NULL, ring2)
# returns OGRERR_UNSUPPORTED_GEOMETRY_TYPE = 3
# GDAL.ogr_g_isvalid(ring2[])  # doesn't work since ring2 is still C_NULL

This is probably related to GDAL itself and how it handles this geometry type, see also #314.

I added https://github.com/JuliaGeo/GDAL.jl/pull/141 to double check that we build GDAL with GEOS support.

yeesian commented 2 years ago

Hmm yeah I checked the implementation in https://github.com/OSGeo/gdal/blob/82401cb09515840dfde480ac4aa10cac7286983a/ogr/ogrgeometry.cpp#L2191-L2213 and don't know what else to try here. If we might be able to translate it into the corresponding GDAL function calls, we can file an issue with https://github.com/OSGeo/gdal/issues

yeesian commented 2 years ago

It seems it might have been reported in the past but never addressed: https://trac.osgeo.org/gdal/ticket/5006. To fix it, GDAL should support bidirectional translation to/from WKB for LinearRings.

However, I don't think LibGEOS properly supports WKB for LinearRings. To verify this, running the following code will return a linestring:

using LibGEOS
wkbwriter = LibGEOS.WKBWriter(LibGEOS._context)
ring = LibGEOS.readgeom("LINEARRING (0 0,0 1,1 1,0 0)")
wkb = LibGEOS.writegeom(ring, wkbwriter)
LibGEOS.readgeom(wkb)

On the other hand, running the following code will return a linearring:

using LibGEOS
wktwriter = LibGEOS.WKTWriter(LibGEOS._context)
ring = LibGEOS.readgeom("LINEARRING (0 0,0 1,1 1,0 0)")
wkt = LibGEOS.writegeom(ring, wktwriter)
LibGEOS.readgeom(wkt)

I suspect this is due to WKBs not having first-class support as one of the original 7 simple features defined in https://libgeos.org/specifications/wkb/#geometry-types.