Plotting maps using GeoTables.jl, Meshes.jl and MeshViz.jl

[greimel]

This is just to report that there might be an alternative way to support plotting maps.

Advantage:

• GeoTables.jl can read map data from .shp files but can also load map data for any country from the GADM website.

Issues:

• Unfortunately, we cannot handle the case that the Meshes.Viz recipe uses elementcolor instead of color (it does so because it also supports boundarycolor, facetcolor, vertexcolor). So AoG doesn't show a legend.
• The recipe doesn't reduce the fineness of the polygon area, so one needs to set the decimation keyword manually. Otherwise plotting maps can take very long (I haven't come across this when using Makie.poly). (A single Brazilian state took more than a minute for me without using decimation)

import GeoTables
using MeshViz: Viz
using DataFrames: DataFrame
using CairoMakie, AlgebraOfGraphics

import Meshes
const AoG = AlgebraOfGraphics
AoG.scientific_type(::Type{<:Meshes.Geometry}) = AoG.geometrical


t = GeoTables.gadm("AUT", children=true)
df = DataFrame(t)
df.color = 1:size(df, 1)

data(df) * visual(Viz, decimation = 0.01) * mapping(:geometry, elementcolor = :color) |> draw

cc @juliohm

[juliohm]

Thank you for pinging @greimel , you are correct that we need decimation with the GADM dataset because it is extremely refined. Most other polygons in the wild are simplified versions already. What we plan to do in MeshViz.jl is set a more sensible default for decimation instead of no decimation. I may actually fix this in the following weeks because of a workshop we are doing at our institute.

The default I have in mind is a binary search until the number of vertices is smaller than a pre-specified threshhold. I am open to suggestions.

[greimel]

Automatic decimation would be great!

Could you make it dependent on the size of the bounding box? Say, skip points unless they are more than box_size .* 0.01 away from the previous one?

Ideally you would use the size of the union of bounding boxes. (Is that possible?)

[juliohm]

Ideally we would set the decimation as a function of the visualization zoom level as well so that the geometries within the visualization window at a given zoom level get the necessary level of detail. That requires more advanced Makie knowledge, which I don't have unfortunately.

If you have ideas on the default decimation method, please feel free to submit a PR, the code should be super simple to read. We need a function default_decimation(geometry) that does something clever and returns a reasonable threshold. Maybe in the future it becomes default_decimation(geometry, scene). This threshold is passed to the simplify function from Meshes.jl before plotting.