Dispersal on a spatial graph

Table of Contents

• Distance-Based Kernels
  • Model's Basics
  • Kernels from Spatial Graphs
  • Using Kernels with Spatial Graphs Edges
    • Default Initialization and Later Update
    • Direct Initialization

Distance-Based Kernels

A straightforward approach to determine the dispersal mode across a complete spatial graph involves correlating the geographic distance between two locations with the likelihood of dispersal between them. The objective of this section is to highlight this model choice by parametrizing a simple Dispersal Location Kernel (in the sense of Nathan et al. 2012). Additionally, we will calculate essential metrics, such as the distance between two coordinates, the dispersal probability between them, and the average dispersal distance anticipated under this distribution.

Model's Basics

The dispersal location kernel represents the statistical pattern of dispersal distances within a population. In this context, it essentially serves as the probability density function (pdf) that outlines the distribution of locations after dispersal in relation to the original location. The expressions have been adjusted to incorporate a scale parameter, \(a\), which is consistent with a distance unit, and a shape parameter, \(b\), that dictates the curvature of the distribution curve, emphasizing the influence of long-distance events.

For a source \((x_0,y_0)\), the dispersal location kernel denoted as \(k_L(r)\) provides the density of the probability of the dispersal end point in the 2D space. In this case, \(k_L(r)dA\) is the probability of a dispersal end point to be within a small 2D area \(dA\) around the location \((x,y)\). Since a probability is unitless and \(dA\) is an area, \(k_L(r)\) is expressed in per unit area in a 2D space.

Quetzal incorporates various kernel types available in the quetzal::demography::dispersal_kernel namespace, and streamlines compile-time dimensional analysis and units conversion using the mp-units library.

Input

#include "quetzal/quetzal.hpp"
#include <iostream>
#include <mp-units/ostream.h>
 
using namespace quetzal;
 
int main()
{
    // Shorten notation
    using namespace mp_units::si::unit_symbols; // SI system: km, m, s
    using quetzal::demography::dispersal_kernel::gaussian;
    using quetzal::geography::lonlat;
 
    // Dispersal kernels operate on distances between geographic coordinates
    constexpr auto source = lonlat(2.25, 48.9176);
    constexpr auto target = lonlat(2.25, 48.95);
    constexpr auto distance = source.great_circle_distance_to(target);
 
    // Used to parametrize the kernel
    constexpr auto param = gaussian<>::param_type{.a = 1000. * m};
 
    // Compute quantities
    auto p = gaussian<>::pdf(distance, param);
    auto d = gaussian<>::mean_dispersal_distance(param);
 
    std::cout << "Dispersal from " << source << " to " << target << " :\n"
              << "distance is " << distance << "\n"
              << "probability is " << p << "\n"
              << "mean dispersal distance is " << d << std::endl;
}

Output

Dispersal from (Lat: 48.9176, Lon: 2.25) to (Lat: 48.95, Lon: 2.25) :
distance is 3602.72 m
probability is 7.34343e-13 1/m²
mean dispersal distance is 886.227 m

---

Kernels from Spatial Graphs

Users can create a spatial graph from landscape data, customizing assumptions about connectivity and topology. By defining a dispersal kernel and calculating dispersal probabilities for each edge based on distance metrics, users can uncover insights into potential dispersal patterns within the spatial graph. This approach offers a powerful means for exploring and understanding spatial dynamics and connectivity, facilitating deeper analysis of ecological or geographic phenomena.

Here we show how to simply print out the probabilities computed at each edge of the spatial graph.

Input

#include "quetzal/quetzal.hpp"
 
#include <cassert>
#include <filesystem>
#include <ranges>
 
namespace geo = quetzal::geography;
using namespace mp_units::si::unit_symbols; // SI system: km, m, s
 
int main()
{
    auto file1 = std::filesystem::current_path() / "data/bio1.tif";
    auto file2 = std::filesystem::current_path() / "data/bio12.tif";
 
    // The raster have 10 bands that we will assign to 2001 ... 2010.
    std::vector<int> times(10);
    std::iota(times.begin(), times.end(), 2001);
 
    // Initialize the landscape: for each var a key and a file, for all a time series.
    using landscape_type = geo::landscape<>;
    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
    // Our graph will not store any useful information
    using vertex_info = geo::no_property;
    using edge_info = geo::no_property;
 
    // We convert the grid to a graph, passing our assumptions
    auto graph = geo::from_grid(land, vertex_info(), edge_info(), geo::connect_fully(), geo::isotropy(),
                                 geo::mirror());
 
    // Define a helper function to compute distance on Earth between two location_descriptors
    auto sphere_distance = [&](auto point1, auto point2)
    {
        return land.to_latlon(point1).great_circle_distance_to(land.to_latlon(point2));
    };
 
    // Define the type of distance-based kernel to use
    using kernel = quetzal::demography::dispersal_kernel::exponential_power<>;
 
    // Transform the edges of the graph to a vector of probability to disperse from source to target
    auto probabilities = graph.edges()
        | std::views::transform( [&](auto const& e){ return sphere_distance( graph.source( e ), graph.target( e ) );} )
        | std::views::transform( [&](auto distance){ return kernel::pdf( distance, { .a=200.*km , .b=5.5 } ) ;} );
 
    // Print the result
    for (auto const& i : probabilities)
        std::cout << i << '\n';
}

Output

5.89588e-23 1/m²
0 1/m²
5.89588e-23 1/m²
5.58082e-132 1/m²
0 1/m²
0 1/m²
0 1/m²
5.58082e-132 1/m²
0 1/m²
1.69142e-30 1/m²
0 1/m²
0 1/m²
5.58082e-132 1/m²
0 1/m²
1.69142e-30 1/m²
0 1/m²
0 1/m²
0 1/m²
5.58082e-132 1/m²
0 1/m²
0 1/m²
0 1/m²
0 1/m²
0 1/m²
0 1/m²
5.58082e-132 1/m²
0 1/m²
1.22468e-40 1/m²
0 1/m²
0 1/m²
0 1/m²
0 1/m²
0 1/m²
5.58082e-132 1/m²
0 1/m²
1.22468e-40 1/m²

---

Using Kernels with Spatial Graphs Edges

In a simulation context, it can be more efficient to compute the dispersal probabilities along the graph's edges once, especially if these probabilities remain constant over time. This saves repeated computations during each simulation step.

There are two main approaches for setting edge properties:

1. Default Initialization and Later Update: You can create the graph with default edge properties and then update those properties later by iterating over the edges. It's simple to implement, but you pay the cost of iterating twice over the edges.
2. Direct Initialization: If you want to initialize the edge properties directly, the property itself should be a user-defined class with a proper default constructor, allowing for efficient initialization during graph creation. This approach helps optimize performance, particularly in simulations with many edges or static dispersal probabilities.

Default Initialization and Later Update

Here, we explain how to create a weighted graph with default edge properties and update those properties later by iterating over the edges.

First, when constructing the graph, each edge is initialized with a default set of properties (e.g., weights, dispersal probabilities). This allows you to quickly set up the graph structure without needing to provide specific values for every edge upfront.

Once the graph is created, you can update the edge properties by looping through each edge and assigning more accurate or context-specific values. This method is particularly useful when the final edge properties are derived from complex calculations or external data that may not be immediately available at graph creation.

By initializing the graph with default values, you maintain flexibility, allowing you to efficiently adjust edge properties later in the simulation or process without needing to recreate the entire graph structure.

Input

#include "quetzal/quetzal.hpp"
 
#include <cassert>
#include <filesystem>
#include <ranges>
 
namespace geo = quetzal::geography;
using namespace mp_units::si::unit_symbols; // SI system: km, m, s
 
int main()
{
    auto file1 = std::filesystem::current_path() / "data/bio1.tif";
    auto file2 = std::filesystem::current_path() / "data/bio12.tif";
 
    // The raster have 10 bands that we will assign to 2001 ... 2010.
    std::vector<int> times(10);
    std::iota(times.begin(), times.end(), 2001);
 
    // Initialize the landscape: for each var a key and a file, for all a time series.
    using landscape_type = geo::landscape<>;
    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
    // Define the type of distance-based kernel to use
    using kernel = quetzal::demography::dispersal_kernel::exponential_power<>;
 
    // Edges will store the result of the kernel's probability distribution function
    using vertex_info = geo::no_property;
    using edge_info = kernel::pdf_result_type;
 
    // We convert the grid to a graph, passing our assumptions
    auto graph = geo::from_grid(land, vertex_info(), edge_info(), geo::connect_fully(), geo::isotropy(),
                                 geo::mirror());
 
    // Define a helper function to compute distance on Earth between two location_descriptors
    auto sphere_distance = [&](auto point1, auto point2){
        return land.to_latlon(point1).great_circle_distance_to(land.to_latlon(point2));
    };
 
    // Transform the edges of the graph to a vector of probability to disperse from source to target
    for (const auto& e : graph.edges() ) {
        auto r = sphere_distance( graph.source( e ), graph.target( e ) );
        graph[e] = kernel::pdf( r, { .a=200.*km , .b=5.5 } );
    }
 
    // Print the result
    for (const auto& e : graph.edges() ) {
        std::cout << e << " : " << graph[e] << '\n';
    }
}

Output

(1,0) : 5.89588e-23 1/m²
(2,0) : 0 1/m²
(2,1) : 5.89588e-23 1/m²
(3,0) : 5.58082e-132 1/m²
(3,1) : 0 1/m²
(3,2) : 0 1/m²
(4,0) : 0 1/m²
(4,1) : 5.58082e-132 1/m²
(4,2) : 0 1/m²
(4,3) : 1.69142e-30 1/m²
(5,0) : 0 1/m²
(5,1) : 0 1/m²
(5,2) : 5.58082e-132 1/m²
(5,3) : 0 1/m²
(5,4) : 1.69142e-30 1/m²
(6,0) : 0 1/m²
(6,1) : 0 1/m²
(6,2) : 0 1/m²
(6,3) : 5.58082e-132 1/m²
(6,4) : 0 1/m²
(6,5) : 0 1/m²
(7,0) : 0 1/m²
(7,1) : 0 1/m²
(7,2) : 0 1/m²
(7,3) : 0 1/m²
(7,4) : 5.58082e-132 1/m²
(7,5) : 0 1/m²
(7,6) : 1.22468e-40 1/m²
(8,0) : 0 1/m²
(8,1) : 0 1/m²
(8,2) : 0 1/m²
(8,3) : 0 1/m²
(8,4) : 0 1/m²
(8,5) : 5.58082e-132 1/m²
(8,6) : 0 1/m²
(8,7) : 1.22468e-40 1/m²

Direct Initialization

Here, we demonstrate how to directly initialize the edge properties during graph creation. This approach is particularly useful when working with large graphs, where manually updating each edge later would be inefficient. By initializing the properties in place, you can ensure that each edge has the correct values from the start.

This method requires defining a user-specific Edge Property class, which includes a constructor to allow for proper instantiation. Doing so enables you to efficiently assign properties such as dispersal probabilities or weights directly to the graph’s edges during the initial setup, streamlining the overall process.

Input

#include "quetzal/quetzal.hpp"
 
#include <cassert>
#include <filesystem>
#include <ranges>
 
namespace geo = quetzal::geography;
using namespace mp_units::si::unit_symbols; // SI system: km, m, s
 
using kernel_type = quetzal::demography::dispersal_kernel::exponential_power<>;
 
struct edge_info
{
    using probability_type = kernel_type::pdf_result_type;
    probability_type weight;
    auto compute(auto source, auto target, const auto & grid)
    {
        auto r = grid.to_latlon(source).great_circle_distance_to(grid.to_latlon(target));
        return kernel_type::pdf( r, { .a=200.*km , .b=5.5 } );
    }
    edge_info(auto source, auto target, const auto & grid): weight(compute(source, target, grid)) {}
    edge_info() = default;
};
 
int main()
{
    auto file1 = std::filesystem::current_path() / "data/bio1.tif";
    auto file2 = std::filesystem::current_path() / "data/bio12.tif";
 
    // The raster have 10 bands that we will assign to 2001 ... 2010.
    std::vector<int> times(10);
    std::iota(times.begin(), times.end(), 2001);
 
    // Initialize the landscape: for each var a key and a file, for all a time series.
    using landscape_type = geo::landscape<>;
    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
    // Edges will store the result of the kernel's probability distribution function
    using vertex_info = geo::no_property;
 
    // We convert the grid to a graph, passing our assumptions
    auto graph = geo::from_grid(land, 
                                vertex_info(), 
                                edge_info(), 
                                geo::connect_fully(), 
                                geo::isotropy(),
                                geo::mirror());
 
    // Print the result
    for (const auto& e : graph.edges() ) {
        std::cout << e << " : " << graph[e].weight << '\n';
    }
}

Output

(1,0) : 0 1/m²
(2,0) : 0 1/m²
(2,1) : 0 1/m²
(3,0) : 0 1/m²
(3,1) : 0 1/m²
(3,2) : 0 1/m²
(4,0) : 0 1/m²
(4,1) : 0 1/m²
(4,2) : 0 1/m²
(4,3) : 0 1/m²
(5,0) : 0 1/m²
(5,1) : 0 1/m²
(5,2) : 0 1/m²
(5,3) : 0 1/m²
(5,4) : 0 1/m²
(6,0) : 0 1/m²
(6,1) : 0 1/m²
(6,2) : 0 1/m²
(6,3) : 0 1/m²
(6,4) : 0 1/m²
(6,5) : 0 1/m²
(7,0) : 0 1/m²
(7,1) : 0 1/m²
(7,2) : 0 1/m²
(7,3) : 0 1/m²
(7,4) : 0 1/m²
(7,5) : 0 1/m²
(7,6) : 0 1/m²
(8,0) : 0 1/m²
(8,1) : 0 1/m²
(8,2) : 0 1/m²
(8,3) : 0 1/m²
(8,4) : 0 1/m²
(8,5) : 0 1/m²
(8,6) : 0 1/m²
(8,7) : 0 1/m²