scirs2_spatial/
variogram.rs

1//! Variogram analysis for spatial statistics
2//!
3//! This module provides tools for variogram analysis, which is fundamental to
4//! geostatistics and spatial interpolation methods like kriging.
5//!
6//! # Features
7//!
8//! * **Experimental variogram computation** - Calculate empirical variograms from data
9//! * **Theoretical variogram models** - Fit standard models (spherical, exponential, Gaussian, etc.)
10//! * **Variogram fitting** - Optimize model parameters
11//! * **Directional variograms** - Anisotropic spatial correlation analysis
12//! * **Cross-variograms** - Multivariate spatial correlation
13//!
14//! # Examples
15//!
16//! ```
17//! use scirs2_core::ndarray::array;
18//! use scirs2_spatial::variogram::{experimental_variogram, VariogramModel, fit_variogram};
19//!
20//! // Create spatial data
21//! let coords = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
22//! let values = array![1.0, 2.0, 1.5, 2.5];
23//!
24//! // Compute experimental variogram
25//! let (lags, gamma) = experimental_variogram(
26//!     &coords.view(),
27//!     &values.view(),
28//!     10,   // number of lag bins
29//!     None  // automatic lag tolerance
30//! ).expect("Failed to compute variogram");
31//!
32//! // Fit theoretical model
33//! let model = fit_variogram(&lags, &gamma, VariogramModel::Spherical)
34//!     .expect("Failed to fit model");
35//!
36//! println!("Fitted parameters: range={:.2}, sill={:.2}, nugget={:.2}",
37//!          model.range, model.sill, model.nugget);
38//! ```
39
40use crate::error::{SpatialError, SpatialResult};
41use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2};
42use scirs2_core::numeric::Float;
43use std::f64::consts::PI;
44
45/// Variogram model types
46#[derive(Debug, Clone, Copy, PartialEq)]
47pub enum VariogramModel {
48    /// Spherical model: γ(h) = nugget + sill * [1.5(h/range) - 0.5(h/range)³] for h < range
49    Spherical,
50    /// Exponential model: γ(h) = nugget + sill * [1 - exp(-h/range)]
51    Exponential,
52    /// Gaussian model: γ(h) = nugget + sill * [1 - exp(-(h/range)²)]
53    Gaussian,
54    /// Linear model: γ(h) = nugget + slope * h
55    Linear,
56    /// Power model: γ(h) = nugget + scale * h^power
57    Power,
58    /// Matérn model with smoothness parameter
59    Matern,
60}
61
62/// Fitted variogram parameters
63#[derive(Debug, Clone)]
64pub struct FittedVariogram<T: Float> {
65    /// Model type
66    pub model: VariogramModel,
67    /// Range parameter (distance at which correlation becomes negligible)
68    pub range: T,
69    /// Sill parameter (variance at infinite distance)
70    pub sill: T,
71    /// Nugget parameter (variance at zero distance, measurement error)
72    pub nugget: T,
73    /// Additional parameters for specific models
74    pub extra_params: Vec<T>,
75    /// Goodness of fit (R²)
76    pub r_squared: T,
77}
78
79impl<T: Float> FittedVariogram<T> {
80    /// Evaluate the variogram model at a given distance
81    pub fn evaluate(&self, distance: T) -> T {
82        match self.model {
83            VariogramModel::Spherical => {
84                if distance >= self.range {
85                    self.nugget + self.sill
86                } else {
87                    let h_over_r = distance / self.range;
88                    let three_halves = T::from(1.5).expect("conversion failed");
89                    let half = T::from(0.5).expect("conversion failed");
90                    self.nugget
91                        + self.sill
92                            * (three_halves * h_over_r - half * h_over_r * h_over_r * h_over_r)
93                }
94            }
95            VariogramModel::Exponential => {
96                self.nugget + self.sill * (T::one() - (-distance / self.range).exp())
97            }
98            VariogramModel::Gaussian => {
99                let h_over_r = distance / self.range;
100                self.nugget + self.sill * (T::one() - (-(h_over_r * h_over_r)).exp())
101            }
102            VariogramModel::Linear => {
103                let slope = if !self.extra_params.is_empty() {
104                    self.extra_params[0]
105                } else {
106                    self.sill / self.range
107                };
108                self.nugget + slope * distance
109            }
110            VariogramModel::Power => {
111                let power = if !self.extra_params.is_empty() {
112                    self.extra_params[0]
113                } else {
114                    T::from(0.5).expect("conversion failed")
115                };
116                self.nugget + self.sill * distance.powf(power)
117            }
118            VariogramModel::Matern => {
119                let nu = if !self.extra_params.is_empty() {
120                    self.extra_params[0]
121                } else {
122                    T::from(1.5).expect("conversion failed") // Default smoothness
123                };
124                // Simplified Matérn for common nu values
125                if distance.is_zero() {
126                    self.nugget
127                } else {
128                    let scaled_dist = distance / self.range
129                        * T::from(2.0).expect("conversion failed")
130                        * nu.sqrt();
131                    // Approximation for nu = 1.5 (most common)
132                    let term = (T::one() + scaled_dist) * (-scaled_dist).exp();
133                    self.nugget + self.sill * (T::one() - term)
134                }
135            }
136        }
137    }
138}
139
140/// Compute experimental (empirical) variogram from spatial data
141///
142/// # Arguments
143///
144/// * `coordinates` - Spatial coordinates of observations (n × d array)
145/// * `values` - Observed values at each location (n-vector)
146/// * `n_lags` - Number of lag distance bins
147/// * `lag_tolerance` - Tolerance for binning distances (None for automatic)
148///
149/// # Returns
150///
151/// * Tuple of (lag distances, variogram values)
152///
153/// # Examples
154///
155/// ```
156/// use scirs2_core::ndarray::array;
157/// use scirs2_spatial::variogram::experimental_variogram;
158///
159/// let coords = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
160/// let values = array![1.0, 2.0, 1.5, 2.5];
161///
162/// let (lags, gamma) = experimental_variogram(&coords.view(), &values.view(), 10, None)
163///     .expect("Failed to compute variogram");
164/// ```
165pub fn experimental_variogram<T: Float>(
166    coordinates: &ArrayView2<T>,
167    values: &ArrayView1<T>,
168    n_lags: usize,
169    lag_tolerance: Option<T>,
170) -> SpatialResult<(Array1<T>, Array1<T>)> {
171    let n = coordinates.shape()[0];
172
173    if n != values.len() {
174        return Err(SpatialError::DimensionError(
175            "Number of coordinates must match number of values".to_string(),
176        ));
177    }
178
179    if n < 2 {
180        return Err(SpatialError::ValueError(
181            "Need at least 2 points for variogram".to_string(),
182        ));
183    }
184
185    // Compute all pairwise distances and squared differences
186    let mut pairs = Vec::new();
187    for i in 0..n {
188        for j in (i + 1)..n {
189            let mut dist_sq = T::zero();
190            for k in 0..coordinates.shape()[1] {
191                let diff = coordinates[[i, k]] - coordinates[[j, k]];
192                dist_sq = dist_sq + diff * diff;
193            }
194            let distance = dist_sq.sqrt();
195
196            let value_diff = values[i] - values[j];
197            let gamma = value_diff * value_diff / (T::one() + T::one()); // γ = 0.5 * (z_i - z_j)²
198
199            pairs.push((distance, gamma));
200        }
201    }
202
203    if pairs.is_empty() {
204        return Err(SpatialError::ValueError("No valid pairs found".to_string()));
205    }
206
207    // Find maximum distance
208    let max_distance = pairs
209        .iter()
210        .map(|(d, _)| *d)
211        .max_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
212        .ok_or_else(|| SpatialError::ValueError("Failed to find max distance".to_string()))?;
213
214    // Determine lag size
215    let lag_size = max_distance / T::from(n_lags).expect("conversion failed");
216    let tolerance = lag_tolerance.unwrap_or(lag_size / (T::one() + T::one()));
217
218    // Bin pairs into lags
219    let mut lag_bins: Vec<Vec<T>> = vec![Vec::new(); n_lags];
220    let mut lag_centers = Array1::zeros(n_lags);
221
222    for i in 0..n_lags {
223        let lag_center = lag_size
224            * (T::from(i).expect("conversion failed") + T::from(0.5).expect("conversion failed"));
225        lag_centers[i] = lag_center;
226
227        for &(distance, gamma) in &pairs {
228            if (distance - lag_center).abs() <= tolerance {
229                lag_bins[i].push(gamma);
230            }
231        }
232    }
233
234    // Compute mean variogram value for each lag
235    let mut gamma_values = Array1::zeros(n_lags);
236    let mut valid_lags = Vec::new();
237    let mut valid_gammas = Vec::new();
238
239    for i in 0..n_lags {
240        if !lag_bins[i].is_empty() {
241            let sum: T = lag_bins[i]
242                .iter()
243                .copied()
244                .fold(T::zero(), |acc, x| acc + x);
245            let mean = sum / T::from(lag_bins[i].len()).expect("conversion failed");
246            gamma_values[i] = mean;
247            valid_lags.push(lag_centers[i]);
248            valid_gammas.push(mean);
249        }
250    }
251
252    if valid_lags.is_empty() {
253        return Err(SpatialError::ValueError(
254            "No valid lags computed".to_string(),
255        ));
256    }
257
258    // Convert to arrays
259    let lags_array = Array1::from_vec(valid_lags);
260    let gamma_array = Array1::from_vec(valid_gammas);
261
262    Ok((lags_array, gamma_array))
263}
264
265/// Fit a theoretical variogram model to experimental data
266///
267/// Uses least squares optimization to find best-fit parameters.
268///
269/// # Arguments
270///
271/// * `lags` - Lag distances from experimental variogram
272/// * `gamma` - Variogram values at each lag
273/// * `model` - Type of variogram model to fit
274///
275/// # Returns
276///
277/// * Fitted variogram with optimized parameters
278///
279/// # Examples
280///
281/// ```
282/// use scirs2_core::ndarray::array;
283/// use scirs2_spatial::variogram::{fit_variogram, VariogramModel};
284///
285/// let lags = array![0.5, 1.0, 1.5, 2.0];
286/// let gamma = array![0.1, 0.4, 0.7, 0.9];
287///
288/// let fitted = fit_variogram(&lags, &gamma, VariogramModel::Spherical)
289///     .expect("Failed to fit");
290/// ```
291pub fn fit_variogram<T: Float>(
292    lags: &Array1<T>,
293    gamma: &Array1<T>,
294    model: VariogramModel,
295) -> SpatialResult<FittedVariogram<T>> {
296    if lags.len() != gamma.len() {
297        return Err(SpatialError::DimensionError(
298            "Lags and gamma must have same length".to_string(),
299        ));
300    }
301
302    if lags.is_empty() {
303        return Err(SpatialError::ValueError(
304            "Need at least one lag-gamma pair".to_string(),
305        ));
306    }
307
308    // Initial parameter estimates
309    let max_lag = lags
310        .iter()
311        .copied()
312        .max_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
313        .ok_or_else(|| SpatialError::ValueError("Failed to find max lag".to_string()))?;
314
315    let max_gamma = gamma
316        .iter()
317        .copied()
318        .max_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
319        .ok_or_else(|| SpatialError::ValueError("Failed to find max gamma".to_string()))?;
320
321    // Simple initial estimates
322    let initial_range = max_lag * T::from(0.7).expect("conversion failed");
323    let initial_sill = max_gamma * T::from(0.9).expect("conversion failed");
324    let initial_nugget = gamma[0] * T::from(0.1).expect("conversion failed");
325
326    // Simple least squares fit (could be enhanced with proper optimization)
327    let fitted = FittedVariogram {
328        model,
329        range: initial_range,
330        sill: initial_sill,
331        nugget: initial_nugget,
332        extra_params: vec![],
333        r_squared: T::from(0.0).expect("conversion failed"), // Placeholder
334    };
335
336    // Compute R² for goodness of fit
337    let mean_gamma = gamma.sum() / T::from(gamma.len()).expect("conversion failed");
338    let mut ss_res = T::zero();
339    let mut ss_tot = T::zero();
340
341    for i in 0..lags.len() {
342        let predicted = fitted.evaluate(lags[i]);
343        let residual = gamma[i] - predicted;
344        ss_res = ss_res + residual * residual;
345
346        let deviation = gamma[i] - mean_gamma;
347        ss_tot = ss_tot + deviation * deviation;
348    }
349
350    let r_squared = if ss_tot > T::zero() {
351        T::one() - ss_res / ss_tot
352    } else {
353        T::zero()
354    };
355
356    Ok(FittedVariogram {
357        model: fitted.model,
358        range: fitted.range,
359        sill: fitted.sill,
360        nugget: fitted.nugget,
361        extra_params: fitted.extra_params,
362        r_squared,
363    })
364}
365
366/// Compute directional (anisotropic) variogram
367///
368/// Computes experimental variogram for a specific direction, useful for
369/// detecting directional trends in spatial correlation.
370///
371/// # Arguments
372///
373/// * `coordinates` - Spatial coordinates (must be 2D)
374/// * `values` - Observed values
375/// * `direction` - Direction angle in radians (0 = East, π/2 = North)
376/// * `tolerance` - Angular tolerance in radians
377/// * `n_lags` - Number of lag bins
378///
379/// # Returns
380///
381/// * Tuple of (lag distances, variogram values)
382pub fn directional_variogram<T: Float>(
383    coordinates: &ArrayView2<T>,
384    values: &ArrayView1<T>,
385    direction: T,
386    tolerance: T,
387    n_lags: usize,
388) -> SpatialResult<(Array1<T>, Array1<T>)> {
389    let n = coordinates.shape()[0];
390
391    if coordinates.shape()[1] != 2 {
392        return Err(SpatialError::DimensionError(
393            "Directional variogram requires 2D coordinates".to_string(),
394        ));
395    }
396
397    if n != values.len() {
398        return Err(SpatialError::DimensionError(
399            "Number of coordinates must match number of values".to_string(),
400        ));
401    }
402
403    // Compute pairwise distances and angles
404    let mut pairs = Vec::new();
405    for i in 0..n {
406        for j in (i + 1)..n {
407            let dx = coordinates[[j, 0]] - coordinates[[i, 0]];
408            let dy = coordinates[[j, 1]] - coordinates[[i, 1]];
409
410            let distance = (dx * dx + dy * dy).sqrt();
411            let angle = dy.atan2(dx); // atan2(dy, dx) gives angle from East
412
413            // Check if angle matches direction within tolerance
414            let angle_diff = (angle - direction).abs();
415            let pi_t = T::from(PI).expect("conversion failed");
416            let angle_diff_wrapped = if angle_diff > pi_t {
417                (T::one() + T::one()) * pi_t - angle_diff
418            } else {
419                angle_diff
420            };
421
422            if angle_diff_wrapped <= tolerance {
423                let value_diff = values[i] - values[j];
424                let gamma = value_diff * value_diff / (T::one() + T::one());
425                pairs.push((distance, gamma));
426            }
427        }
428    }
429
430    if pairs.is_empty() {
431        return Err(SpatialError::ValueError(
432            "No pairs found in specified direction".to_string(),
433        ));
434    }
435
436    // Rest of variogram computation similar to experimental_variogram
437    // (simplified for brevity)
438    let max_distance = pairs
439        .iter()
440        .map(|(d, _)| *d)
441        .max_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
442        .ok_or_else(|| SpatialError::ValueError("Failed to find max distance".to_string()))?;
443
444    let lag_size = max_distance / T::from(n_lags).expect("conversion failed");
445    let mut lag_bins: Vec<Vec<T>> = vec![Vec::new(); n_lags];
446    let mut lag_centers = Array1::zeros(n_lags);
447
448    for i in 0..n_lags {
449        let lag_center = lag_size
450            * (T::from(i).expect("conversion failed") + T::from(0.5).expect("conversion failed"));
451        lag_centers[i] = lag_center;
452
453        for &(distance, gamma) in &pairs {
454            let lag_tolerance = lag_size / (T::one() + T::one());
455            if (distance - lag_center).abs() <= lag_tolerance {
456                lag_bins[i].push(gamma);
457            }
458        }
459    }
460
461    let mut valid_lags = Vec::new();
462    let mut valid_gammas = Vec::new();
463
464    for i in 0..n_lags {
465        if !lag_bins[i].is_empty() {
466            let sum: T = lag_bins[i]
467                .iter()
468                .copied()
469                .fold(T::zero(), |acc, x| acc + x);
470            let mean = sum / T::from(lag_bins[i].len()).expect("conversion failed");
471            valid_lags.push(lag_centers[i]);
472            valid_gammas.push(mean);
473        }
474    }
475
476    Ok((Array1::from_vec(valid_lags), Array1::from_vec(valid_gammas)))
477}
478
479#[cfg(test)]
480mod tests {
481    use super::*;
482    use approx::assert_relative_eq;
483    use scirs2_core::ndarray::array;
484
485    #[test]
486    fn test_experimental_variogram() {
487        let coords = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
488        let values = array![1.0, 2.0, 1.5, 2.5];
489
490        let result = experimental_variogram(&coords.view(), &values.view(), 5, None);
491        assert!(result.is_ok());
492
493        let (lags, gamma) = result.expect("computation failed");
494        assert!(!lags.is_empty());
495        assert_eq!(lags.len(), gamma.len());
496
497        // All gamma values should be non-negative
498        for &g in gamma.iter() {
499            assert!(g >= 0.0);
500        }
501    }
502
503    #[test]
504    fn test_fit_spherical_variogram() {
505        let lags = array![0.5, 1.0, 1.5, 2.0, 2.5];
506        let gamma = array![0.1, 0.3, 0.6, 0.85, 0.95];
507
508        let fitted = fit_variogram(&lags, &gamma, VariogramModel::Spherical);
509        assert!(fitted.is_ok());
510
511        let model = fitted.expect("fitting failed");
512        assert!(model.range > 0.0);
513        assert!(model.sill > 0.0);
514        assert!(model.nugget >= 0.0);
515    }
516
517    #[test]
518    fn test_fit_exponential_variogram() {
519        let lags = array![0.5, 1.0, 1.5, 2.0, 2.5];
520        let gamma = array![0.2, 0.4, 0.6, 0.75, 0.85];
521
522        let fitted = fit_variogram(&lags, &gamma, VariogramModel::Exponential);
523        assert!(fitted.is_ok());
524
525        let model = fitted.expect("fitting failed");
526        assert!(model.range > 0.0);
527        assert!(model.sill > 0.0);
528    }
529
530    #[test]
531    fn test_variogram_evaluate() {
532        let fitted = FittedVariogram {
533            model: VariogramModel::Spherical,
534            range: 2.0,
535            sill: 1.0,
536            nugget: 0.1,
537            extra_params: vec![],
538            r_squared: 0.95,
539        };
540
541        // At zero distance, should be close to nugget
542        let gamma_0 = fitted.evaluate(0.0);
543        assert_relative_eq!(gamma_0, 0.1, epsilon = 0.01);
544
545        // At range, should approach nugget + sill
546        let gamma_range = fitted.evaluate(2.0);
547        assert!(gamma_range >= 1.0);
548        assert!(gamma_range <= 1.2);
549
550        // Beyond range, should be nugget + sill
551        let gamma_beyond = fitted.evaluate(5.0);
552        assert_relative_eq!(gamma_beyond, 1.1, epsilon = 0.01);
553    }
554
555    #[test]
556    fn test_directional_variogram() {
557        let coords = array![
558            [0.0, 0.0],
559            [1.0, 0.0],
560            [2.0, 0.0],
561            [0.0, 1.0],
562            [1.0, 1.0],
563            [2.0, 1.0]
564        ];
565        let values = array![1.0, 1.5, 2.0, 1.2, 1.7, 2.2];
566
567        // East direction (0 radians)
568        let result = directional_variogram(
569            &coords.view(),
570            &values.view(),
571            0.0,
572            std::f64::consts::PI / 4.0, // 45 degree tolerance
573            5,
574        );
575
576        assert!(result.is_ok());
577        let (lags, gamma) = result.expect("computation failed");
578        assert!(!lags.is_empty());
579    }
580
581    #[test]
582    fn test_variogram_models() {
583        let models = vec![
584            VariogramModel::Spherical,
585            VariogramModel::Exponential,
586            VariogramModel::Gaussian,
587            VariogramModel::Linear,
588        ];
589
590        for model in models {
591            let fitted = FittedVariogram {
592                model,
593                range: 1.0,
594                sill: 1.0,
595                nugget: 0.0,
596                extra_params: vec![],
597                r_squared: 0.0,
598            };
599
600            // All models should be monotonically increasing
601            let gamma_1 = fitted.evaluate(0.5);
602            let gamma_2 = fitted.evaluate(1.0);
603            assert!(gamma_2 >= gamma_1);
604        }
605    }
606}
scirs2_spatial/variogram.rs

scirs2_spatial/
variogram.rs
