Kriging

This crate supports ordinary and universal kriging. Additionaly, it provides tools to create and fit variograms.

For now, all data must have EPSG:4326 (WGS84) as Coordinate Reference System (CRS).

lib.rs:

Kriging

This crate supports ordinary and universal kriging. Additionaly, it provides tools to create and fit variograms.

Ordinary kriging example

use std::path::PathBuf;
use polars::prelude::*;

use kriging::variogram::*;
use kriging::kriging::*;
use kriging::raster::*;

fn main() {
    let df = CsvReadOptions::default()
        .with_infer_schema_length(None)
        .with_has_header(true)
        .try_into_reader_with_file_path(Some("data.csv"))?
        .finish()?;

    let longitudes = df.column("lon")?.f64()?;
    let latitudes = df.column("lat")?.f64()?;
    let temps = df.column("temperature")?.f64()?;

    let longitudes = longitudes.to_vec().into_iter().filter_map(|x| x).collect();
    let latitudes = latitudes.to_vec().into_iter().filter_map(|x| x).collect();
    let temps = temps.to_vec().into_iter().filter_map(|x| x).collect();

    // Create an empirical variogram
    let empirical_variogram = EmpiricalVariogram::new(&longitudes, &latitudes, &temps);

    // Calculate the max values (optional for a first variogram estimation)
    let max_x_data = empirical_variogram.get_max_distance();
    let max_y_data = empirical_variogram.get_max_semivariance();

    // Define initial parameters for the model
    let mut variogram = Variogram::new(
        0.0, 
        max_y_data, 
        max_x_data / 3.0, 
        VariogramModel::Spherical
    );

    // Perform gradient descent to fit variogram
    variogram.fit(
        &empirical_variogram, 
        1e-2, 
        1000, 
        1e-6
    );

    println!("\nFitted Parameters:");
    println!("Nugget: {:.4}", variogram.nugget);
    println!("Sill: {:.4}", variogram.sill);
    println!("Range: {:.2}", variogram.range);

    // Plot the points (empirical variogram) and the fitted curve (variogram)
    let mut semivariogram_file_path = PathBuf::from("semivariogram.png");
    plot(&variogram, &empirical_variogram, &semivariogram_file_path);

    // Prepare coords and values for ordinary kriging
    let known_coords: Vec<(f64, f64)> = longitudes
        .into_iter()
        .zip(latitudes)
        .collect();

    let known_values: Vec<f64> = temps;

    // Read prediction points from TIF file
    let mut template_raster_file_path = PathBuf::from("sky_view_factor.tif");
    let prediction_points = read_raster_coords(&template_raster_file_path)?;

    // Perform ordinary kriging
    let predictions = ordinary_kriging(&known_coords, &known_values, &prediction_points, &variogram);

    // Export prediction
    let mut output_raster_file_path = PathBuf::from("interpolation.tif");
    export_predictions_to_raster(
        &template_raster_file_path,
        &output_raster_file_path, 
        &predictions
    )?;
}

Universal kriging example

use std::path::PathBuf;
use polars::prelude::*;

use kriging::variogram::*;
use kriging::kriging::*;
use kriging::raster::*;

fn main() {
    // Read observations
    let df = CsvReadOptions::default()
        .with_infer_schema_length(None)
        .with_has_header(true)
        .try_into_reader_with_file_path(Some("data.csv"))?
        .finish()?;

    let longitudes = df.column("lon")?.f64()?;
    let latitudes = df.column("lat")?.f64()?;
    let temps = df.column("temperature")?.f64()?;

    let longitudes = longitudes.to_vec().into_iter().filter_map(|x| x).collect();
    let latitudes = latitudes.to_vec().into_iter().filter_map(|x| x).collect();
    let temps = temps.to_vec().into_iter().filter_map(|x| x).collect();

    // This variogram will be fitted with the residuals inside universal_kriging()
    let initial_estimation_variogram = Variogram::new(
        0.1, 
        0.4, 
        1000, 
        VariogramModel::Spherical
    );

    // Prepare coords and values for ordinary kriging
    let known_coords: Vec<(f64, f64)> = longitudes
        .into_iter()
        .zip(latitudes)
        .collect();

    let known_values: Vec<f64> = temps;

    // Read prediction points from TIF file
    let mut template_raster_file_path = PathBuf::from("sky_view_factor.tif");
    let prediction_points = read_raster_coords(&template_raster_file_path)?;

    // Perform ordinary kriging
    // Read covariates
    let gli_path = PathBuf::from(env::var("CARGO_MANIFEST_DIR").unwrap()).join("data/gli.tif");
    let svf_path = PathBuf::from(env::var("CARGO_MANIFEST_DIR").unwrap()).join("data/svf.tif");
    let drift = RasterDrift::new(&[gli_path.to_str().unwrap(), svf_path.to_str().unwrap()]);

    // Perform ordinary kriging
    let (predictions, variances) = universal_kriging(
        &known_coords, &known_values, &prediction_points, drift, initial_estimation_variogram
    );

    // Export prediction
    let mut output_raster_file_path = PathBuf::from("interpolation.tif");
    export_predictions_to_raster(
        &template_raster_file_path,
        &output_raster_file_path,
        &predictions
    )?;

    // Export variances
    let mut output_raster_file_path = PathBuf::from(&root_dir);
    output_raster_file_path.push("variances.tif");
    export_predictions_to_raster(
        &template_raster_file_path,
        &output_raster_file_path, 
        &variances
    )?;
}