gistools/readers/grib2/sections/_3/templates.rs
1// TEMPLATE INFO: https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-1.shtml
2use crate::{
3 parsers::Reader,
4 readers::{Grib2Table3_1, Grib2Table3_2, Grib2Table3_3, Grib2Table3_4, Grib2Table3_5},
5};
6use alloc::vec;
7use s2json::{VectorMultiPoint, VectorPoint};
8
9/// # Grid Units
10///
11/// ## Links
12/// - [Docs](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect3.shtml)
13#[derive(Debug, Clone, Copy, PartialEq, Eq)]
14pub enum Grib2GridUnits {
15 /// degrees
16 Degrees,
17 /// meters
18 Meters,
19}
20
21/// Returns a template generator for the given template number
22/// All templates are listed [here](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-1.shtml)
23#[derive(Debug, Clone, PartialEq)]
24pub enum Grib2Template3 {
25 /// Latitude/Longitude (or equidistant cylindrical, or Plate Carree)
26 EquatorialTemplate(EquatorialTemplate),
27 /// Polar Stereographic Projection (Can be North or South)
28 PolarTemplate(PolarTemplate),
29}
30impl Grib2Template3 {
31 /// Create a new instance of Grib2Template3
32 ///
33 /// - `template`: template number parse block
34 /// - `section`: byte block
35 ///
36 /// ## Returns
37 /// Template generator
38 pub fn new<T: Reader>(template: Grib2Table3_1, section: &T) -> Self {
39 // TODO: Addd all Grib2Table3_1 options and set correct transform
40 match template {
41 Grib2Table3_1::LatitudeLongitude => {
42 Grib2Template3::EquatorialTemplate(EquatorialTemplate::new(section))
43 }
44 Grib2Table3_1::PolarStereographicProjection => {
45 Grib2Template3::PolarTemplate(PolarTemplate::new(section))
46 }
47 _ => panic!("Template 3.{template} not defined"),
48 }
49 }
50 /// Convert this section into grid data
51 pub fn build_grid(&mut self) -> VectorMultiPoint {
52 match self {
53 Grib2Template3::EquatorialTemplate(template) => template.build_grid(),
54 Grib2Template3::PolarTemplate(template) => template.build_grid(),
55 }
56 }
57}
58
59/// # GRIB2 - GRID DEFINITION TEMPLATE 3.0
60///
61/// ## Latitude/Longitude (or equidistant cylindrical, or Plate Carree)
62///
63/// ## Links
64///
65/// - [Read more...](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-0.shtml)
66///
67/// ## Notes
68///
69/// - Basic angle of the initial production domain and subdivisions of this basic angle are provided
70/// to manage cases where the recommended unit of 10-6 degrees is not applicable to describe the
71/// extreme longitudes and latitudes, and direction increments. For these last six descriptors, the
72/// unit is equal to the ratio of the basic angle and the subdivisions number. For ordinary cases,
73/// zero and missing values should be coded, equivalent to respective values of 1 and 106 (10-6
74/// degrees unit).
75///
76/// - For data on a quasi-regular grid, in which all the rows or columns do not necessarily have the
77/// same number of grid points either Ni (octets 31-34) of Nj (octets 35-38) and the corresponding Di
78/// (octets 64-67) or Dj (octets 68-71) shall be coded with all bits set to 1 (missing). The actual
79/// number of points along each parallel or meridian shall be coded in the octets immediately following
80/// the grid definition template (octets [xx+1]-nn), as described in the description of the grid
81/// definition section.
82///
83/// - A quasi-regular grid is only defined for appropriate grid scanning modes. Either rows or columns,
84/// but not both simultaneously, may have variable numbers of points or variable spacing. The first
85/// point in each row (column) shall be positioned at the meridian (parallel) indicted by octets 47-54.
86/// The grid points shall be evenly spaced in latitude (longitude).
87///
88/// A scale value of radius of spherical Earth, or major axis of oblate spheroid Earth is delivered
89/// from applying appropriate scale factor to the value expressed in meters.
90///
91/// - It is recommended to use unsigned direction increments.
92///
93/// - In most cases, multiplying Ni (octets 31-34) by Nj (octets 35-38) yields the total number of
94/// points in the grid. However, this may not be true if bit 8 of the scanning mode flags (octet 72)
95/// is set to 1.
96#[derive(Debug, Clone, PartialEq)]
97pub struct EquatorialTemplate {
98 /// Shape of Earth [Table 3.2](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-2.shtml)
99 pub shape: Grib2Table3_2,
100 /// Scale Factor of radius of spherical Earth
101 pub radius_scale_factor: u8,
102 /// Scale value of radius of spherical Earth
103 pub radius_scale_value: u32,
104 /// Scale factor of major axis of oblate spheroid Earth
105 pub major_axis_scale_factor: u8,
106 /// Scale value of major axis of oblate spheroid Earth
107 pub major_axis_scale_value: u32,
108 /// Scale factor of minor axis of oblate spheroid Earth
109 pub minor_axis_scale_factor: u8,
110 /// Scale value of minor axis of oblate spheroid Earth
111 pub minor_axis_scale_value: u32,
112 /// Number of points along a parallel (W-E)
113 pub nx: u32,
114 /// Number of points along a meridian (N-S)
115 pub ny: u32,
116 /// Basic angle of the initial production domain
117 pub basic_angle: f64,
118 /// Subdivisions of basic angle used to define extreme longitudes and latitudes, and direction increments
119 pub subdivisions: f64,
120 /// Latitude of first grid point
121 pub lat1: f64,
122 /// Longitude of first grid point
123 pub lon1: f64,
124 /// Resolution and component flags [Table 3.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-3.shtml)
125 pub resolution: Grib2Table3_3,
126 /// Latitude of last grid point
127 pub lat2: f64,
128 /// Longitude of last grid point
129 pub lon2: f64,
130 /// i direction increment
131 pub dx: f64,
132 /// j direction increment
133 pub dy: f64,
134 /// Scanning mode [Table 3.4](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-4.shtml)
135 pub scan_mode: Grib2Table3_4,
136 /// Grid Units
137 pub grid_units: Grib2GridUnits,
138}
139impl EquatorialTemplate {
140 /// Create a new instance of EquatorialTemplate
141 ///
142 /// ## Parameters
143 /// - `section`: byte block for template 3.0
144 ///
145 /// ## Returns
146 /// The parsed template
147 pub fn new<T: Reader>(section: &T) -> Self {
148 let shape = section.uint8(Some(14));
149 let basic_angle = section.uint32_be(Some(38)) as f64;
150 let subdivisions = section.uint32_be(Some(42)) as f64;
151 let lat1 = section.int32_be(Some(46));
152 let lat2 = section.int32_be(Some(55));
153 // build resolution values
154 let resolution_code = section.uint8(Some(54));
155 // build scan_mode values
156 let scan_mode_code = section.uint8(Some(71));
157
158 let ratio = if basic_angle == 0. { 1.0e-6 } else { basic_angle / subdivisions };
159
160 Self {
161 shape: shape.into(),
162 radius_scale_factor: section.uint8(Some(15)),
163 radius_scale_value: section.uint32_be(Some(16)),
164 major_axis_scale_factor: section.uint8(Some(20)),
165 major_axis_scale_value: section.uint32_be(Some(21)),
166 minor_axis_scale_factor: section.uint8(Some(25)),
167 minor_axis_scale_value: section.uint32_be(Some(26)),
168 nx: section.uint32_be(Some(30)),
169 ny: section.uint32_be(Some(34)),
170 basic_angle,
171 subdivisions,
172 lat1: (if lat1 < 0 { -lat1 ^ 0x80000000u32 as i32 } else { lat1 }) as f64 * ratio,
173 lon1: section.int32_be(Some(50)) as f64 * ratio,
174 resolution: resolution_code.into(),
175 lat2: (if lat2 < 0 { -lat2 ^ 0x80000000u32 as i32 } else { lat2 }) as f64 * ratio,
176 lon2: section.int32_be(Some(59)) as f64 * ratio,
177 dx: section.int32_be(Some(63)) as f64 * ratio,
178 dy: section.int32_be(Some(67)) as f64 * ratio,
179 scan_mode: scan_mode_code.into(),
180 grid_units: Grib2GridUnits::Degrees,
181 }
182 }
183 /// Convert this section into grid data
184 pub fn build_grid(&mut self) -> VectorMultiPoint {
185 // for now let's just follow the most basic scan mode
186 let Self { lat1, lat2, lon1, lon2, nx, ny, .. } = *self;
187 // Step sizes for interpolation
188 let lon_step = (lon2 - lon1) / (nx as f64 - 1.);
189 let lat_step = (lat2 - lat1) / (ny as f64 - 1.);
190
191 let mut res = vec![];
192
193 for y in 0..ny {
194 let y = y as f64;
195 for x in 0..nx {
196 let x = x as f64;
197 // Interpolate longitude and latitude
198 let lon = lon1 + x * lon_step;
199 let lat = lat1 + y * lat_step;
200 // create point and apply transform if provided (this grid is already in the correct projection)
201 res.push(VectorPoint::new_xy(lon, lat, None));
202 }
203 }
204
205 res
206 }
207}
208
209/// # GRIB2 - GRID DEFINITION TEMPLATE 3.20
210///
211/// ## Polar Stereographic Projection (Can be North or South)
212///
213/// ## Links
214/// - [Read more...](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-20.shtml)
215///
216/// ## Notes
217/// - The orientation of the grid is given by the longitude of the meridian along which the
218/// y-axis increases, LoV.
219/// - The projection is defined by the latitude at which Dx and Dy are specified, LaD.
220/// - Grid lengths Dx and Dy are in meters at the latitude LaD.
221/// - Bit 3 of the resolution and component flags should be set to 1 to indicate that Dx and Dy
222/// are given in meters.
223#[derive(Debug, Clone, PartialEq)]
224pub struct PolarTemplate {
225 /// Shape of Earth [Table 3.2](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-2.shtml)
226 pub shape: Grib2Table3_2,
227 /// Scale Factor of radius of spherical Earth
228 pub radius_scale_factor: u8,
229 /// Scale value of radius of spherical Earth
230 pub radius_scale_value: u32,
231 /// Scale factor of major axis of oblate spheroid Earth
232 pub major_axis_scale_factor: u8,
233 /// Scale value of major axis of oblate spheroid Earth
234 pub major_axis_scale_value: u32,
235 /// Scale factor of minor axis of oblate spheroid Earth
236 pub minor_axis_scale_factor: u8,
237 /// Scale value of minor axis of oblate spheroid Earth
238 pub minor_axis_scale_value: u32,
239 /// Number of points along the x-axis
240 pub nx: u32,
241 /// Number of points along the y-axis
242 pub ny: u32,
243 /// Latitude of first grid point
244 pub lat1: f64,
245 /// Longitude of first grid point
246 pub lon1: f64,
247 /// Latitude where Dx and Dy are specified
248 pub lat_d: f64,
249 /// Orientation of the grid (LoV)
250 pub lon_v: f64,
251 /// Resolution and component flags [Table 3.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-3.shtml)
252 pub resolution: Grib2Table3_3,
253 /// x-direction grid length (meters at LaD)
254 pub dx: f64,
255 /// y-direction grid length (meters at LaD)
256 pub dy: f64,
257 /// Projection center flag [Table 3.5](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-5.shtml)
258 pub proj_center: Grib2Table3_5,
259 /// Scanning mode [Table 3.4](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-4.shtml)
260 pub scan_mode: Grib2Table3_4,
261 /// Grid Units
262 pub grid_units: Grib2GridUnits,
263}
264impl PolarTemplate {
265 /// Create a new instance of PolarTemplate
266 ///
267 /// ## Parameters
268 /// - `section`: byte block for template 3.20
269 ///
270 /// ## Returns
271 /// The parsed template
272 pub fn new<T: Reader>(section: &T) -> Self {
273 let shape = section.uint8(Some(14));
274 let lat1 = section.int32_be(Some(38)) as f64;
275 let lon1 = section.int32_be(Some(42)) as f64;
276 let lat_d = section.int32_be(Some(47)) as f64;
277 let lon_v = section.int32_be(Some(51)) as f64;
278 let dx = section.int32_be(Some(55)) as f64;
279 let dy = section.int32_be(Some(59)) as f64;
280 let proj_center = section.uint8(Some(63));
281 // build resolution values
282 let resolution_code = section.uint8(Some(54));
283 // build scan mode
284 let scan_mode_code = section.uint8(Some(64));
285
286 Self {
287 shape: shape.into(),
288 radius_scale_factor: section.uint8(Some(15)),
289 radius_scale_value: section.uint32_be(Some(16)),
290 major_axis_scale_factor: section.uint8(Some(20)),
291 major_axis_scale_value: section.uint32_be(Some(21)),
292 minor_axis_scale_factor: section.uint8(Some(25)),
293 minor_axis_scale_value: section.uint32_be(Some(26)),
294 nx: section.uint32_be(Some(30)),
295 ny: section.uint32_be(Some(34)),
296 lat1,
297 lon1,
298 lat_d,
299 lon_v,
300 resolution: resolution_code.into(),
301 dx,
302 dy,
303 proj_center: proj_center.into(),
304 scan_mode: scan_mode_code.into(),
305 grid_units: Grib2GridUnits::Meters,
306 }
307 }
308 /// Convert this section into grid data
309 pub fn build_grid(&mut self) -> VectorMultiPoint {
310 // for now let's just follow the most basic scan mode
311 let Self { lat1, dx, lon1, dy, nx, ny, .. } = *self;
312 let mut res = vec![];
313
314 for y in 0..ny {
315 let y = y as f64;
316 for x in 0..nx {
317 let x = x as f64;
318 // Interpolate longitude and latitude
319 let lon = lon1 + x * dx;
320 let lat = lat1 + y * dy;
321 // create point
322 let point = VectorPoint::new_xy(lon, lat, None);
323 // apply transform if provided
324 // if (transformer !== undefined) point = transformer.forward(point);
325 res.push(point);
326 }
327 }
328
329 res
330 }
331}