GeoInterface.jl icon indicating copy to clipboard operation
GeoInterface.jl copied to clipboard

Makie plotting empty LibGEOS polygon fails

Open jaakkor2 opened this issue 2 years ago • 7 comments

using GLMakie, GeoInterfaceMakie, LibGEOS
GeoInterfaceMakie.@enable(LibGEOS.AbstractGeometry)
p1 = readgeom("POLYGON((0 0, 1 0, 0 1, 0 0))")
p2 = readgeom("POLYGON EMPTY")
plot(p1) # ok
plot(p2) # fails

gives

ERROR: MethodError: no method matching GeometryBasics.LineString(::Vector{Any})
Closest candidates are:
  GeometryBasics.LineString(::AbstractVector{<:Pair{P, P}}) where {N, T, P<:GeometryBasics.AbstractPoint{N, T}} at C:\Users\jaakkor2\.julia\packages\GeometryBasics\KE3OI\src\basic_types.jl:212
  GeometryBasics.LineString(::AbstractVector{<:GeometryBasics.AbstractPoint}) at C:\Users\jaakkor2\.julia\packages\GeometryBasics\KE3OI\src\basic_types.jl:208
  GeometryBasics.LineString(::AbstractVector{<:GeometryBasics.AbstractPoint}, ::AbstractVector{<:Integer}) at C:\Users\jaakkor2\.julia\packages\GeometryBasics\KE3OI\src\basic_types.jl:241
  ...
Stacktrace:
  [1] convert(#unused#::Type{GeometryBasics.LineString}, type::GeoInterface.LineStringTrait, geom::LinearRing)
    @ GeometryBasics C:\Users\jaakkor2\.julia\packages\GeometryBasics\KE3OI\src\geointerface.jl:90
  [2] convert(#unused#::Type{GeometryBasics.Polygon}, type::GeoInterface.PolygonTrait, geom::Polygon)
    @ GeometryBasics C:\Users\jaakkor2\.julia\packages\GeometryBasics\KE3OI\src\geointerface.jl:95
  [3] basicsgeom
    @ C:\Users\jaakkor2\.julia\packages\GeoInterfaceMakie\1cwTc\src\GeoInterfaceMakie.jl:55 [inlined]
  [4] _convert_arguments
    @ C:\Users\jaakkor2\.julia\packages\GeoInterfaceMakie\1cwTc\src\GeoInterfaceMakie.jl:66 [inlined]
  [5] convert_arguments
    @ C:\Users\jaakkor2\.julia\packages\GeoInterfaceMakie\1cwTc\src\GeoInterfaceMakie.jl:80 [inlined]
  [6] plot!(scene::Scene, P::Type{Any}, attributes::Attributes, args::Polygon; kw_attributes::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Makie C:\Users\jaakkor2\.julia\packages\Makie\LCBYx\src\interfaces.jl:317
  [7] plot!(scene::Scene, P::Type{Any}, attributes::Attributes, args::Polygon)
    @ Makie C:\Users\jaakkor2\.julia\packages\Makie\LCBYx\src\interfaces.jl:302
  [8] plot(P::Type{Any}, args::Polygon; axis::NamedTuple{(), Tuple{}}, figure::NamedTuple{(), Tuple{}}, kw_attributes::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Makie C:\Users\jaakkor2\.julia\packages\Makie\LCBYx\src\figureplotting.jl:48
  [9] plot
    @ C:\Users\jaakkor2\.julia\packages\Makie\LCBYx\src\figureplotting.jl:31 [inlined]
 [10] #plot#13
    @ C:\Users\jaakkor2\.julia\packages\MakieCore\77C6Z\src\recipes.jl:34 [inlined]
 [11] plot(args::Polygon)
    @ MakieCore C:\Users\jaakkor2\.julia\packages\MakieCore\77C6Z\src\recipes.jl:33
 [12] top-level scope
    @ REPL[19]:1

with

⌃ [e9467ef8] GLMakie v0.7.4
  [0edc0954] GeoInterfaceMakie v0.1.0
  [a90b1aa1] LibGEOS v0.7.4

jaakkor2 avatar Dec 07 '22 15:12 jaakkor2

So this is a problem with the GeometryBasics.jl convert method.

https://github.com/JuliaGeometry/GeometryBasics.jl/blob/db65c41fd75f478b4079d8f4825784f307af024d/src/geointerface.jl#L90

Needs to be changed to: LineString(Point{dim, Float64}[Point{dim, Float64}(GeoInterface.coordinates(p)) for p in getgeom(geom)])

So the vector type is set even when its empty. In fact most of the vectors in the file need that treatement, e.g.

https://github.com/JuliaGeometry/GeometryBasics.jl/blob/db65c41fd75f478b4079d8f4825784f307af024d/src/geointerface.jl#L111

rafaqz avatar Dec 07 '22 16:12 rafaqz

I looked into this: we need to think about what to do with empty geometries in getgeom.

It's currently not clear what to return as we don't know the eltype.

We may need to implement at better iterator that has a know eltype even when its empty.

rafaqz avatar Jan 01 '23 22:01 rafaqz

Maybe we could explicitly create an empty geometry for polygon https://github.com/jaakkor2/GeometryBasics.jl/blob/emptypolygon/src/geointerface.jl#L94 and for multipolygon https://github.com/jaakkor2/GeometryBasics.jl/blob/emptypolygon/src/geointerface.jl#L116

But now I am a bit puzzled what the dimension should be since GeometryBasics.ncoord returns zero for empty polygons.

jaakkor2 avatar Jan 02 '23 06:01 jaakkor2

Yeah, there are a lot of issues with zero length polygons and other geometries.

Base julia uses some dark arts like return value inference to make zero length iterators make sense. We probably don't want to do that. But I think this might need some serious design work to solve this in a general way.

But: does plotting zero length polygons even work in Makie.jl with GeometryBasics.Polygon??

rafaqz avatar Jan 02 '23 21:01 rafaqz

using GLMakie
using GLMakie.Makie.GeometryBasics
p1 = Polygon(Point2{Float64}[])
p2 = Polygon([Point2(NaN, NaN)])

plot(p1) errors, but plot(p2) works.

jaakkor2 avatar Jan 03 '23 08:01 jaakkor2

Ok so what you want isn't possible in Makie unless we do that NaN point hack.

@SimonDanisch any ideas here? should Makie be able to plot zero length polygons?

rafaqz avatar Jan 17 '23 12:01 rafaqz