GeoInterface.jl
GeoInterface.jl copied to clipboard
Makie plotting empty LibGEOS polygon fails
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
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
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.
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.
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
??
using GLMakie
using GLMakie.Makie.GeometryBasics
p1 = Polygon(Point2{Float64}[])
p2 = Polygon([Point2(NaN, NaN)])
plot(p1)
errors, but plot(p2)
works.
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?