Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

broken GeoInterface.convert(GeometryBasics.MultiPolygon) #182

Open
KingBoomie opened this issue Oct 19, 2022 · 2 comments
Open

broken GeoInterface.convert(GeometryBasics.MultiPolygon) #182

KingBoomie opened this issue Oct 19, 2022 · 2 comments

Comments

@KingBoomie
Copy link

KingBoomie commented Oct 19, 2022

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
Copy link
Member

sjkelly commented Oct 19, 2022

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
Copy link
Contributor

rafaqz commented Oct 24, 2022

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants