Multivariate maps

using SimpleSDMLayers
using StatsPlots
using StatsPlots.PlotMeasures
using Statistics

Justification for this use case: We can show more than one (specifically, two or three) variables on a single map, using a bivariate or trivariate color scale. In order to illustrate these mappings, we will look at the joint distribution of three measures: eveness of land use, terrain roughness, and proportion of urbanized land.

boundaries = (left=-12.0, right=30.0, bottom=36.0, top=72.0)
layer1 = convert(
    SimpleSDMPredictor(EarthEnv, HabitatHeterogeneity, 2; resolution=5, boundaries...),
layer2 = convert(
    Float16, SimpleSDMPredictor(EarthEnv, Topography, 7; resolution=5, boundaries...)
SDM response → 864×1008 grid with 870912 Float16-valued cells
  Latitudes	36.0 ⇢ 72.0
  Longitudes	-12.0 ⇢ 30.0

The landcover data are finer than the other layers, so we will coarsen them. Because of rounding issues, the left and right cordinates need to be rounded.

layer3 = coarsen(
    convert(Float16, SimpleSDMPredictor(EarthEnv, LandCover, 9; boundaries...)),
    (5, 5),
layer3.left = round(layer3.left; digits=0)
layer3.right = round(layer3.right; digits=0)

We finally mask everything according to the first layer

layer2 = mask(layer1, layer2);
layer3 = mask(layer1, layer3);

Note that bivariate maps usually work best when used with 9 classes in total (so 3 for each side). The next decision is to take a bivaraite color palette, and the combinations below are commonly used. Note that you can definitely use diverging colors if you want. If you use colors in the RGBA format (e.g. colorant"#ef0ce8c4"), the color map will account for transparency.

p0 = colorant"#e8e8e8"
bv_pal_1 = (p0=p0, p1=colorant"#64acbe", p2=colorant"#c85a5a")
bv_pal_2 = (p0=p0, p1=colorant"#73ae80", p2=colorant"#6c83b5")
bv_pal_3 = (p0=p0, p1=colorant"#9972af", p2=colorant"#c8b35a")
bv_pal_4 = (p0=p0, p1=colorant"#be64ac", p2=colorant"#5ac8c8")
(p0 = RGB{N0f8}(0.91,0.91,0.91), p1 = RGB{N0f8}(0.745,0.392,0.675), p2 = RGB{N0f8}(0.353,0.784,0.784))

The bivariate map itself is a call to plot. Internally, this will transform the layers into quantiles (determined by the classes keyword, defaults to 3):

plot(layer1, layer3; st=:bivariate, bv_pal_3...)

Note that you can use the bivariate shorthand as well:

pl1 = bivariate(layer1, layer3; classes=3, frame=:box, bv_pal_4...)
xaxis!(pl1, "Longitude")
yaxis!(pl1, "Latitude")

We can repeat essentially the same process for the legend:

pl2 = bivariatelegend(layer1, layer3; classes=3, bv_pal_4...)
xaxis!(pl2, layernames(EarthEnv, HabitatHeterogeneity, 2))
yaxis!(pl2, layernames(EarthEnv, LandCover, 9))

And now, we can plot the legend next to the map - future releases of the package will hopefully offer this in a far more user friendly way.

plot(pl1, pl2; layout=@layout [a{0.75w} b])