1use crate::error::{Error, Result};
4use crate::util::grib_i32;
5
6#[derive(Debug, Clone, PartialEq)]
13#[non_exhaustive]
14pub enum GridDefinition {
15 LatLon(LatLonGrid),
17 LambertConformal(LambertConformalGrid),
19 Unsupported(u16),
21}
22
23#[derive(Debug, Clone, PartialEq)]
25pub struct LatLonGrid {
26 pub ni: u32,
27 pub nj: u32,
28 pub lat_first: i32,
29 pub lon_first: i32,
30 pub lat_last: i32,
31 pub lon_last: i32,
32 pub di: u32,
33 pub dj: u32,
34 pub scanning_mode: u8,
35}
36
37#[derive(Debug, Clone, PartialEq)]
39pub struct LambertConformalGrid {
40 pub number_of_points: u32,
41 pub shape_of_earth: u8,
42 pub scale_factor_radius: u8,
43 pub scaled_value_radius: u32,
44 pub scale_factor_major_axis: u8,
45 pub scaled_value_major_axis: u32,
46 pub scale_factor_minor_axis: u8,
47 pub scaled_value_minor_axis: u32,
48 pub nx: u32,
49 pub ny: u32,
50 pub lat_first: i32,
51 pub lon_first: u32,
52 pub resolution_and_component_flags: u8,
53 pub lat_d: i32,
54 pub lon_v: u32,
55 pub dx: u32,
56 pub dy: u32,
57 pub projection_center_flag: u8,
58 pub scanning_mode: u8,
59 pub latin1: i32,
60 pub latin2: i32,
61 pub lat_southern_pole: i32,
62 pub lon_southern_pole: u32,
63}
64
65impl GridDefinition {
66 pub fn template_number(&self) -> u16 {
71 match self {
72 Self::LatLon(_) => 0,
73 Self::LambertConformal(_) => 30,
74 Self::Unsupported(template) => *template,
75 }
76 }
77
78 pub fn as_lat_lon(&self) -> Option<&LatLonGrid> {
79 match self {
80 Self::LatLon(grid) => Some(grid),
81 _ => None,
82 }
83 }
84
85 pub fn as_lambert_conformal(&self) -> Option<&LambertConformalGrid> {
86 match self {
87 Self::LambertConformal(grid) => Some(grid),
88 _ => None,
89 }
90 }
91
92 pub fn unsupported_template(&self) -> Option<u16> {
93 match self {
94 Self::Unsupported(template) => Some(*template),
95 _ => None,
96 }
97 }
98
99 pub fn shape(&self) -> (usize, usize) {
100 match self {
101 Self::LatLon(g) => (g.ni as usize, g.nj as usize),
102 Self::LambertConformal(g) => (g.nx as usize, g.ny as usize),
103 Self::Unsupported(_) => (0, 0),
104 }
105 }
106
107 pub fn ndarray_shape(&self) -> Vec<usize> {
108 let (ni, nj) = self.shape();
109 match self {
110 Self::LatLon(_) | Self::LambertConformal(_) if ni > 0 && nj > 0 => vec![nj, ni],
111 _ => Vec::new(),
112 }
113 }
114
115 pub fn num_points(&self) -> usize {
116 match self {
117 Self::LatLon(_) => {
118 let (ni, nj) = self.shape();
119 ni.saturating_mul(nj)
120 }
121 Self::LambertConformal(g) => g.number_of_points as usize,
122 Self::Unsupported(_) => 0,
123 }
124 }
125
126 pub fn validate_supported_scan_order(&self) -> Result<()> {
127 match self {
128 Self::LatLon(grid) => grid.validate_supported_scan_order(),
129 Self::LambertConformal(grid) => grid.validate_supported_scan_order(),
130 Self::Unsupported(template) => Err(Error::UnsupportedGridTemplate(*template)),
131 }
132 }
133
134 pub fn reorder_for_ndarray_in_place<T>(&self, values: &mut [T]) -> Result<()> {
135 match self {
136 Self::LatLon(grid) => grid.reorder_for_ndarray_in_place(values),
137 Self::LambertConformal(grid) => grid.reorder_for_ndarray_in_place(values),
138 Self::Unsupported(template) => Err(Error::UnsupportedGridTemplate(*template)),
139 }
140 }
141
142 pub fn parse(section_bytes: &[u8]) -> Result<Self> {
143 if section_bytes.len() < 14 {
144 return Err(Error::InvalidSection {
145 section: 3,
146 reason: format!("expected at least 14 bytes, got {}", section_bytes.len()),
147 });
148 }
149 if section_bytes[4] != 3 {
150 return Err(Error::InvalidSection {
151 section: section_bytes[4],
152 reason: "not a grid definition section".into(),
153 });
154 }
155
156 let template = u16::from_be_bytes(section_bytes[12..14].try_into().unwrap());
157 match template {
158 0 => parse_latlon(section_bytes),
159 30 => parse_lambert_conformal(section_bytes),
160 _ => Ok(Self::Unsupported(template)),
161 }
162 }
163}
164
165impl LatLonGrid {
166 pub fn longitudes(&self) -> Vec<f64> {
167 let step = self.di as f64 / 1_000_000.0;
168 let signed_step = if self.i_scans_positive() { step } else { -step };
169 let start = self.lon_first as f64 / 1_000_000.0;
170 (0..self.ni)
171 .map(|index| start + signed_step * index as f64)
172 .collect()
173 }
174
175 pub fn latitudes(&self) -> Vec<f64> {
176 let step = self.dj as f64 / 1_000_000.0;
177 let signed_step = if self.j_scans_positive() { step } else { -step };
178 let start = self.lat_first as f64 / 1_000_000.0;
179 (0..self.nj)
180 .map(|index| start + signed_step * index as f64)
181 .collect()
182 }
183
184 pub fn reorder_for_ndarray<T>(&self, mut values: Vec<T>) -> Result<Vec<T>> {
185 self.reorder_grib_scan_to_ndarray_in_place(&mut values)?;
186 Ok(values)
187 }
188
189 pub fn reorder_for_ndarray_in_place<T>(&self, values: &mut [T]) -> Result<()> {
190 self.reorder_grib_scan_to_ndarray_in_place(values)
191 }
192
193 pub fn reorder_grib_scan_to_ndarray<T>(&self, mut values: Vec<T>) -> Result<Vec<T>> {
194 self.reorder_grib_scan_to_ndarray_in_place(&mut values)?;
195 Ok(values)
196 }
197
198 pub fn reorder_grib_scan_to_ndarray_in_place<T>(&self, values: &mut [T]) -> Result<()> {
199 transform_supported_scan_order_in_place(
200 values,
201 self.ni as usize,
202 self.nj as usize,
203 self.scanning_mode,
204 )
205 }
206
207 pub fn reorder_ndarray_to_grib_scan<T>(&self, mut values: Vec<T>) -> Result<Vec<T>> {
208 self.reorder_ndarray_to_grib_scan_in_place(&mut values)?;
209 Ok(values)
210 }
211
212 pub fn reorder_ndarray_to_grib_scan_in_place<T>(&self, values: &mut [T]) -> Result<()> {
213 transform_supported_scan_order_in_place(
214 values,
215 self.ni as usize,
216 self.nj as usize,
217 self.scanning_mode,
218 )
219 }
220
221 pub fn validate_supported_scan_order(&self) -> Result<()> {
222 validate_supported_scan_order(self.scanning_mode)
223 }
224
225 fn i_scans_positive(&self) -> bool {
226 i_scans_positive(self.scanning_mode)
227 }
228
229 fn j_scans_positive(&self) -> bool {
230 j_scans_positive(self.scanning_mode)
231 }
232}
233
234impl LambertConformalGrid {
235 pub fn reorder_for_ndarray_in_place<T>(&self, values: &mut [T]) -> Result<()> {
236 transform_supported_scan_order_in_place(
237 values,
238 self.nx as usize,
239 self.ny as usize,
240 self.scanning_mode,
241 )
242 }
243
244 pub fn validate_supported_scan_order(&self) -> Result<()> {
245 validate_supported_scan_order(self.scanning_mode)
246 }
247}
248
249fn transform_supported_scan_order_in_place<T>(
250 values: &mut [T],
251 ni: usize,
252 nj: usize,
253 scanning_mode: u8,
254) -> Result<()> {
255 validate_supported_scan_order(scanning_mode)?;
256 if values.len() != ni * nj {
257 return Err(Error::DataLengthMismatch {
258 expected: ni * nj,
259 actual: values.len(),
260 });
261 }
262
263 if adjacent_rows_alternate_direction(scanning_mode) {
264 reverse_alternating_rows(values, ni, nj, i_scans_positive(scanning_mode));
265 }
266
267 Ok(())
268}
269
270fn validate_supported_scan_order(scanning_mode: u8) -> Result<()> {
271 if i_points_are_consecutive(scanning_mode) {
272 Ok(())
273 } else {
274 Err(Error::UnsupportedScanningMode(scanning_mode))
275 }
276}
277
278fn i_scans_positive(scanning_mode: u8) -> bool {
279 scanning_mode & 0b1000_0000 == 0
280}
281
282fn j_scans_positive(scanning_mode: u8) -> bool {
283 scanning_mode & 0b0100_0000 != 0
284}
285
286fn i_points_are_consecutive(scanning_mode: u8) -> bool {
287 scanning_mode & 0b0010_0000 == 0
288}
289
290fn adjacent_rows_alternate_direction(scanning_mode: u8) -> bool {
291 scanning_mode & 0b0001_0000 != 0
292}
293
294fn reverse_alternating_rows<T>(values: &mut [T], ni: usize, nj: usize, i_scans_positive: bool) {
295 for row in 0..nj {
296 let reverse = if i_scans_positive {
297 row % 2 == 1
298 } else {
299 row % 2 == 0
300 };
301 if reverse {
302 values[row * ni..(row + 1) * ni].reverse();
303 }
304 }
305}
306
307fn parse_latlon(data: &[u8]) -> Result<GridDefinition> {
308 if data.len() < 72 {
309 return Err(Error::InvalidSection {
310 section: 3,
311 reason: format!("template 3.0 requires 72 bytes, got {}", data.len()),
312 });
313 }
314
315 let ni = u32::from_be_bytes(data[30..34].try_into().unwrap());
316 let nj = u32::from_be_bytes(data[34..38].try_into().unwrap());
317 let lat_first = grib_i32(&data[46..50]).unwrap();
318 let lon_first = grib_i32(&data[50..54]).unwrap();
319 let lat_last = grib_i32(&data[55..59]).unwrap();
320 let lon_last = grib_i32(&data[59..63]).unwrap();
321 let di = u32::from_be_bytes(data[63..67].try_into().unwrap());
322 let dj = u32::from_be_bytes(data[67..71].try_into().unwrap());
323 let scanning_mode = data[71];
324
325 Ok(GridDefinition::LatLon(LatLonGrid {
326 ni,
327 nj,
328 lat_first,
329 lon_first,
330 lat_last,
331 lon_last,
332 di,
333 dj,
334 scanning_mode,
335 }))
336}
337
338fn parse_lambert_conformal(data: &[u8]) -> Result<GridDefinition> {
339 if data.len() < 81 {
340 return Err(Error::InvalidSection {
341 section: 3,
342 reason: format!("template 3.30 requires 81 bytes, got {}", data.len()),
343 });
344 }
345
346 Ok(GridDefinition::LambertConformal(LambertConformalGrid {
347 number_of_points: u32::from_be_bytes(data[6..10].try_into().unwrap()),
348 shape_of_earth: data[14],
349 scale_factor_radius: data[15],
350 scaled_value_radius: u32::from_be_bytes(data[16..20].try_into().unwrap()),
351 scale_factor_major_axis: data[20],
352 scaled_value_major_axis: u32::from_be_bytes(data[21..25].try_into().unwrap()),
353 scale_factor_minor_axis: data[25],
354 scaled_value_minor_axis: u32::from_be_bytes(data[26..30].try_into().unwrap()),
355 nx: u32::from_be_bytes(data[30..34].try_into().unwrap()),
356 ny: u32::from_be_bytes(data[34..38].try_into().unwrap()),
357 lat_first: grib_i32(&data[38..42]).unwrap(),
358 lon_first: u32::from_be_bytes(data[42..46].try_into().unwrap()),
359 resolution_and_component_flags: data[46],
360 lat_d: grib_i32(&data[47..51]).unwrap(),
361 lon_v: u32::from_be_bytes(data[51..55].try_into().unwrap()),
362 dx: u32::from_be_bytes(data[55..59].try_into().unwrap()),
363 dy: u32::from_be_bytes(data[59..63].try_into().unwrap()),
364 projection_center_flag: data[63],
365 scanning_mode: data[64],
366 latin1: grib_i32(&data[65..69]).unwrap(),
367 latin2: grib_i32(&data[69..73]).unwrap(),
368 lat_southern_pole: grib_i32(&data[73..77]).unwrap(),
369 lon_southern_pole: u32::from_be_bytes(data[77..81].try_into().unwrap()),
370 }))
371}
372
373#[cfg(test)]
374mod tests {
375 use super::{GridDefinition, LambertConformalGrid, LatLonGrid};
376 use crate::binary::encode_wmo_i32;
377
378 #[test]
379 fn reports_latlon_shape() {
380 let grid = GridDefinition::LatLon(LatLonGrid {
381 ni: 3,
382 nj: 2,
383 lat_first: 50_000_000,
384 lon_first: -120_000_000,
385 lat_last: 49_000_000,
386 lon_last: -118_000_000,
387 di: 1_000_000,
388 dj: 1_000_000,
389 scanning_mode: 0,
390 });
391
392 assert_eq!(grid.shape(), (3, 2));
393 assert_eq!(grid.ndarray_shape(), vec![2, 3]);
394 assert_eq!(grid.template_number(), 0);
395 assert!(grid.as_lat_lon().is_some());
396 assert!(grid.as_lambert_conformal().is_none());
397 assert_eq!(grid.unsupported_template(), None);
398 match grid {
399 GridDefinition::LatLon(ref latlon) => {
400 assert_eq!(latlon.longitudes(), vec![-120.0, -119.0, -118.0]);
401 assert_eq!(latlon.latitudes(), vec![50.0, 49.0]);
402 }
403 other => panic!("expected lat/lon grid, got {other:?}"),
404 }
405 }
406
407 #[test]
408 fn parses_lambert_conformal_template() {
409 let section = build_lambert_section();
410 let grid = GridDefinition::parse(§ion).unwrap();
411
412 assert_eq!(grid.shape(), (3, 2));
413 assert_eq!(grid.ndarray_shape(), vec![2, 3]);
414 assert_eq!(grid.num_points(), 6);
415 assert_eq!(grid.template_number(), 30);
416 assert!(grid.as_lat_lon().is_none());
417 assert!(grid.as_lambert_conformal().is_some());
418 assert_eq!(grid.unsupported_template(), None);
419 match grid {
420 GridDefinition::LambertConformal(lambert) => {
421 assert_eq!(lambert.number_of_points, 6);
422 assert_eq!(lambert.shape_of_earth, 1);
423 assert_eq!(lambert.scaled_value_radius, 6_371_200);
424 assert_eq!(lambert.nx, 3);
425 assert_eq!(lambert.ny, 2);
426 assert_eq!(lambert.lat_first, 12_190_000);
427 assert_eq!(lambert.lon_first, 226_541_000);
428 assert_eq!(lambert.resolution_and_component_flags, 0x08);
429 assert_eq!(lambert.lat_d, 25_000_000);
430 assert_eq!(lambert.lon_v, 265_000_000);
431 assert_eq!(lambert.dx, 2_539_703);
432 assert_eq!(lambert.dy, 2_539_703);
433 assert_eq!(lambert.projection_center_flag, 0);
434 assert_eq!(lambert.scanning_mode, 0);
435 assert_eq!(lambert.latin1, 25_000_000);
436 assert_eq!(lambert.latin2, 25_000_000);
437 assert_eq!(lambert.lat_southern_pole, -90_000_000);
438 assert_eq!(lambert.lon_southern_pole, 0);
439 }
440 other => panic!("expected Lambert conformal grid, got {other:?}"),
441 }
442 }
443
444 #[test]
445 fn reports_unsupported_template_helpers() {
446 let grid = GridDefinition::Unsupported(3_276);
447
448 assert_eq!(grid.template_number(), 3_276);
449 assert!(grid.as_lat_lon().is_none());
450 assert!(grid.as_lambert_conformal().is_none());
451 assert_eq!(grid.unsupported_template(), Some(3_276));
452 }
453
454 #[test]
455 fn normalizes_alternating_row_scan() {
456 let grid = LatLonGrid {
457 ni: 3,
458 nj: 2,
459 lat_first: 0,
460 lon_first: 0,
461 lat_last: 0,
462 lon_last: 0,
463 di: 1,
464 dj: 1,
465 scanning_mode: 0b0001_0000,
466 };
467
468 let ordered = grid
469 .reorder_for_ndarray(vec![1.0, 2.0, 3.0, 6.0, 5.0, 4.0])
470 .unwrap();
471 assert_eq!(ordered, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
472 }
473
474 #[test]
475 fn converts_ndarray_order_to_alternating_scan_order() {
476 let grid = LatLonGrid {
477 ni: 3,
478 nj: 2,
479 lat_first: 0,
480 lon_first: 0,
481 lat_last: 0,
482 lon_last: 0,
483 di: 1,
484 dj: 1,
485 scanning_mode: 0b0001_0000,
486 };
487
488 let scan_order = grid
489 .reorder_ndarray_to_grib_scan(vec![1, 2, 3, 4, 5, 6])
490 .unwrap();
491 assert_eq!(scan_order, vec![1, 2, 3, 6, 5, 4]);
492 }
493
494 #[test]
495 fn normalizes_lambert_alternating_row_scan() {
496 let grid = LambertConformalGrid {
497 number_of_points: 6,
498 shape_of_earth: 1,
499 scale_factor_radius: 0,
500 scaled_value_radius: 6_371_200,
501 scale_factor_major_axis: 0,
502 scaled_value_major_axis: 0,
503 scale_factor_minor_axis: 0,
504 scaled_value_minor_axis: 0,
505 nx: 3,
506 ny: 2,
507 lat_first: 0,
508 lon_first: 0,
509 resolution_and_component_flags: 0,
510 lat_d: 0,
511 lon_v: 0,
512 dx: 1,
513 dy: 1,
514 projection_center_flag: 0,
515 scanning_mode: 0b0001_0000,
516 latin1: 0,
517 latin2: 0,
518 lat_southern_pole: 0,
519 lon_southern_pole: 0,
520 };
521
522 let mut values = vec![1.0, 2.0, 3.0, 6.0, 5.0, 4.0];
523 grid.reorder_for_ndarray_in_place(&mut values).unwrap();
524 assert_eq!(values, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
525 }
526
527 #[test]
528 fn preserves_non_alternating_scan_modes_in_current_reader_order() {
529 for scanning_mode in [0b0000_0000, 0b1000_0000, 0b0100_0000, 0b1100_0000] {
530 let grid = LatLonGrid {
531 ni: 3,
532 nj: 2,
533 lat_first: 0,
534 lon_first: 0,
535 lat_last: 0,
536 lon_last: 0,
537 di: 1,
538 dj: 1,
539 scanning_mode,
540 };
541
542 let values = vec![1, 2, 3, 4, 5, 6];
543 assert_eq!(
544 grid.reorder_grib_scan_to_ndarray(values.clone()).unwrap(),
545 values
546 );
547 assert_eq!(
548 grid.reorder_ndarray_to_grib_scan(values.clone()).unwrap(),
549 values
550 );
551 }
552 }
553
554 #[test]
555 fn rejects_j_consecutive_scan_order() {
556 let grid = LatLonGrid {
557 ni: 3,
558 nj: 2,
559 lat_first: 0,
560 lon_first: 0,
561 lat_last: 0,
562 lon_last: 0,
563 di: 1,
564 dj: 1,
565 scanning_mode: 0b0010_0000,
566 };
567
568 let err = grid
569 .reorder_ndarray_to_grib_scan(vec![1, 2, 3, 4, 5, 6])
570 .unwrap_err();
571 assert!(matches!(
572 err,
573 crate::Error::UnsupportedScanningMode(0b0010_0000)
574 ));
575 }
576
577 fn build_lambert_section() -> Vec<u8> {
578 let mut section = vec![0u8; 81];
579 section[..4].copy_from_slice(&81u32.to_be_bytes());
580 section[4] = 3;
581 section[6..10].copy_from_slice(&6u32.to_be_bytes());
582 section[12..14].copy_from_slice(&30u16.to_be_bytes());
583 section[14] = 1;
584 section[16..20].copy_from_slice(&6_371_200u32.to_be_bytes());
585 section[30..34].copy_from_slice(&3u32.to_be_bytes());
586 section[34..38].copy_from_slice(&2u32.to_be_bytes());
587 section[38..42].copy_from_slice(&encode_wmo_i32(12_190_000).unwrap());
588 section[42..46].copy_from_slice(&226_541_000u32.to_be_bytes());
589 section[46] = 0x08;
590 section[47..51].copy_from_slice(&encode_wmo_i32(25_000_000).unwrap());
591 section[51..55].copy_from_slice(&265_000_000u32.to_be_bytes());
592 section[55..59].copy_from_slice(&2_539_703u32.to_be_bytes());
593 section[59..63].copy_from_slice(&2_539_703u32.to_be_bytes());
594 section[65..69].copy_from_slice(&encode_wmo_i32(25_000_000).unwrap());
595 section[69..73].copy_from_slice(&encode_wmo_i32(25_000_000).unwrap());
596 section[73..77].copy_from_slice(&encode_wmo_i32(-90_000_000).unwrap());
597 section
598 }
599}