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(
    Float16,
    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...)),
    mean,
    (5, 5),
)
layer3.left = round(layer3.left; digits=0)
layer3.right = round(layer3.right; digits=0)

30.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])

Using the subplot and inset arguments of Plots.jl, we can have the legend within the figure. Note how in this example we expand the limits on the x axis to make the legend fit, but also use more classes in the map to have a smoother result.

p1 = bivariate(
    layer1,
    layer3;
    classes=6,
    bv_pal_2...,
    frame=:box,
    xlim=(-24, maximum(longitudes(layer1))),
)
xaxis!(p1, "Longitude")
yaxis!(p1, "Latitude")
p2 = bivariatelegend!(
    layer1,
    layer3;
    bv_pal_2...,
    inset=(1, bbox(0.04, 0.05, 0.28, 0.28, :top, :left)),
    subplot=2,
    xlab=layernames(EarthEnv, HabitatHeterogeneity, 2),
    ylab=layernames(EarthEnv, LandCover, 9),
    guidefontsize=7,
)

Using a trivariate mapping follows the same process, with layers representing the red, green, and blue channel respectively.

plot(layer1, layer2, layer3; st=:trivariate, frame=:grid)

There are two options for this type of plots. The first is quantiles=true (which maps quantiles rather than raw values), and the second is simplex=false, which makes all values sum to 1 within a pixel. For example:

trivariate(layer1, layer2, layer3; quantiles=true, simplex=true, frame=:grid)

It is a good idea to question whether using simplex is appropriate. The legend can also be plotted using trivariatelegend:

trivariatelegend(layer1, layer2, layer3; quantiles=true, simplex=true)

The legend function admits three additional arguments for the names of the red, green, and blue channels:

trivariatelegend(
    layer1,
    layer2,
    layer3;
    quantiles=true,
    simplex=true,
    red="Heterogeneous",
    green="Rough",
    blue="Urbanized",
)

We can also combine the two elements:

tri1 = trivariate(
    layer1,
    layer2,
    layer3;
    xlim=(-24, maximum(longitudes(layer1))),
    frame=:grid,
    grid=false,
    simplex=false,
)
xaxis!(tri1, "Longitude")
yaxis!(tri1, "Latitude")
tri2 = trivariatelegend!(
    layer1,
    layer2,
    layer3;
    inset=(1, bbox(0.01, 0.03, 0.35, 0.35, :top, :left)),
    subplot=2,
    red="Heterogeneity",
    green="Roughness",
    blue="Urban",
    annotationfontsize=6,
    simplex=false,
)