nexrad_process/filter/
clutter.rs

1use crate::result::{Error, Result};
2use crate::SweepProcessor;
3use nexrad_model::data::{GateStatus, SweepField};
4
5/// Correlation coefficient-based clutter removal.
6///
7/// Masks gates in a target field (e.g., reflectivity) where the corresponding
8/// gate in a correlation coefficient (CC) field falls below a threshold. This is
9/// a standard technique for removing non-meteorological echoes such as ground
10/// clutter, biological targets, and debris.
11///
12/// Typical CC thresholds range from 0.85 to 0.95. Pure precipitation has
13/// CC values near 1.0, while clutter and biological targets tend to have
14/// CC well below 0.9.
15///
16/// # Example
17///
18/// ```ignore
19/// use nexrad_process::filter::CorrelationCoefficientFilter;
20///
21/// let cc_field = SweepField::from_radials(radials, Product::CorrelationCoefficient).unwrap();
22/// let filter = CorrelationCoefficientFilter::new(0.90, cc_field);
23/// let cleaned = filter.process(&reflectivity_field)?;
24/// ```
25pub struct CorrelationCoefficientFilter {
26    /// Gates with CC below this value are masked in the target field.
27    threshold: f32,
28    /// The correlation coefficient field to compare against.
29    cc_field: SweepField,
30}
31
32impl CorrelationCoefficientFilter {
33    /// Create a new CC filter with the given threshold and CC field.
34    ///
35    /// # Errors
36    ///
37    /// Returns an error if the threshold is not in (0, 1].
38    pub fn new(threshold: f32, cc_field: SweepField) -> Result<Self> {
39        if threshold <= 0.0 || threshold > 1.0 {
40            return Err(Error::InvalidParameter(
41                "CC threshold must be in (0, 1]".to_string(),
42            ));
43        }
44        Ok(Self {
45            threshold,
46            cc_field,
47        })
48    }
49}
50
51impl SweepProcessor for CorrelationCoefficientFilter {
52    fn name(&self) -> &str {
53        "CorrelationCoefficientFilter"
54    }
55
56    fn process(&self, input: &SweepField) -> Result<SweepField> {
57        if input.azimuth_count() != self.cc_field.azimuth_count() {
58            return Err(Error::InvalidGeometry(format!(
59                "CC field azimuth count ({}) does not match input field ({})",
60                self.cc_field.azimuth_count(),
61                input.azimuth_count(),
62            )));
63        }
64
65        // When gate counts differ (e.g., REF at 460 km vs CC at 300 km), apply
66        // the filter over the shared range and leave remaining gates unchanged.
67        let shared_gates = input.gate_count().min(self.cc_field.gate_count());
68
69        let mut output = input.clone();
70
71        for az_idx in 0..input.azimuth_count() {
72            for gate_idx in 0..shared_gates {
73                let (_, input_status) = input.get(az_idx, gate_idx);
74                if input_status != GateStatus::Valid {
75                    continue;
76                }
77
78                let (cc_val, cc_status) = self.cc_field.get(az_idx, gate_idx);
79
80                // If CC is invalid or below threshold, mask the target gate
81                if cc_status != GateStatus::Valid || cc_val < self.threshold {
82                    output.set(az_idx, gate_idx, 0.0, GateStatus::NoData);
83                }
84            }
85        }
86
87        Ok(output)
88    }
89}
90
91#[cfg(test)]
92mod tests {
93    use super::*;
94
95    fn make_fields(gate_count: usize) -> (SweepField, SweepField) {
96        let azimuths = vec![0.0, 1.0, 2.0];
97
98        let mut target = SweepField::new_empty(
99            "Reflectivity",
100            "dBZ",
101            0.5,
102            azimuths.clone(),
103            1.0,
104            2.0,
105            0.25,
106            gate_count,
107        );
108
109        let mut cc = SweepField::new_empty("CC", "", 0.5, azimuths, 1.0, 2.0, 0.25, gate_count);
110
111        // Fill target with valid data
112        for az in 0..3 {
113            for gate in 0..gate_count {
114                target.set(az, gate, 30.0, GateStatus::Valid);
115                cc.set(az, gate, 0.98, GateStatus::Valid);
116            }
117        }
118
119        (target, cc)
120    }
121
122    #[test]
123    fn test_cc_filter_preserves_high_cc() {
124        let (target, cc) = make_fields(5);
125        let filter = CorrelationCoefficientFilter::new(0.90, cc).unwrap();
126        let result = filter.process(&target).unwrap();
127
128        for az in 0..3 {
129            for gate in 0..5 {
130                let (val, status) = result.get(az, gate);
131                assert_eq!(val, 30.0);
132                assert_eq!(status, GateStatus::Valid);
133            }
134        }
135    }
136
137    #[test]
138    fn test_cc_filter_masks_low_cc() {
139        let (target, mut cc) = make_fields(5);
140        // Set low CC for some gates
141        cc.set(1, 2, 0.5, GateStatus::Valid);
142        cc.set(1, 3, 0.8, GateStatus::Valid);
143
144        let filter = CorrelationCoefficientFilter::new(0.90, cc).unwrap();
145        let result = filter.process(&target).unwrap();
146
147        // Low CC gates should be masked
148        let (_, status) = result.get(1, 2);
149        assert_eq!(status, GateStatus::NoData);
150
151        let (_, status) = result.get(1, 3);
152        assert_eq!(status, GateStatus::NoData);
153
154        // High CC gate should be preserved
155        let (val, status) = result.get(1, 1);
156        assert_eq!(val, 30.0);
157        assert_eq!(status, GateStatus::Valid);
158    }
159
160    #[test]
161    fn test_cc_filter_masks_invalid_cc() {
162        let (target, mut cc) = make_fields(5);
163        cc.set(0, 0, 0.0, GateStatus::NoData);
164
165        let filter = CorrelationCoefficientFilter::new(0.90, cc).unwrap();
166        let result = filter.process(&target).unwrap();
167
168        let (_, status) = result.get(0, 0);
169        assert_eq!(status, GateStatus::NoData);
170    }
171
172    #[test]
173    fn test_cc_filter_invalid_threshold() {
174        let (_, cc) = make_fields(5);
175        assert!(CorrelationCoefficientFilter::new(0.0, cc).is_err());
176
177        let (_, cc) = make_fields(5);
178        assert!(CorrelationCoefficientFilter::new(1.5, cc).is_err());
179    }
180
181    #[test]
182    fn test_cc_filter_different_gate_counts_ok() {
183        // REF with more gates than CC (common: REF 460 km vs CC 300 km)
184        let (target, _) = make_fields(10);
185        let (_, cc_shorter) = make_fields(5);
186
187        let filter = CorrelationCoefficientFilter::new(0.90, cc_shorter).unwrap();
188        let result = filter.process(&target).unwrap();
189
190        // Gates in the shared range should be preserved (high CC)
191        let (val, status) = result.get(0, 2);
192        assert_eq!(val, 30.0);
193        assert_eq!(status, GateStatus::Valid);
194
195        // Gates beyond CC range should be left unchanged
196        let (val, status) = result.get(0, 7);
197        assert_eq!(val, 30.0);
198        assert_eq!(status, GateStatus::Valid);
199    }
200
201    #[test]
202    fn test_cc_filter_azimuth_mismatch_error() {
203        // Different azimuth counts should still error
204        let azimuths_3 = vec![0.0, 1.0, 2.0];
205        let azimuths_4 = vec![0.0, 1.0, 2.0, 3.0];
206
207        let mut target =
208            SweepField::new_empty("Reflectivity", "dBZ", 0.5, azimuths_3, 1.0, 2.0, 0.25, 5);
209        let mut cc = SweepField::new_empty("CC", "", 0.5, azimuths_4, 1.0, 2.0, 0.25, 5);
210
211        for az in 0..3 {
212            for gate in 0..5 {
213                target.set(az, gate, 30.0, GateStatus::Valid);
214            }
215        }
216        for az in 0..4 {
217            for gate in 0..5 {
218                cc.set(az, gate, 0.98, GateStatus::Valid);
219            }
220        }
221
222        let filter = CorrelationCoefficientFilter::new(0.90, cc).unwrap();
223        assert!(filter.process(&target).is_err());
224    }
225}
nexrad_process/filter/clutter.rs

nexrad_process/filter/
clutter.rs
