Writing spatial variables to rasters

Once users have simulated various continuous geospatial values such as friction coefficients or population sizes across a landscape, they will often desire a visual confirmation of these values. One approach is to export these quantities as a multiband raster file and utilize different software, such as Quetzal-CRUMBS or the raster package in R, for visualization.

To export a specific quantity, it is necessary to create a callable function that can be invoked at each location within the raster model during an arbitrary time interval. Each point of the time interval parameter will become a new band in the new raster file.

Input

#include "quetzal/geography.hpp"
#include <cassert>
#include <filesystem>
#include <iostream>
 
using namespace quetzal;
 
int main()
{
    using time_type = int;
    using raster_type = geography::raster<time_type>;
    using location_descriptor = typename raster_type::location_descriptor;
 
    // The raster has 10 bands that we will assign to 2001 ... 2011.
    std::vector<time_type> times(10);
    std::iota(times.begin(), times.end(), 2001);
 
    auto input_file = std::filesystem::current_path() / "data/bio1.tif";
 
    // Read the raster
    auto bio1 = raster_type::from_file(input_file, times);
 
    // Define the quantity we want to print (it could be any function of space and time)
    auto transform = [](location_descriptor x, time_type t) { return std::optional(2 * x + t); };
 
    // Output file name
    auto output_file = std::filesystem::current_path() / "my_quantity.tif";
 
    // Write only between these times
    bio1.to_geotiff(transform, 2001, 2005, output_file);
 
    std::cout << (std::filesystem::exists("my_quantity.tif") ? "success" : "error") << std::endl;
 
    // Clean up
    std::filesystem::remove_all(output_file);
}

Output

success