nexrad_process/derived/vertical.rs
1use crate::result::{Error, Result};
2use nexrad_model::data::{GateStatus, SweepField, VerticalField};
3
4/// Effective earth radius in km using the standard 4/3 refraction model.
5const EFFECTIVE_EARTH_RADIUS_KM: f64 = 6371.0 * 4.0 / 3.0;
6
7/// Vertical cross-section (RHI-style) — assembles a range-height display from
8/// PPI scan data at a fixed azimuth.
9///
10/// For each cell in the output grid, the beam-height equation (standard 4/3
11/// earth-radius model) is used to map range and altitude back to each
12/// elevation tilt's polar coordinates. The maximum valid value across all
13/// tilts is retained, producing a pseudo-RHI from PPI scan data.
14///
15/// Unlike [`CompositeReflectivity`](super::CompositeReflectivity), this type does not
16/// implement [`ScanDerivedProduct`](crate::ScanDerivedProduct) because its output is a
17/// [`VerticalField`] (range vs. altitude) rather than a
18/// [`CartesianField`](nexrad_model::data::CartesianField). Its
19/// construction parameters (azimuth, range, altitude, dimensions) fully define the
20/// output grid, so no coordinate system or geographic extent is needed at compute time.
21///
22/// # Example
23///
24/// ```ignore
25/// use nexrad_process::derived::VerticalCrossSection;
26/// use nexrad_model::data::{SweepField, Product};
27///
28/// // Extract reflectivity fields from all sweeps
29/// let ref_fields: Vec<SweepField> = scan.sweeps().iter()
30/// .filter_map(|s| SweepField::from_radials(s.radials(), Product::Reflectivity))
31/// .collect();
32///
33/// // Create a cross-section at 200° azimuth, 0-230 km range, 0-18 km altitude
34/// let vcs = VerticalCrossSection::new(200.0, 230.0, 18000.0, 600, 300)?;
35/// let vertical_field = vcs.compute(&ref_fields)?;
36///
37/// // Render with nexrad-render
38/// use nexrad_render::{render_vertical, default_color_scale, RenderOptions};
39/// let scale = default_color_scale(Product::Reflectivity);
40/// let result = render_vertical(&vertical_field, &scale, &RenderOptions::new(1200, 600))?;
41/// result.save("vertical_cross_section.png")?;
42/// ```
43pub struct VerticalCrossSection {
44 /// Target azimuth in degrees (0-360).
45 azimuth_degrees: f32,
46 /// Maximum range in km for the horizontal axis.
47 max_range_km: f64,
48 /// Maximum altitude in meters for the vertical axis.
49 max_altitude_m: f64,
50 /// Output grid width (range bins).
51 width: usize,
52 /// Output grid height (altitude bins).
53 height: usize,
54}
55
56impl VerticalCrossSection {
57 /// Create a new vertical cross-section builder.
58 ///
59 /// # Parameters
60 ///
61 /// - `azimuth_degrees` — target azimuth (0-360 degrees)
62 /// - `max_range_km` — horizontal extent in km
63 /// - `max_altitude_m` — vertical extent in meters
64 /// - `width` — number of horizontal bins (columns)
65 /// - `height` — number of vertical bins (rows)
66 ///
67 /// # Errors
68 ///
69 /// Returns an error if any parameter is non-positive.
70 pub fn new(
71 azimuth_degrees: f32,
72 max_range_km: f64,
73 max_altitude_m: f64,
74 width: usize,
75 height: usize,
76 ) -> Result<Self> {
77 if max_range_km <= 0.0 {
78 return Err(Error::InvalidParameter(
79 "max_range_km must be positive".to_string(),
80 ));
81 }
82 if max_altitude_m <= 0.0 {
83 return Err(Error::InvalidParameter(
84 "max_altitude_m must be positive".to_string(),
85 ));
86 }
87 if width == 0 || height == 0 {
88 return Err(Error::InvalidParameter(
89 "output dimensions must be > 0".to_string(),
90 ));
91 }
92 Ok(Self {
93 azimuth_degrees,
94 max_range_km,
95 max_altitude_m,
96 width,
97 height,
98 })
99 }
100
101 /// Compute the vertical cross-section from sweep fields at multiple elevations.
102 ///
103 /// Each input field contributes data at the altitude determined by its elevation
104 /// angle and the beam-height equation. Where multiple tilts overlap the same
105 /// output cell, the maximum valid value is retained. The label and unit of the
106 /// output field are taken from the first input field.
107 ///
108 /// # Errors
109 ///
110 /// Returns an error if no fields are provided.
111 pub fn compute(&self, fields: &[SweepField]) -> Result<VerticalField> {
112 if fields.is_empty() {
113 return Err(Error::MissingData("no sweep fields provided".to_string()));
114 }
115
116 let label = format!("{} RHI", fields[0].label());
117 let unit = fields[0].unit().to_string();
118
119 let mut output = VerticalField::new(
120 label,
121 unit,
122 (0.0, self.max_range_km),
123 (0.0, self.max_altitude_m),
124 self.width,
125 self.height,
126 );
127
128 let re = EFFECTIVE_EARTH_RADIUS_KM;
129
130 for field in fields {
131 let elev_rad = (field.elevation_degrees() as f64).to_radians();
132
133 for col in 0..self.width {
134 let range_km = (col as f64 + 0.5) / self.width as f64 * self.max_range_km;
135
136 // Beam height using 4/3 earth radius model
137 let altitude_km = range_km * elev_rad.sin() + (range_km * range_km) / (2.0 * re);
138 let altitude_m = altitude_km * 1000.0;
139
140 if altitude_m < 0.0 || altitude_m > self.max_altitude_m {
141 continue;
142 }
143
144 // Row 0 = top = highest altitude
145 let row = ((self.max_altitude_m - altitude_m) / self.max_altitude_m
146 * self.height as f64) as usize;
147 if row >= self.height {
148 continue;
149 }
150
151 if let Some((val, status)) = field.value_at_polar(self.azimuth_degrees, range_km) {
152 if status == GateStatus::Valid {
153 let (existing_val, existing_status) = output.get(row, col);
154 if existing_status != GateStatus::Valid || val > existing_val {
155 output.set(row, col, val, GateStatus::Valid);
156 }
157 }
158 }
159 }
160 }
161
162 Ok(output)
163 }
164}
165
166#[cfg(test)]
167mod tests {
168 use super::*;
169
170 fn make_uniform_field(elevation: f32, value: f32) -> SweepField {
171 let mut azimuths = Vec::new();
172 for i in 0..360 {
173 azimuths.push(i as f32);
174 }
175
176 let gate_count = 100;
177 let mut field = SweepField::new_empty(
178 "Reflectivity",
179 "dBZ",
180 elevation,
181 azimuths,
182 1.0,
183 2.0,
184 0.25,
185 gate_count,
186 );
187
188 for az in 0..360 {
189 for gate in 0..gate_count {
190 field.set(az, gate, value, GateStatus::Valid);
191 }
192 }
193
194 field
195 }
196
197 #[test]
198 fn test_vertical_basic() {
199 let fields = vec![make_uniform_field(0.5, 30.0), make_uniform_field(3.5, 40.0)];
200
201 let vcs = VerticalCrossSection::new(180.0, 25.0, 5000.0, 50, 25).unwrap();
202 let result = vcs.compute(&fields).unwrap();
203
204 assert_eq!(result.width(), 50);
205 assert_eq!(result.height(), 25);
206
207 // Should have some valid data (low-elevation beam at close range)
208 let mut found_valid = false;
209 for row in 0..result.height() {
210 for col in 0..result.width() {
211 let (_, status) = result.get(row, col);
212 if status == GateStatus::Valid {
213 found_valid = true;
214 }
215 }
216 }
217 assert!(found_valid, "expected at least some valid data");
218 }
219
220 #[test]
221 fn test_vertical_takes_max() {
222 let fields = vec![make_uniform_field(0.5, 20.0), make_uniform_field(0.5, 40.0)];
223
224 let vcs = VerticalCrossSection::new(180.0, 25.0, 5000.0, 50, 25).unwrap();
225 let result = vcs.compute(&fields).unwrap();
226
227 // Overlapping beams should take the max value
228 for row in 0..result.height() {
229 for col in 0..result.width() {
230 let (val, status) = result.get(row, col);
231 if status == GateStatus::Valid {
232 assert_eq!(val, 40.0);
233 }
234 }
235 }
236 }
237
238 #[test]
239 fn test_vertical_empty_fields_error() {
240 let vcs = VerticalCrossSection::new(180.0, 25.0, 5000.0, 50, 25).unwrap();
241 assert!(vcs.compute(&[]).is_err());
242 }
243
244 #[test]
245 fn test_vertical_invalid_params() {
246 assert!(VerticalCrossSection::new(180.0, 0.0, 5000.0, 50, 25).is_err());
247 assert!(VerticalCrossSection::new(180.0, 25.0, 0.0, 50, 25).is_err());
248 assert!(VerticalCrossSection::new(180.0, 25.0, 5000.0, 0, 25).is_err());
249 assert!(VerticalCrossSection::new(180.0, 25.0, 5000.0, 50, 0).is_err());
250 }
251}