Convert Shapefile.PolygonZ to ArchGDAL geometry

[felixcremer]

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]

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]

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!:

ArchGDAL.jl/src/ogr/geometry.jl

Lines 1263 to 1272 in dc19553

	function setpoint!(
	geom::G,
	i::Integer,
	x::Real,
	y::Real,
	z::Real,
	)::G where {G<:AbstractGeometry}
	GDAL.ogr_g_setpoint(geom, i, x, y, z)
	return geom
	end

[yeesian]

Someone is calling GI.coordinates.

Is it this line?

ArchGDAL.jl/src/geointerface.jl

Line 352 in dc19553

return f(GeoInterface.coordinates(geom))

[yeesian]

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]

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]

@rafaqz following up on #416 (comment), 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

ArchGDAL.jl/src/ogr/geometry.jl

Lines 1792 to 1805 in dc19553

	# TODO make measures work
	# @eval function $f1(coords::Tuple{<:Real,<:Real}, m)
	# geom = $f1(Val{wkbPointM}())
	# addpoint!(geom, coords...)
	# return geom
	# end
	# @eval function $f1(coords::Tuple{<:Real,<:Real,<:Real}, m)
	# geom = $f1(Val{wkbPointZM}())
	# addpoint!(geom, coords...)
	# return geom
	# end
	# Coordinates can be Vector of Real or
	# TODO handle M and 25D

[rafaqz]

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]

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