Plotting polygons defined in degrees

Specific Domains Visualization
1

I’m trying to plot the polygons from GADM.jl. I do the following:

using GeoTables
using GeoArtifacts
using Countries
using DataFrames
using Meshes
using GLMakie 
using GeoMakie

> gdf = GADM.get("FIN"; level=0)
> p = reduce(vcat, gdf.geometry)
> typeof(p.geoms[1])
 PolyArea{🌐, CoordRefSystems.GeodeticLatLon{CoordRefSystems.WGS84Latest, Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}, Ring{🌐, CoordRefSystems.GeodeticLatLon{CoordRefSystems.WGS84Latest, Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}, CircularArrays.CircularVector{Meshes.Point{🌐, CoordRefSystems.GeodeticLatLon{CoordRefSystems.WGS84Latest, Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}}, Vector{Meshes.Point{🌐, CoordRefSystems.GeodeticLatLon{CoordRefSystems.WGS84Latest, Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}}}}}, Vector{Ring{🌐, CoordRefSystems.GeodeticLatLon{CoordRefSystems.WGS84Latest, Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}, CircularArrays.CircularVector{Meshes.Point{🌐, CoordRefSystems.GeodeticLatLon{CoordRefSystems.WGS84Latest, Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}}, Vector{Meshes.Point{🌐, CoordRefSystems.GeodeticLatLon{CoordRefSystems.WGS84Latest, Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}}}}}}}

This is not plottable by viz!() (from Meshes.jl I believe). Simply nothing is plotted. I can plot PolyAreas when the coordinates are plain numbers, not (:lat, :lon) tuples in degree units.

2

Check Optimize visualization of non-Euclidean geometries and domains · Issue #1367 · JuliaGeometry/Meshes.jl · GitHub for additional information.

The current workaround is to project the spherical geometries into “flat” Euclidean geometries before plotting, and that is already optimized:

julia> using GeoStats

julia> using GeoArtifacts

julia> import GLMakie as Mke

julia> gtb = GADM.get("FIN")
               1×3 GeoTable over 1 GeometrySet
┌─────────────┬─────────────┬───────────────────────────────┐
│    GID_0    │   COUNTRY   │           geometry            │
│ Categorical │ Categorical │         MultiPolygon          │
│  [NoUnits]  │  [NoUnits]  │ 🖈 GeodeticLatLon{WGS84Latest} │
├─────────────┼─────────────┼───────────────────────────────┤
│     FIN     │   Finland   │     Multi(2879×PolyArea)      │
└─────────────┴─────────────┴───────────────────────────────┘

julia> ptb = gtb |> Proj(Robinson)
            1×3 GeoTable over 1 GeometrySet
┌─────────────┬─────────────┬─────────────────────────┐
│    GID_0    │   COUNTRY   │        geometry         │
│ Categorical │ Categorical │      MultiPolygon       │
│  [NoUnits]  │  [NoUnits]  │ 🖈 Robinson{WGS84Latest} │
├─────────────┼─────────────┼─────────────────────────┤
│     FIN     │   Finland   │  Multi(2879×PolyArea)   │
└─────────────┴─────────────┴─────────────────────────┘

julia> viewer(ptb)

image
image730×608 35.8 KB

If you have the time and skill to contribute, please reach out in our community channels.

3

Ok this works. Although the color is not what I expected.

image
image331×497 15.2 KB

4

Please share the script used, otherwise it is impossible to debug. The color is assigned by the viewer or you can specify it with the lower level viz function.

I highly recommend reading the GDSJL book:

5

Another simple way to plot this is

using DataFrames
using GADM
using GeoMakie
using CairoMakie  # Or GLMakie should work
gdf = GADM.get("FIN"; depth=0) |> DataFrame

fig = Figure()
ga = GeoAxis(fig[1, 1]; dest="+proj=eqdc +datum=WGS84 +lat_1=60 +lat_2=71")
poly!(ga, gdf.geom; color=:blue, strokecolor=:darkblue, strokewidth=1, shading=NoShading)
fig
6

Here is the code:

using GeoTables
using GeoArtifacts
using Countries
using DataFrames
using Meshes
using GLMakie
using GeoMakie
using GeoStats
p = GADM.get(“FIN”, level = 0)
viz(p.geometry |> Meshes.Proj(Mercator), color = :blue)

7

I cannot reproduce the issue:

image
image298×436 12.9 KB

Notice that GADM.get doesn’t have a level option. You’re probably confusing the GADM.jl package that we used to maintain in the past with the new GeoArtifacts.GADM submodule.

Please remove all imports from other stacks to avoid confusion. All you need is GeoStats.jl, GeoArtifacts.jl and a Makie.jl backend.

8

And another even simpler is

using GMT

data = gadm("IND");
viz(data, proj=:guess)

GMTjl_j
GMTjl_j1912×3411 292 KB

9

gadm() is different from GADM.get()? Why we need not the Proj() function in this case?

10

Ok, it seems to be depth instead. The color problem seems to be specific to GLMakie, not CairoMakie. And colored fills, not grey ones.

11

Are you using llvmpipe? That’s known to cause color issues

12

Notice that I used GLMakie.jl in my machine, so the problem is not there.

13

It could be machine-specific as @jules pointed out.

14

This is a different function from the GMT.jl package, which is a wrapper to an external library (not written in Julia). Unfortunately they decided to call their command viz after we picked that name, so users get even more confused.

Notice that they call proj=:guess as well. GMT doesn’t support interactive 3D scenes nor spherical geometries as we do.

15

Yes I am.