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

Makie plotting empty LibGEOS polygon fails #81

Open
jaakkor2 opened this issue Dec 7, 2022 · 7 comments
Open

Makie plotting empty LibGEOS polygon fails #81

jaakkor2 opened this issue Dec 7, 2022 · 7 comments

Comments

@jaakkor2
Copy link
Contributor

jaakkor2 commented Dec 7, 2022

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

rafaqz commented Dec 7, 2022

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

rafaqz commented Jan 1, 2023

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

jaakkor2 commented Jan 2, 2023

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

rafaqz commented Jan 2, 2023

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

jaakkor2 commented Jan 3, 2023

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

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

@rafaqz
Copy link
Member

rafaqz commented Jan 17, 2023

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

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.

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