NeutralLandscapes.jl Walkthrough

This document provides a tutorial in using NeutralLandscapes.jl, the julia library for simulating neutral landscapes. Neutral landscapes have become a common tool in landscape ecology, primarily as a null model of spatial data (Gardner et al. 1987). Here we present a package for neutral-landscape generation in julia, which both is considerably faster than both existing R and python packages, but also implements novel methods of simulating temporal change in landscapes with fixed statistical properties.

Gardner, Robert H., Bruce T. Milne, Monica G. Turnei, and Robert V. O’Neill. 1987. “Neutral Models for the Analysis of Broad-Scale Landscape Pattern.” Landscape Ecology 1 (1): 19–28. https://doi.org/10.1007/BF02275262.

In this walkthrough we’ll cover (1) installing julia and NeutralLandscapes.jl, (2) using different neutral landscape generators using the rand function, (3) simulating temporal change in neutral landscapes (a feature unique to NL.jl), and (4) simulating the occurrence of a species tracking its niche over time, producing the following output.

The final product at the of this walkthrough

1. Getting Started

In this section we’ll install julia (if you don’t already have it) and the NeutralLandscapes.jl package.

Installing julia

First, if you don’t have julia installed on your machine, install it from https://julialang.org/downloads/. (Note if you are using an M1 Mac, v1.8 is recommended. Otherwise either v1.7 or v1.8 are both fine).

There are two ways to interact with julia. You can either open the Read-Evaluate-Print-Loop (REPL) by trying julia in a terminal, e.g.

user@machine> julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.0 (2021-11-30)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia>

The REPL acts just like a console in R. The second method is to use julia to run a .jl script, e.g.

user@machine> julia script.jl

Installing NeutralLandscapes.jl

The easiest way to install NeutralLandscapes.jl, or any julia package is through the REPL. Open the REPL by typing julia in a terminal . The julia shell can operate in a few different modes. By default, the shell runs in julia mode, where the prompt looks like

julia>

You can enter the package-management mode by press the square-bracket key (]), which changes to prompt to the package mode:

(v1.7) pkg>

To install NeutralLandscapes, we run

(v1.7) pkg> add NeutralLandscapes

and then once done installing, hit backspace to exit Pkg mode and to return to the julia console. Read more about the Julia REPL modes here and Pkg, the package-management system for julia, here.

2. Making Neutral Landscapes

In this section we’ll cover the syntax for generating neutral landscapes, starting with a simple gradient and then showing the various other types of landscape types we can generate.

Simple example

We start by loading the package using

using NeutralLandscapes

Now let’s create our first neutral landscape. We’ll create a linear gradient using the PlanarGradient generator. We can search for the documentation by entering the help mode of the REPL by pressing ?, or by using the ? function.

?(PlanarGradient)

search: PlanarGradient

PlanarGradient <: NeutralLandscapeMaker

PlanarGradient(; direction=360rand())
PlanarGradient(direction)

This type is used to generate a planar gradient landscape, where values change as a bilinear function of the x and y coordinates. The direction is expressed as a floating point value, which will be in [0,360]. The inner constructor takes the mod of the value passed and 360, so that a value that is out of the correct interval will be corrected.

Using the above process we can understand what arguments to each neutral landscape generator are, and what they do.

Using rand to generate landscapes

To use PlanarGradient to generate a landscape, we use the rand function. The rand function is used for all sorts of random generation in julia. For example, just rand() generates a random number between 0 and 1, e.g.

rand()

0.46110243791553074

Although rand can be called with different arguments to run a different version of the function that corresponds to the type of the input arguments. For example, using the Distributions package allows us to draw random samples from various distributions, e.g. a Poisson distribution with mean 7.

using Distributions
rand(Poisson(7))

8

or a set of 10 independent draws from a standard Normal distribution.

rand(Normal(0,1), 10)

10-element Vector{Float64}:
 -1.0750868899665584
  1.4726498920951105
  0.08334169367764765
  1.001596634832882
  0.41825145919788725
  0.3409662777991979
  0.968971026996396
 -0.7741079564132628
  1.9897432432718793
 -0.2856305633137191

This is possible through a feature called multiple dispatch, where a function with the same name (rand, in this case) will behave differently depending on the arguments passed to it (similar to overloading in C/C++).

The syntax for generating random neutral landscapes is rand(NEUTRAL_LANDSCAPE_GENERATOR, X_SIZE, Y_SIZE). For example, to generate a PlanarGradient landscape with dimensions 200 pixels x 100 pixels, we run

neutrallandscape = rand(PlanarGradient(), (200,100));

which we plot with

using Plots
heatmap(neutrallandscape)

More complicated NeutralLandscapeMakers

NeutralLandscapes.jl supports a large number of methods for generating neutral landscapes. Feature parity is with the Python package for neutral landscapes (Etherington, Holland, and O’Sullivan 2015).

Etherington, Thomas R., E. Penelope Holland, and David O’Sullivan. 2015. “NLMpy: A Python Software Package for the Creation of Neutral Landscape Models Within a General Numerical Framework.” Edited by Timothée Poisot. Methods in Ecology and Evolution 6 (2): 164–68. https://doi.org/10.1111/2041-210X.12308.

Gradients

The simplest models are gradients. We implement NoGradient, which is just random, as well as PlanarGradient (a linear gradient), EdgeGradient (a symetric gradient around a line), and DistanceGradient (gradient as a function of distance from a given point).

siz = (100,100)
random = rand(NoGradient(), siz)
planar = rand(PlanarGradient(), siz)
edge = rand(EdgeGradient(), siz)

obj = falses(siz)
obj[10:25, 10:25] .= true
dist = rand(DistanceGradient(findall(vec(obj))), siz)

plot(heatmap(random),
heatmap(edge),
heatmap(dist),
heatmap(planar))

Clustering

Clustering methods (RectangularCluster, NearestNeighborElement, and NearestNeighborCluster) produce discrete boundaries. These are useful for simulating landcover data, or other data with hard boundaries.

rect = rand(RectangularCluster(4, 8), siz)
nne = rand(NearestNeighborElement(200), siz)
nnc = rand(NearestNeighborCluster(0.4), siz)
plot(heatmap(rect), heatmap(nne), heatmap(nnc))

Spatially Autocorrelation

Fractal Brownian motion has a long history in ecology (Milne 1992). NL.jl implements two similar but slightly different method for FBM generation, DiamondSquare and MidpointDisplacement. Each take a parameter H, between 0 and 1, which describes the level of autocorrelation.

Milne, Bruce T. 1992. “Spatial Aggregation and Neutral Models in Fractal Landscapes.” The American Naturalist 139 (1): 32–57. https://doi.org/10.1086/285312.

mpd = rand(MidpointDisplacement(0.75), siz)
ds = rand(DiamondSquare(0.25), siz)
plot(heatmap(mpd), heatmap(ds), size=(700,300))

Classification

Using the classify function, we can bin values into a set of discrete categories. This can be useful for simulating discrete patches in otherwise uninhabitable landscape.

clas1 = classify(ds, [1-0.5, 0.5])
clas2 = classify(ds, [1-0.75, 0.75])
clas3 = classify(mpd, ones(3))
clas4 = classify(mpd, [0.2i for i in 0:5])
plot(heatmap(clas1), heatmap(clas2), heatmap(clas3), heatmap(clas4))

3. Landscapes that change over time

Now we get to the novel part in NeutralLandscapes.jl, which is the ability to simulate change using 3 different methods.

function makeplot(updater)
    init = rand(PlanarGradient(π/2), 100, 100)
    seq = update(updater, init, 30)
    print(length(seq))
    plot(heatmap(init, cbar=:none), heatmap(seq[10], cbar=:none), heatmap(seq[20], cbar=:none), heatmap(seq[30], cbar=:none), size=(700,200), layout=(1,4))
end

makeplot (generic function with 1 method)

Temporally variable change

?(TemporallyVariableUpdater)

search: TemporallyVariableUpdater

TemporallyVariableUpdater{D,S} <: NeutralLandscapeUpdater

A NeutralLandscapeUpdater that has a prescribed level of temporal variation (variability) and rate of change (rate), but no spatial correlation in where change is distributed.

tvu = TemporallyVariableUpdater()
makeplot(tvu)

30

Spatially autocorrelated change

?(SpatiallyAutocorrelatedUpdater)

search: SpatiallyAutocorrelatedUpdater SpatiotemporallyAutocorrelatedUpdater

SpatiallyAutocorrelatedUpdater{SU,R,V}

A NeutralLandscapeUpdater that has a prescribed level of spatial variation (variability) and rate of change (rate), and where the spatial distribution of this change is proportional to a neutral landscape generated with spatialupdater at every time step.

TODO: make it possible to fix a given spatial updater at each timestep.

sau = SpatiallyAutocorrelatedUpdater(spatialupdater=MidpointDisplacement(0.5))
makeplot(sau)

30

Spatiotemporally autocorrelated change

?(SpatiotemporallyAutocorrelatedUpdater)

search: SpatiotemporallyAutocorrelatedUpdater

SpatiotemporallyAutocorrelatedUpdater{SU,R,V}

A NeutralLandscapeUpdater that has a prescribed level of spatial and temporal variation (variability) and rate of change (rate), and where the spatial distribution of this change is proportional to a neutral landscape generated with spatialupdater at every time step.

TODO: perhaps spatial and temporal should each have their own variability param

stau = SpatiotemporallyAutocorrelatedUpdater(spatialupdater=MidpointDisplacement(0.5), variability=1.0)
makeplot(stau)

30

4. Simulating occurrence data across Quebec

Now we’ll use NL.jl to simulate the occurrence of a species tracking its niche over time across Quebec.

Masking a NeutralLandscape across Quebec

First we’ll show how NeutralLandscapes.jl can integrate with SimpleSDMLayers.jl, the julia package for species distribution modeling.

using SimpleSDMLayers
using Plots
ENV["SDMLAYERS_PATH"] = "/home/michael/data/";

First we use SimpleSDMLayers.jl to download BioClim data for Quebec

quebec = SimpleSDMPredictor(WorldClim, BioClim; left=-90., right=-55., top=63., bottom=45.)

SDM predictor → 108×210 grid with 15058 Float32-valued cells
  Latitudes 45.0 ⇢ 63.0
  Longitudes    -90.0 ⇢ -55.0

Now we’ll construst a mask to use with NL.jl

qcmask = fill(true, size(quebec))
qcmask[findall(isnothing, quebec.grid)] .= false;

And now we’ll plot several sample NeutralLandscapes masked over Quebec.

pltsettings = (cbar=:none, frame=:box)
plot(
    heatmap(rand(MidpointDisplacement(0.8), size(quebec), mask=qcmask); pltsettings...),
    heatmap(rand(PlanarGradient(), size(quebec), mask=qcmask); pltsettings...),
    heatmap(rand(PerlinNoise((4,4)), size(quebec), mask=qcmask); pltsettings...),
    heatmap(rand(NearestNeighborElement(30), size(quebec), mask=qcmask); pltsettings...),
    dpi=400
)

Defining a niche

Now we’re going to define a niche function, which returns the occurrence probability of a species with a niche that is a Gaussian with mean nichemean and standard deviation nichewidth.

niche(env; nichemean = 0.7, nichewidth = 0.3) = exp(-1*((env-nichemean)/nichewidth)^2)

niche (generic function with 1 method)

And now we’ll write a function to visualize these niches

function nicheplot(; nichemean = 0.7, nichewidth = 0.3)
    envrange = 0:10^-3:1
    plot(
        envrange,
        label=:none,
        lw=3,
        title="Niche center: $nichemean, Niche width: $nichewidth",
        [niche(x; nichemean=nichemean, nichewidth=nichewidth) for x in envrange],
        xlabel="Environmental value",
        ylabel="Occurrence Probability",
        titlefontsize=9,
        labelfontsize=9,
        frame=:box
    )
end

nicheplot (generic function with 1 method)

Now let’s plot some sample niches to get a visual sense of how they work.

envrange = 0:10^-3:1
plot(
    nicheplot(;nichemean=0.5, nichewidth=0.2),
    nicheplot(;nichemean=0.7, nichewidth=0.1),
    nicheplot(;nichemean=0.1, nichewidth=0.05),
    nicheplot(;nichemean=0.5, nichewidth=0.5),
    size=(600,600)
)

Simulating occurrence

Finally, we’ll write the occurrence function which returns a set of numpoints coordinates drawn with probabilities proportional to the occurrence probability provided by niche at each cell in the landscape.

function occurrence(landscape, qcmask, numpoints = 50; kwargs...)
    probabilitymap = map(e-> isnan(e) ? 0 : niche(e; kwargs...), landscape)
    normalizedprob = probabilitymap ./ sum(probabilitymap)

    occmat = zeros(size(landscape))
    occind = rand(Categorical(vec(normalizedprob)), numpoints)
    occmat[occind] .= 1
    occpoints = findall(x->x==1, occmat)
    longs, lats = [x[2] for x in occpoints], [x[1] for x in occpoints]
end

occurrence (generic function with 2 methods)

Now we set an initial condition that is moderately spatially autocorrelated and that is generally higher in the south than the north (think temperature).

grad = rand(PlanarGradient(0), size(quebec), mask=qcmask)
mpd = rand(MidpointDisplacement(0.3), size(quebec), mask=qcmask)

init = blend(grad,mpd)
heatmap(init)

Now we simulate the temporally changing landscape

svu = SpatiotemporallyAutocorrelatedUpdater(
    rate = 0.03,
    variability = 0.03,
    spatialupdater = MidpointDisplacement(0.5)
)
seq = normalize(update(svu, init, 30));

and finally simulate occurrence and plot, to produce the final output.

anim = @animate for (t,s) in enumerate(seq)
    longs, lats = occurrence(s, qcmask, 30; nichemean=0.5, nichewidth=0.05)
    heatmap(s, dpi=150, title = "Year $t", frame=:box, ticks=:none, aspectratio=1, size=(400,300), xlim=(0,220), ylim=(0,110), clim=(0,1))
    scatter!(longs, lats, c=:dodgerblue, msw=0, label="")
end
gif(anim, "result.gif", fps=5)

┌ Info: Saved animation to 
│   fn = /home/michael/code/NLIntroductions/result.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114