Generating virtual species

In this vignette, we provide a demonstration of how the different SpeciesDistributionToolkit functions can be chained together to rapidly create a virtual species, generate its range map, and sample points from it according to the predicted suitability.

using SpeciesDistributionToolkit
using CairoMakie
using Statistics

We start by defining the extent in which we want to create the virtual species. For the purpose of this example, we will use the country of Paraguay. Note that the boundingbox function returns the coordinates in WGS84.

place = SpeciesDistributionToolkit.openstreetmap("Paraguay")
extent = SpeciesDistributionToolkit.boundingbox(place)

(left = -62.64066696166992, right = -54.257999420166016, bottom = -27.606393814086914, top = -19.287647247314453)

We then download some environmental data. In this example, we use the BioClim variables as distributed by CHELSA. In order to simplify the code, we will only use BIO1 (mean annual temperature) and BIO12 (total annual precipitation). Note that we collect these layers in a vector typed as SDMLayer{Float32}, in order to ensure that future operations already recevie floating point values.

provider = RasterData(CHELSA2, BioClim)
L = SDMLayer{Float32}[SDMLayer(provider; layer = l, extent...) for l in ["BIO1", "BIO12"]]

2-element Vector{SDMLayer{Float32}}:
 🗺️  A 999 × 1007 layer (1005993 Float32 cells)
 🗺️  A 999 × 1007 layer (1005993 Float32 cells)

We now mask the layers using the polygons we downloaded initially. Here, this is done in two steps, first the masking of the first layer, and second the masking of all other layers. Currently unreleased versions of the package have a shortcut for this operation.

rescale!.(mask!(L, place))

2-element Vector{SDMLayer{Float32}}:
 🗺️  A 999 × 1007 layer (508477 Float32 cells)
 🗺️  A 999 × 1007 layer (508477 Float32 cells)

In the next steps, we will generate some virtual species. These are defined by an environmental response to each layer, linking the value of the layer at a point to the suitability score. For the sake of expediency, we only use logistic responses, and generate one function for each layer (drawing α from a normal distribution, and β uniformly).

logistic(x, α, β) = 1 / (1 + exp((x - β) / α))
logistic(α, β) = (x) -> logistic(x, α, β)

logistic (generic function with 2 methods)

We will next write a function that will sample a value in each layer to serve as the center of the logistic response, draw a shape parameter at random, and then return the suitability of the environment for a virtual species under the assumption of a set prevalence:

function virtualspecies(L::Vector{<:SDMLayer}; prevalence::Float64 = 0.25)
    prevalence = clamp(prevalence, eps(), 1.0)
    invalid = true
    while invalid
        centers = [rand(values(l)) for l in L]
        shapes = randn(length(L))
        predictors = [logistic(centers[i], shapes[i]) for i in eachindex(L)]
        predictions = [predictors[i].(L[i]) for i in eachindex(L)]
        rescale!.(predictions)
        global prediction = prod(predictions)
        rescale!(prediction)
        invalid = any(isnan, prediction)
    end
    cutoff = quantile(prediction, 1 - prevalence)
    return prediction, cutoff, prediction .>= cutoff
end

virtualspecies (generic function with 1 method)

We can now apply this new function to create a simulated species distribution:

vsp, τ, vrng = virtualspecies(L; prevalence = 0.33)

(🗺️  A 999 × 1007 layer (508477 Float64 cells), 0.24334662924170145, 🗺️  A 999 × 1007 layer (508477 Bool cells))

We can plot the output of this function to inspect what the generated range looks like:

Code for the figure

f = Figure(; size = (700, 300))
ax1 = Axis(f[1, 1]; aspect = DataAspect(), title = "Score")
ax2 = Axis(f[1, 2]; aspect = DataAspect(), title = "Range")
heatmap!(ax1, vsp; colormap = :navia)
heatmap!(ax2, vrng; colormap = :Greens)
for ax in [ax1, ax2]
    tightlimits!(ax)
    hidedecorations!(ax)
    hidespines!(ax)
    lines!(ax, place; color = :black)
end

Random observations for the virtual species are generated by setting the probability of inclusion to 0 for all values above the cutoff, and then sampling proportionally to the suitability for all remaining points. Note that the method is called backgroundpoints, as it is normally used for pseudo-absences. The second argument of this method is the number of points to generate.

presencelayer = backgroundpoints(mask(vsp, nodata(vrng, false)), 100)

🗺️  A 999 × 1007 layer with 167798 Bool cells
   Projection: +proj=longlat +datum=WGS84 +no_defs

We can finally plot the result:

Code for the figure

f = Figure(; size = (700, 700))
ax = Axis(f[1, 1]; aspect = DataAspect())
heatmap!(ax, vrng; colormap = :Greens)
lines!(ax, place; color = :black)
scatter!(
    ax,
    presencelayer;
    color = :white,
    strokecolor = :black,
    strokewidth = 2,
    markersize = 10,
    label = "Virtual presences",
)
tightlimits!(ax)
hidespines!(ax)
hidedecorations!(ax)
axislegend(ax; position = :lb, framevisible = false)

These data could, for example, be used to benchmark species distribution models. For the analysis presented in the manuscript, we are interested in applying the simulation of virtual species to a large number, in order to say something about potential patterns of biodiversity. For this reason, we will now simulate 100 species, with prevalences drawn uniformly in the unit interval.

ranges = [virtualspecies(L; prevalence = rand())[3] for _ in 1:100]
first(ranges)

🗺️  A 999 × 1007 layer with 508477 Bool cells
   Projection: +proj=longlat +datum=WGS84 +no_defs

We can sum these layers to obtain a measurement of the simulated species richness:

Code for the figure

f = Figure(; size = (700, 700))
richness = mosaic(sum, ranges)
ax = Axis(f[1, 1]; aspect = DataAspect())
hm = heatmap!(ax, richness; colorrange = (0, length(ranges)), colormap = :navia)
Colorbar(
    f[1, 1],
    hm;
    label = "Species richness",
    alignmode = Inside(),
    width = Relative(0.4),
    valign = :bottom,
    halign = :left,
    tellheight = false,
    tellwidth = false,
    vertical = false,
)
hidedecorations!(ax)
lines!(ax, place; color = :black)
hidespines!(ax)

We can now transform these data into a partition of the contribution of each species and location to the total beta-diversity:

function LCBD(ranges::Vector{SDMLayer{Bool}}; transformation::Function = identity)
    Y = transformation(hcat(values.(ranges)...))
    S = (Y .- mean(Y; dims = 1)) .^ 2.0
    SStotal = sum(S)
    BDtotal = SStotal / (size(Y, 1) - 1)
    SSj = sum(S; dims = 1)
    SCBDj = SSj ./ SStotal
    SSi = sum(S; dims = 2)
    LCBDi = SSi ./ SStotal
    betadiv = similar(first(ranges), Float32)
    betadiv.grid[findall(betadiv.indices)] .= vec(LCBDi)
    return betadiv, vec(SCBDj), BDtotal
end

LCBD (generic function with 1 method)

For good measure, we will also add a function to perform the Hellinger transformation:

function hellinger(Y::AbstractMatrix{T}) where {T <: Number}
    yi = sum(Y; dims = 2) .+ 1
    return sqrt.(Y ./ yi)
end

hellinger (generic function with 1 method)

We can now apply this function to our simulated ranges. The first output is a layer with the local contributions to beta diversity (LCBD), the second is a vector with the contribution of species, and the last element is the total beta diversity.

βl, βs, βt = LCBD(ranges; transformation = hellinger)

(🗺️  A 999 × 1007 layer (508477 Float32 cells), [0.01178323642784345, 0.026675539059340166, 0.012231283085205097, 0.0022990213034466246, 0.0011023945559805444, 0.01069322347663175, 0.010327236427132654, 0.0032294613578439974, 0.011505484686242538, 0.008245509826359257, 0.020776058427294086, 0.014555426002917286, 0.009553924244099652, 0.008138519087395324, 0.011443642902892254, 0.008923018269658163, 0.004385458806021812, 0.0006932748173793108, 0.014750643065308667, 0.008764337739414782, 0.010465708023939388, 0.0032356812296776482, 0.012849101865606219, 0.004059855325529815, 0.02481454882451712, 0.007262483308043095, 0.011477960989853563, 0.003673295909780829, 0.0048005535229008155, 0.009495366614010627, 0.013099877246023536, 0.010262852865087265, 0.008147710195024456, 0.005915855522721753, 0.011140155452291684, 0.011491863695109352, 0.006579985965477292, 0.026559894980905174, 0.0039019760124425094, 0.011057244444408292, 0.015424085483571383, 0.00845632958473842, 0.015970123311793608, 0.009423498011370896, 0.0012097273351671896, 0.01136634320314775, 0.013616759594375718, 0.008417744680479131, 0.010862994228055619, 0.0010281371339471725, 0.010727302546632293, 0.011127018731780153, 0.0022341862658990805, 0.009612446053336995, 0.011885582020373823, 0.007308538014713277, 0.010323950708152703, 0.011321425163965513, 0.009513173404649923, 0.012034954346008693, 0.0075649558956571005, 0.005983461431067319, 0.002352889026806826, 0.014264464758642278, 0.013867491767127864, 0.020811544799726887, 0.010649334293756111, 0.010987149405423407, 0.0016013200038222926, 0.007917301470563699, 0.00944265161484676, 0.010242268924089293, 0.015760559632993845, 0.01332314808268688, 0.011601147015340728, 0.009278209817378128, 0.02561630190717479, 0.010764805648372141, 0.010009038550711198, 0.021513922436260626, 0.002662434540212842, 0.010885578271187815, 0.004976748855254721, 0.012180846041250562, 0.011116946637204597, 0.007805234102414247, 0.015912742350539046, 0.009044733833884667, 0.010350188444304534, 0.005501189538980812, 0.013159332527035246, 0.010858261584869066, 0.0045728164936392835, 0.011350198342258585, 0.00098968380159947, 0.011498698810812183, 0.010339360035150753, 0.010582200214576043, 0.0005090338503854197, 0.00991881985807468], 0.367896468146942)

We can now plot the various elements (most of the code below is actually laying out the sub-panels for the plot):

Code for the figure

f = Figure(; size = (700, 700))
ax = Axis(f[1, 1]; aspect = DataAspect())
hm = heatmap!(
    ax,
    (βl - mean(βl)) / std(βl);
    colormap = Reverse(:RdYlBu),
    colorrange = (-2, 2),
)
lines!(ax, place; color = :black)
Colorbar(
    f[1, 1],
    hm;
    label = "Normalized LCBD",
    alignmode = Inside(),
    width = Relative(0.4),
    valign = :bottom,
    halign = :left,
    tellheight = false,
    tellwidth = false,
    vertical = false,
)
tightlimits!(ax)
hidespines!(ax)
hidedecorations!(ax)
ax_inset = Axis(f[1, 1];
    width = Relative(0.36),
    height = Relative(0.1),
    halign = 1.0,
    valign = 0.9,
    xlabel = "Rel. SCBD")
hist!(ax_inset, βs ./ maximum(βs); color = :lightgrey, bins = 40)
xlims!(ax_inset, 0, 1)
hidespines!(ax_inset, :t, :r, :l)
hideydecorations!(ax_inset)
hidexdecorations!(ax_inset; ticks = false, ticklabels = false, label = false)
tightlimits!(ax_inset)