Makie plotting empty LibGEOS polygon fails

[jaakkor2]

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

[rafaqz]

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]

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.

[jaakkor2]

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.

[rafaqz]

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??

[jaakkor2]

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

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

[rafaqz]

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?

[asinghvi17]

I've implemented the NaN-point hack for missing geometries in GeoInterfaceMakie (see #139). I guess we could also check whether the geometry is empty or not and then convert appropriately? That's not going in #139 but it could be a separate PR.