search

yeesian / ArchGDAL.jl

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

Convert Shapefile.PolygonZ to ArchGDAL geometry #416

Open felixcremer opened 4 months ago

felixcremer commented 4 months ago

I am trying to reproject a shapefile geometry with ArchGDAL and I fail at converting the geometry to an ArchGDAL geometry.

Here b is an Shapefile.PolygonZ This conversion works for flat Polygons. 


julia> GI.convert(AG, b)
ERROR: MethodError: no method matching addpoint!(::ArchGDAL.Geometry{ArchGDAL.wkbLineString}, ::Float64, ::Float64, ::Float64, ::Float64)

Closest candidates are:
  addpoint!(::G, ::Real, ::Real, ::Real) where G<:ArchGDAL.AbstractGeometry
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1298
  addpoint!(::G, ::Real, ::Real) where G<:ArchGDAL.AbstractGeometry
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1308

Stacktrace:
 [1] unsafe_createlinearring(coords::Vector{Vector{Float64}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1815
 [2] unsafe_createpolygon(coords::Vector{Vector{Vector{Float64}}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1848
 [3] createmultipolygon(coords::Vector{Vector{Vector{Vector{Float64}}}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1848
 [4] convert(::Type{ArchGDAL.IGeometry}, type::GeoInterface.MultiPolygonTrait, geom::Shapefile.PolygonZ)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/geointerface.jl:352
 [5] convert(package::Module, geom::Shapefile.PolygonZ)
   @ GeoInterface ~/.julia/packages/GeoInterface/tuF7K/src/fallbacks.jl:149
 [6] top-level scope
   @ REPL[129]:1
rafaqz commented 4 months ago

Also pretty bad that it converts to a vector. Instead of calling GI.x, GI.y and GI.z on it directly in the addpoint! constructor. Someone is calling GI.coordinates.

yeesian commented 4 months ago

What are the 4 Float64s in addpoint!(::ArchGDAL.Geometry{ArchGDAL.wkbLineString}, ::Float64, ::Float64, ::Float64, ::Float64) -- x, y, z and?

I was wondering if it might be "M", but it's a PolygonZ, so is it i (which is also weird since it's a Float)? But if it is i, that seems to map to setpoint! (below) rather than addpoint!: https://github.com/yeesian/ArchGDAL.jl/blob/dc19553f4d9e9eb03dca60a6467d2fae0a659f2a/src/ogr/geometry.jl#L1263-L1272

yeesian commented 4 months ago

Someone is calling GI.coordinates.

Is it this line? https://github.com/yeesian/ArchGDAL.jl/blob/dc19553f4d9e9eb03dca60a6467d2fae0a659f2a/src/geointerface.jl#L352

yeesian commented 4 months ago

I am trying to reproject a shapefile geometry with ArchGDAL and I fail at converting the geometry to an ArchGDAL geometry.

Here b is an Shapefile.PolygonZ This conversion works for flat Polygons.

julia> GI.convert(AG, b)
ERROR: MethodError: no method matching addpoint!(::ArchGDAL.Geometry{ArchGDAL.wkbLineString}, ::Float64, ::Float64, ::Float64, ::Float64)

Closest candidates are:
  addpoint!(::G, ::Real, ::Real, ::Real) where G<:ArchGDAL.AbstractGeometry
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1298
  addpoint!(::G, ::Real, ::Real) where G<:ArchGDAL.AbstractGeometry
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1308

Stacktrace:
 [1] unsafe_createlinearring(coords::Vector{Vector{Float64}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1815
 [2] unsafe_createpolygon(coords::Vector{Vector{Vector{Float64}}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1848
 [3] createmultipolygon(coords::Vector{Vector{Vector{Vector{Float64}}}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1848
 [4] convert(::Type{ArchGDAL.IGeometry}, type::GeoInterface.MultiPolygonTrait, geom::Shapefile.PolygonZ)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/geointerface.jl:352
 [5] convert(package::Module, geom::Shapefile.PolygonZ)
   @ GeoInterface ~/.julia/packages/GeoInterface/tuF7K/src/fallbacks.jl:149
 [6] top-level scope
   @ REPL[129]:1

Can we have a way of reproducing b::Shapefile.PolygonZ here?

felixcremer commented 4 months ago

Thanks for looking into it. I attached a reproducible example with a polygon that fails for me here.

julia> shpz = Shapefile.Handle("/home/fcremer/Daten/EDL_tests/shpz_test.shp")
Shapefile.Handle{Union{Missing, Shapefile.PolygonZ}}(Shapefile.Header(9994, 238, 1000, 15, Shapefile.Rect(668839.125, 5.177321e6, 668913.125, 5.177382e6), Shapefile.Interval(0.0, 0.0), Shapefile.Interval(-1.7976931348623157e308, -1.7976931348623157e308)), Union{Missing, Shapefile.PolygonZ}[Shapefile.PolygonZ(Shapefile.Rect(668839.125, 5.177321e6, 668913.125, 5.177382e6), Int32[0], Shapefile.Point[Shapefile.Point(668891.625, 5.177382e6), Shapefile.Point(668913.125, 5.177334e6), Shapefile.Point(668900.0, 5.1773235e6), Shapefile.Point(668879.0, 5.177321e6), Shapefile.Point(668841.75, 5.177328e6), Shapefile.Point(668839.125, 5.177347e6), Shapefile.Point(668846.3125, 5.177371e6), Shapefile.Point(668862.125, 5.177382e6), Shapefile.Point(668891.625, 5.177382e6)], Shapefile.Interval(0.0, 0.0), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], Shapefile.Interval(-1.7976931348623157e308, -1.7976931348623157e308), [-1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308])], ESRIWellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"ETRS_1989_UTM_Zone_32N\",GEOGCS[\"GCS_ETRS_1989\",DATUM[\"D_ETRS_1989\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",500000.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",9.0],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0]]"))

julia> GI.convert(AG, first(shpz.shapes))
ERROR: MethodError: no method matching addpoint!(::ArchGDAL.Geometry{ArchGDAL.wkbLineString}, ::Float64, ::Float64, ::Float64, ::Float64)

Closest candidates are:
  addpoint!(::G, ::Real, ::Real, ::Real) where G<:ArchGDAL.AbstractGeometry
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1298
  addpoint!(::G, ::Real, ::Real) where G<:ArchGDAL.AbstractGeometry
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1308

Stacktrace:
 [1] unsafe_createlinearring(coords::Vector{Vector{Float64}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1815
 [2] unsafe_createpolygon(coords::Vector{Vector{Vector{Float64}}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1848
 [3] createmultipolygon(coords::Vector{Vector{Vector{Vector{Float64}}}})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/ogr/geometry.jl:1848
 [4] convert(::Type{ArchGDAL.IGeometry}, type::GeoInterface.MultiPolygonTrait, geom::Shapefile.PolygonZ)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/geointerface.jl:352
 [5] convert(package::Module, geom::Shapefile.PolygonZ)
   @ GeoInterface ~/.julia/packages/GeoInterface/tuF7K/src/fallbacks.jl:149
 [6] top-level scope
   @ REPL[65]:1

shpz_test.zip

yeesian commented 4 months ago

@rafaqz following up on https://github.com/yeesian/ArchGDAL.jl/issues/416#issuecomment-1956066092, can I check if it's working-as-intended for GeoInterface.coordinates to be returning points with 4 dimensions for Shapefile.PolygonZ?

julia> GI.coordinates(first(shpz.shapes))
1-element Vector{Vector{Vector{Vector{Float64}}}}:
 [[[668891.625, 5.177382e6, 0.0, -1.7976931348623157e308], [668913.125, 5.177334e6, 0.0, -1.7976931348623157e308], [668900.0, 5.1773235e6, 0.0, -1.7976931348623157e308], [668879.0, 5.177321e6, 0.0, -1.7976931348623157e308], [668841.75, 5.177328e6, 0.0, -1.7976931348623157e308], [668839.125, 5.177347e6, 0.0, -1.7976931348623157e308], [668846.3125, 5.177371e6, 0.0, -1.7976931348623157e308], [668862.125, 5.177382e6, 0.0, -1.7976931348623157e308], [668891.625, 5.177382e6, 0.0, -1.7976931348623157e308]]]

Update 1: On reading https://www.esri.com/content/dam/esrisites/sitecore-archive/Files/Pdfs/library/whitepapers/pdfs/shapefile.pdf it seems like the spec for PolygonZ does include the buffer for M values, and that GeoInterface.coordinates(...) considers the M values as part of the coordinates.

julia> GI.is3d(shpz)
true

julia> GI.ismeasured(shpz)
true

Update 2: Seems like the answer is yes, and that this is hitting the TODOs in https://github.com/yeesian/ArchGDAL.jl/blob/dc19553f4d9e9eb03dca60a6467d2fae0a659f2a/src/ogr/geometry.jl#L1792-L1805

rafaqz commented 4 months ago

I suspected something like that.

But either way we shouldnt use GeoInterface.coordinates like this anymore.

It allocates hugely then its type unstable to splat the points into a function.

Better to use GI.x, GI.y, GI.z. I guess we throw a warning at the start saying the M is disgarded.

@felixcremer for your immmediate problem GeometryOps.jl has reproject too that uses Proj.jl directly and should work.

You might be able to convert it to GeometryBasics.jl or something first to drop the M values

felixcremer commented 4 months ago

Thanks for digging into it and thanks for the suggestion of GeometryOps. I will try that then.