JuliaGeometry / JuliaGeometry/GeometryBasics.jl

broken GeoInterface.convert(GeometryBasics.MultiPolygon)

Open
#182 2 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
Julia
Stars
215
Forks
64
Avg merge
3d 18h
Merged PRs (30d)
2

Description

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.

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start in src/geointerface.jl at the GeoInterface.convert(::Type{LineString}, ...) implementation referenced by the report, then reproduce the MultiPolygon conversion with GADM, GeoInterface, GeometryBasics, and ArchGDAL. Compare the reported workaround with the existing conversion path and add or run a focused regression test for the MultiPolygon example. Done means the conversion no longer raises the ambiguous Point2 constructor error.

Written by the indexing model from the issue text.

Assessment

Tech stack
julia
Domain
computer-graphics
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Stale
Clarity
Mostly clear
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.