kriging 0.0.1

This crate supports ordinary and universal kriging. Additionaly, it provides tools to create and fit variograms.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

use std::env;
use std::path::PathBuf;
use polars::prelude::*;
use std::error::Error as StdError;
use variogram::*;
use kriging::*;
use raster::*;

mod variogram;
mod kriging;
mod raster;

fn main() -> Result<(), Box<dyn StdError>> {
    let root_dir = env::var("CARGO_MANIFEST_DIR")?;
    let mut csv_file_path = PathBuf::from(&root_dir);
    csv_file_path.push("data/data_netatmo.csv");

    let df = CsvReadOptions::default()
            .with_infer_schema_length(None)
            .with_has_header(true)
            .try_into_reader_with_file_path(Some(csv_file_path.into()))?
            .finish()?;

    println!("{:?}", df);

    let longitudes = df.column("lon")?.f64()?;
    let latitudes = df.column("lat")?.f64()?;
    let temps = df.column("temp_diff")?.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
    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
    println!("\nRunning gradient descent...");
    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(&root_dir);
    semivariogram_file_path.push("data/semivariogram.png");
    plot(&variogram, &empirical_variogram, &semivariogram_file_path);

    // Perform 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(&root_dir);
    template_raster_file_path.push("data/svf.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);

    // Perform universal kriging
    // Initialize and fit theoretical model (Spherical as example)
    let initial_estimation_variogram = Variogram::new(
        0.0,                   // nugget
        max_y_data,              // sill
        max_x_data / 3.0,       // range
        VariogramModel::Spherical,
    );

    // 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()]);

    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(&root_dir);
    output_raster_file_path.push("data/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("data/variances.tif");
    export_predictions_to_raster(
        &template_raster_file_path,
        &output_raster_file_path, 
        &variances
    )?;
    
    Ok(())
}