broken GeoInterface.convert(GeometryBasics.MultiPolygon)

[KingBoomie]

Hey, sorry if this is in the wrong repo, but I'm unsure which package is actually the culprit here. possibilities include GeometryBasics, GeoInterface and ArchGDAL

when trying to plot a country polygon from GADM with makie:

using GADM, GeoInterface
using GeometryBasics: MultiPolygon

gdalgeom = only(GADM.get("CHN").geom)
basicgeom = GeoInterface.convert(MultiPolygon, gdalgeom)
# actual plotting is unnecessary for reproduction, but motivates this example
# poly(basicgeom, background_color=:transparent)

this errors out with

click here to see error

ERROR: MethodError: GeometryBasics.Point2{Float64}(::Base.Generator{UnitRange{Int64}, GeoInterface.var"#21#22"{PointTrait, ArchGDAL.IGeometry{ArchGDAL.wkbPoint}}}) is ambiguous. Candidates:
  (::Type{SA})(gen::Base.Generator) where SA<:StaticArraysCore.StaticArray in StaticArrays at /home/kris/.julia/packages/StaticArrays/PUoe1/src/SArray.jl:54
  GeometryBasics.Point{S, T}(x) where {S, T} in GeometryBasics at /home/kris/.julia/packages/GeometryBasics/6JxlJ/src/fixed_arrays.jl:44
Possible fix, define
  GeometryBasics.Point{S, T}(::Base.Generator) where {S, T}
Stacktrace:
  [1] (::GeometryBasics.var"#218#219"{Int64})(p::ArchGDAL.IGeometry{ArchGDAL.wkbPoint})
    @ GeometryBasics ./none:0
  [2] iterate
    @ ./generator.jl:47 [inlined]
  [3] collect(itr::Base.Generator{Base.Generator{UnitRange{Int64}, GeoInterface.var"#19#20"{LineStringTrait, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}}}, GeometryBasics.var"#218#219"{Int64}})
    @ Base ./array.jl:787
  [4] convert(#unused#::Type{GeometryBasics.LineString}, type::LineStringTrait, geom::ArchGDAL.IGeometry{ArchGDAL.wkbLineString})
    @ GeometryBasics ~/.julia/packages/GeometryBasics/6JxlJ/src/geointerface.jl:113
  [5] convert(#unused#::Type{GeometryBasics.Polygon}, type::PolygonTrait, geom::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon})
    @ GeometryBasics ~/.julia/packages/GeometryBasics/6JxlJ/src/geointerface.jl:119
  [6] (::GeometryBasics.var"#224#225"{PolygonTrait})(poly::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon})
    @ GeometryBasics ./none:0
  [7] iterate
    @ ./generator.jl:47 [inlined]
  [8] collect(itr::Base.Generator{Base.Generator{UnitRange{Int64}, GeoInterface.var"#19#20"{MultiPolygonTrait, ArchGDAL.IGeometry{ArchGDAL.wkbMultiPolygon}}}, GeometryBasics.var"#224#225"{PolygonTrait}})
    @ Base ./array.jl:787
  [9] convert(#unused#::Type{MultiPolygon}, type::MultiPolygonTrait, geom::ArchGDAL.IGeometry{ArchGDAL.wkbMultiPolygon})
    @ GeometryBasics ~/.julia/packages/GeometryBasics/6JxlJ/src/geointerface.jl:140
 [10] convert(T::Type, geom::ArchGDAL.IGeometry{ArchGDAL.wkbMultiPolygon})
    @ GeoInterface ~/.julia/packages/GeoInterface/J298z/src/interface.jl:612
 [11] top-level scope
    @ REPL[4]:1

I'm currently reporting it here, since I fixed it by modifying GeoInterface.convert(::Type{LineString}, type::LineStringTrait, geom) defined here

- function GeoInterface.convert(::Type{LineString}, type::LineStringTrait, geom)
-    dim = Int(ncoord(geom))
-    return LineString([Point{dim, Float64}(GeoInterface.coordinates(p)) for p in getgeom(geom)])
- end
+ function GeoInterface.convert(::Type{LineString}, type::LineStringTrait, geom)
+     dim = Int(ncoord(geom))
+     if dim == 2
+         N = GeoInterface.npoint(geom)
+         points = Array{Point2f}(undef, N)
+         for (i,p) in enumerate(getgeom(geom))
+             x,y = GeoInterface.coordinates(p)
+             points[i] = Point2f(x,y)
+         end
+         return LineString(points)
+ 
+     elseif dim == 3
+         N = GeoInterface.npoint(geom)
+         points = Array{Point3f}(undef, N)
+         for (i,p) in enumerate(getgeom(geom))
+             x,y,z = GeoInterface.coordinates(p)
+             points[i] = Point3f(x,y,z)
+         end
+         return LineString(points)
+ 
+     else
+         error("geom dim needs to be 2 or 3")
+     end
+ end

This is not a PR because I'm pretty sure this isn't the correct way to solve this problem (it just works for the combination of packages I have) and I'm unsure on how to go about testing a whole Ecosystem of Packages. Hoping a maintainer can pick it up and fix it, but if not, it's atleast a documented workaround, should someone stumble upon a similar problem.

[sjkelly]

The error message suggests we define:

  GeometryBasics.Point{S, T}(::Base.Generator) where {S, T}

Because:

 (::Type{SA})(gen::Base.Generator) where SA<:StaticArraysCore.StaticArray in StaticArrays at /home/kris/.julia/packages/StaticArrays/PUoe1/src/SArray.jl:54
  GeometryBasics.Point{S, T}(x) where {S, T} in GeometryBasics at /home/kris/.julia/packages/GeometryBasics/6JxlJ/src/fixed_arrays.jl:44

Are ambiguous.

[rafaqz]

I think there is a simpler fix than either of these, we just need to be more a little more explicit than using coordinates.

@KingBoomie if you have time to make a PR, we can check dimensionality of the geometry using GeoInterface.is3d

Something like this is a little simpler, more often type stable, and should work similarly to your code:

function GeoInterface.convert(::Type{LineString}, type::LineStringTrait, geom)
    points = if GeoInterface.is3d(geom) 
       [Point{3, Float64}(GeoInterface.x(p), GeoInterface.y(p), GeoInterface.z(p)) for p in getgeom(geom)]
    else
       [Point{2, Float64}(GeoInterface.x(p), GeoInterface.y(p)) for p in getgeom(geom)]
    end
    return LineString(points)
end

Probably defining the Array with undef as you are is actually faster than this though, so feel free to use that instead.