diff --git a/geodesy/gsb/5458.gsa b/geodesy/gsb/5458.gsa new file mode 100644 index 0000000..40b1195 --- /dev/null +++ b/geodesy/gsb/5458.gsa @@ -0,0 +1,71 @@ +NUM_OREC 11 +NUM_SREC 11 +NUM_FILE 1 +GS_TYPE SECONDS +VERSION 2.0 +SYSTEM_F INTER +SYSTEM_T GRS80 +MAJOR_F 6378388.0 +MINOR_F 6356911.946127946 +MAJOR_T 6378137.0 +MINOR_T 6356752.314140356 + +SUB_NAME 5458 +PARENT NONE +CREATED 20231110 +UPDATED 20231110 +S_LAT 194400.0 +N_LAT 208800.0 +E_LONG -57600.0 +W_LONG -28800.0 +LAT_INC 3600.0 +LONG_INC 3600.0 +GS_COUNT 45 + +16 54 0 0 +15 54 0 0 +14 54 0 0 +13 54 0 0 +12 54 0 0 +11 54 0 0 +10 54 0 0 +09 54 0 0 +08 54 0 0 +16 55 0 0 +15 55 0 0 +14 55 0 0 +13 55 0 0 +12 55 0 0 +11 55 0 0 +10 55 0 0 +09 55 0 0 +08 55 0 0 +16 56 0 0 +15 56 0 0 +14 56 0 0 +13 56 0 0 +12 56 0 0 +11 56 0 0 +10 56 0 0 +09 56 0 0 +08 56 0 0 +16 57 0 0 +15 57 0 0 +14 57 0 0 +13 57 0 0 +12 57 0 0 +11 57 0 0 +10 57 0 0 +09 57 0 0 +08 57 0 0 +16 58 0 0 +15 58 0 0 +14 58 0 0 +13 58 0 0 +12 58 0 0 +11 58 0 0 +10 58 0 0 +09 58 0 0 +08 58 0 0 + +END diff --git a/geodesy/gsb/5458.gsb b/geodesy/gsb/5458.gsb new file mode 100644 index 0000000..c9c4f2b Binary files /dev/null and b/geodesy/gsb/5458.gsb differ diff --git a/geodesy/gsb/5458_with_subgrid.gsa b/geodesy/gsb/5458_with_subgrid.gsa new file mode 100644 index 0000000..ffc4d4b --- /dev/null +++ b/geodesy/gsb/5458_with_subgrid.gsa @@ -0,0 +1,99 @@ +NUM_OREC 11 +NUM_SREC 11 +NUM_FILE 2 +GS_TYPE SECONDS +VERSION 2.0 +SYSTEM_F INTER +SYSTEM_T GRS80 +MAJOR_F 6378388.0 +MINOR_F 6356911.946127946 +MAJOR_T 6378137.0 +MINOR_T 6356752.314140356 + +SUB_NAME 5458 +PARENT NONE +CREATED 20231110 +UPDATED 20231110 +S_LAT 194400.0 +N_LAT 208800.0 +E_LONG -57600.0 +W_LONG -28800.0 +LAT_INC 3600.0 +LONG_INC 3600.0 +GS_COUNT 45 + +16 54 0 0 +15 54 0 0 +14 54 0 0 +13 54 0 0 +12 54 0 0 +11 54 0 0 +10 54 0 0 +09 54 0 0 +08 54 0 0 +16 55 0 0 +15 55 0 0 +14 55 0 0 +13 55 0 0 +12 55 0 0 +11 55 0 0 +10 55 0 0 +09 55 0 0 +08 55 0 0 +16 56 0 0 +15 56 0 0 +14 56 0 0 +13 56 0 0 +12 56 0 0 +11 56 0 0 +10 56 0 0 +09 56 0 0 +08 56 0 0 +16 57 0 0 +15 57 0 0 +14 57 0 0 +13 57 0 0 +12 57 0 0 +11 57 0 0 +10 57 0 0 +09 57 0 0 +08 57 0 0 +16 58 0 0 +15 58 0 0 +14 58 0 0 +13 58 0 0 +12 58 0 0 +11 58 0 0 +10 58 0 0 +09 58 0 0 +08 58 0 0 + +SUB_NAME 5556 +PARENT 5458 +CREATED 20231110 +UPDATED 20231110 +S_LAT 198000.0 +N_LAT 201600.0 +E_LONG -50400.0 +W_LONG -43200.0 +LAT_INC 1800.0 +LONG_INC 1800.0 +GS_COUNT 15 + +14.0 56.0 0 0 +13.5 56.0 0 0 +13.0 56.0 0 0 +12.5 56.0 0 0 +12.0 56.0 0 0 +14.0 55.5 0 0 +13.5 55.5 0 0 +13.0 55.5 0 0 +12.5 55.5 0 0 +12.0 55.5 0 0 +14.0 55.0 0 0 +13.5 55.0 0 0 +13.0 55.0 0 0 +12.5 55.0 0 0 +12.0 55.0 0 0 + +END diff --git a/geodesy/gsb/5458_with_subgrid.gsb b/geodesy/gsb/5458_with_subgrid.gsb new file mode 100644 index 0000000..d36dd5c Binary files /dev/null and b/geodesy/gsb/5458_with_subgrid.gsb differ diff --git a/src/grid/mod.rs b/src/grid/mod.rs index 13408bf..40a24fd 100644 --- a/src/grid/mod.rs +++ b/src/grid/mod.rs @@ -25,10 +25,10 @@ pub trait Grid: Debug { /// geodetic grids in the Gravsoft format. #[derive(Debug, Default, Clone)] pub struct BaseGrid { - lat_0: f64, // Latitude of the first (typically northernmost) row of the grid - lat_1: f64, // Latitude of the last (typically southernmost) row of the grid - lon_0: f64, // Longitude of the first (typically westernmost) column of each row - lon_1: f64, // Longitude of the last (typically easternmost) column of each row + lat_n: f64, // Latitude of the first (typically northernmost) row of the grid + lat_s: f64, // Latitude of the last (typically southernmost) row of the grid + lon_w: f64, // Longitude of the first (typically westernmost) column of each row + lon_e: f64, // Longitude of the last (typically easternmost) column of each row dlat: f64, // Signed distance between two consecutive rows dlon: f64, // Signed distance between two consecutive columns rows: usize, @@ -47,8 +47,8 @@ impl Grid for BaseGrid { /// "On the border" qualifies as within. fn contains(&self, position: &Coor4D, margin: f64) -> bool { // We start by assuming that the last row (latitude) is the southernmost - let mut min = self.lat_1; - let mut max = self.lat_0; + let mut min = self.lat_s; + let mut max = self.lat_n; // If it's not, we swap if self.dlat > 0. { @@ -61,8 +61,8 @@ impl Grid for BaseGrid { } // The default assumption is the other way round for columns (longitudes) - min = self.lon_0; - max = self.lon_1; + min = self.lon_w; + max = self.lon_e; // If it's not, we swap if self.dlon < 0. { (min, max) = (max, min) @@ -89,49 +89,65 @@ impl Grid for BaseGrid { let grid = &self.grid; + // For now, we support top-to-bottom, left-to-right scan order only. + // This is the common case for most non-block grid formats, with + // NTv2 the odd man out. But since we normalize the NTv2 scan order + // during parsing, we just cruise along here + let dlat = self.dlat.abs(); + let dlon = self.dlon.abs(); + // The interpolation coordinate relative to the grid origin - let rlon = at[0] - self.lon_0; - let rlat = at[1] - self.lat_0; + let rlon = at[0] - self.lon_w; + let rlat = self.lat_n - at[1]; // The (row, column) of the lower left node of the grid cell containing // the interpolation coordinate - or, in the case of extrapolation: // the nearest cell inside the grid. - let row = (rlat / self.dlat).floor() as i64; - let col = (rlon / self.dlon).floor() as i64; + let row = (rlat / dlat).ceil() as i64; + let col = (rlon / dlon).floor() as i64; let col = col.clamp(0_i64, (self.cols - 2) as i64) as usize; let row = row.clamp(1_i64, (self.rows - 1) as i64) as usize; // Index of the first band element of each corner value #[rustfmt::skip] - let (ll, lr, ur, ul) = ( + let (ll, lr, ul, ur) = ( self.offset + self.bands * (self.cols * row + col ), self.offset + self.bands * (self.cols * row + col + 1), - self.offset + self.bands * (self.cols * (row - 1) + col + 1), self.offset + self.bands * (self.cols * (row - 1) + col ), + self.offset + self.bands * (self.cols * (row - 1) + col + 1), ); - let ll_lon = self.lon_0 + col as f64 * self.dlon; - let ll_lat = self.lat_0 + row as f64 * self.dlat; + let ll_lon = self.lon_w + col as f64 * dlon; + let ll_lat = self.lat_n - row as f64 * dlat; // Cell relative, cell unit coordinates in a right handed CS - let rlon = (at[0] - ll_lon) / self.dlon; - let rlat = (at[1] - ll_lat) / -self.dlat; + let rlon = (at[0] - ll_lon) / dlon; + let rlat = (at[1] - ll_lat) / dlat; - // Interpolate + // We cannot return more than 4 bands in a Coor4D, so we ignore + // any exceeding bands + let bands = self.bands.min(4); let mut left = Coor4D::origin(); - for i in 0..self.bands { - left[i] = (1. - rlat) * grid[ll + i] as f64 + rlat * grid[ul + i] as f64; + + // Interpolate (or extrapolate, if we're outside of the physical grid) + for i in 0..bands { + let lower = grid[ll + i] as f64; + let upper = grid[ul + i] as f64; + left[i] = (1. - rlat) * lower + rlat * upper; } let mut right = Coor4D::origin(); - for i in 0..self.bands { - right[i] = (1. - rlat) * grid[lr + i] as f64 + rlat * grid[ur + i] as f64; + for i in 0..bands { + let lower = grid[lr + i] as f64; + let upper = grid[ur + i] as f64; + right[i] = (1. - rlat) * lower + rlat * upper; } let mut result = Coor4D::origin(); - for i in 0..self.bands { + for i in 0..bands { result[i] = (1. - rlon) * left[i] + rlon * right[i]; } + Some(result) } } @@ -143,18 +159,18 @@ impl BaseGrid { offset: Option, ) -> Result { if header.len() < 7 { - return Err(Error::General("Incomplete grid")); + return Err(Error::General("Malformed header")); } - let lat_0 = header[0]; - let lat_1 = header[1]; - let lon_0 = header[2]; - let lon_1 = header[3]; - let dlat = header[4].copysign(lat_1 - lat_0); - let dlon = header[5].copysign(lon_1 - lon_0); + let lat_n = header[0]; + let lat_s = header[1]; + let lon_w = header[2]; + let lon_e = header[3]; + let dlat = header[4].copysign(lat_s - lat_n); + let dlon = header[5].copysign(lon_e - lon_w); let bands = header[6] as usize; - let rows = ((lat_1 - lat_0) / dlat + 1.5).floor() as usize; - let cols = ((lon_1 - lon_0) / dlon + 1.5).floor() as usize; + let rows = ((lat_s - lat_n) / dlat + 1.5).floor() as usize; + let cols = ((lon_e - lon_w) / dlon + 1.5).floor() as usize; let elements = rows * cols * bands; let offset = offset.unwrap_or(0); @@ -166,10 +182,10 @@ impl BaseGrid { } Ok(BaseGrid { - lat_0, - lat_1, - lon_0, - lon_1, + lat_n, + lat_s, + lon_w, + lon_e, dlat, dlon, rows, @@ -258,22 +274,22 @@ fn gravsoft_grid_reader(buf: &[u8]) -> Result<(Vec, Vec), Error> { return Err(Error::General("Incomplete Gravsoft header")); } - // The Gravsoft header has lat_1 before lat_0 + // The Gravsoft header has lat_s before lat_n header.swap(0, 1); // Count the number of bands - let lat_0 = header[0]; - let lat_1 = header[1]; - let lon_0 = header[2]; - let lon_1 = header[3]; + let lat_n = header[0]; + let lat_s = header[1]; + let lon_w = header[2]; + let lon_e = header[3]; // The Gravsoft header has inverted sign for dlat. We force // the two deltas to have signs compatible with the grid // organization - let dlat = header[4].copysign(lat_1 - lat_0); - let dlon = header[5].copysign(lon_1 - lon_0); - let rows = ((lat_1 - lat_0) / dlat + 1.5).floor() as usize; - let cols = ((lon_1 - lon_0) / dlon + 1.5).floor() as usize; + let dlat = header[4].copysign(lat_s - lat_n); + let dlon = header[5].copysign(lon_e - lon_w); + let rows = ((lat_s - lat_n) / dlat + 1.5).floor() as usize; + let cols = ((lon_e - lon_w) / dlon + 1.5).floor() as usize; let bands = grid.len() / (rows * cols); if (rows * cols * bands) > grid.len() || bands < 1 { return Err(Error::General("Incomplete Gravsoft grid")); @@ -304,7 +320,7 @@ fn gravsoft_grid_reader(buf: &[u8]) -> Result<(Vec, Vec), Error> { mod tests { use super::*; - // lat_0, lat_1, lon_0, lon_1, dlat, dlon + // lat_n, lat_s, lon_w, lon_e, dlat, dlon const HEADER: [f64; 6] = [58., 54., 8., 16., -1., 1.]; #[allow(dead_code)] @@ -326,33 +342,6 @@ mod tests { 54.08, 54.09, 54.10, 54.11, 54.12, 54.13, 54.14, 54.15, 54.16, ]; - // A geoid in inverse row order - #[rustfmt::skip] - const UPSIDE_DOWN_GEOID: [f32; 5*9] = [ - 54.08, 54.09, 54.10, 54.11, 54.12, 54.13, 54.14, 54.15, 54.16, - 55.08, 55.09, 55.10, 55.11, 55.12, 55.13, 55.14, 55.15, 55.16, - 56.08, 56.09, 56.10, 56.11, 56.12, 56.13, 56.14, 56.15, 56.16, - 57.08, 57.09, 57.10, 57.11, 57.12, 57.13, 57.14, 57.15, 57.16, - 58.08, 58.09, 58.10, 58.11, 58.12, 58.13, 58.14, 58.15, 58.16, - ]; - - #[rustfmt::skip] - const MIRRORED_GEOID: [f32; 5*9] = [ - 58.16, 58.15, 58.14, 58.13, 58.12, 58.11, 58.10, 58.09, 58.08, - 57.16, 57.15, 57.14, 57.13, 57.12, 57.11, 57.10, 57.09, 57.08, - 56.16, 56.15, 56.14, 56.13, 56.12, 56.11, 56.10, 56.09, 56.08, - 55.16, 55.15, 55.14, 55.13, 55.12, 55.11, 55.10, 55.09, 55.08, - 54.16, 54.15, 54.14, 54.13, 54.12, 54.11, 54.10, 54.09, 54.08, - ]; - - #[rustfmt::skip] - const MIRRORED_UPSIDE_DOWN_GEOID: [f32; 5*9] = [ - 54.16, 54.15, 54.14, 54.13, 54.12, 54.11, 54.10, 54.09, 54.08, - 55.16, 55.15, 55.14, 55.13, 55.12, 55.11, 55.10, 55.09, 55.08, - 56.16, 56.15, 56.14, 56.13, 56.12, 56.11, 56.10, 56.09, 56.08, - 57.16, 57.15, 57.14, 57.13, 57.12, 57.11, 57.10, 57.09, 57.08, - 58.16, 58.15, 58.14, 58.13, 58.12, 58.11, 58.10, 58.09, 58.08, - ]; #[test] fn grid_header() -> Result<(), Error> { // Create a datum correction grid (2 bands) @@ -371,10 +360,6 @@ mod tests { datum_header[4] = -datum_header[4]; let datum = BaseGrid::plain(&datum_header, Some(&datum_grid), None)?; - let c = Coor4D::geo(55.06, 12.03, 0., 0.); - let d = datum.at(&c, 1.0).unwrap(); - assert!(c.default_ellps_dist(&d.to_arcsec().to_radians()) < 1.0); - // Extrapolation let c = Coor4D::geo(100., 50., 0., 0.); // ...with output converted back to arcsec @@ -411,70 +396,6 @@ mod tests { let n = geoid.at(&c, 1.0).unwrap(); assert!((n[0] - (58.75 + 0.0825)).abs() < 0.0001); - - // Create an upside-down geoid grid (1 band) - let mut geoid_header = datum_header.clone(); - geoid_header.swap(0, 1); // lat_0=54, lat_1=58 - geoid_header[6] = 1.0; // 1 band - let geoid_grid = Vec::from(UPSIDE_DOWN_GEOID); - let geoid = BaseGrid::plain(&geoid_header, Some(&geoid_grid), None)?; - - let c = Coor4D::geo(58.75, 08.25, 0., 0.); - assert_eq!(geoid.contains(&c, 0.0), false); - assert_eq!(geoid.contains(&c, 1.0), true); - - let n = geoid.at(&c, 1.0).unwrap(); - assert!((n[0] - 58.83).abs() < 0.1); - - let c = Coor4D::geo(53.25, 8.0, 0., 0.); - assert_eq!(geoid.contains(&c, 0.0), false); - assert_eq!(geoid.contains(&c, 1.0), true); - - let n = geoid.at(&c, 1.0).unwrap(); - assert!((n[0] - (53.25 + 0.08)).abs() < 0.0001); - - // Create a mirrored geoid grid (1 band) - let mut geoid_header = datum_header.clone(); - geoid_header.swap(2, 3); // lon_0=16, lon_1=8 - geoid_header[5] = -geoid_header[5]; - geoid_header[6] = 1.0; // 1 band - let geoid_grid = Vec::from(MIRRORED_GEOID); - let geoid = BaseGrid::plain(&geoid_header, Some(&geoid_grid), None)?; - - let c = Coor4D::geo(58.75, 08.25, 0., 0.); - assert_eq!(geoid.contains(&c, 0.0), false); - assert_eq!(geoid.contains(&c, 1.0), true); - - let n = geoid.at(&c, 1.0).unwrap(); - assert!((n[0] - 58.83).abs() < 0.1); - - let c = Coor4D::geo(53.25, 8.0, 0., 0.); - assert_eq!(geoid.contains(&c, 0.0), false); - assert_eq!(geoid.contains(&c, 1.0), true); - - let n = geoid.at(&c, 1.0).unwrap(); - assert!((n[0] - (53.25 + 0.08)).abs() < 0.001); - - // Create a mirrored upside down geoid grid (1 band) - geoid_header.swap(0, 1); // lon_0=16, lon_1=8 - geoid_header[4] = -geoid_header[4]; - let geoid_grid = Vec::from(MIRRORED_UPSIDE_DOWN_GEOID); - let geoid = BaseGrid::plain(&geoid_header, Some(&geoid_grid), None)?; - - let c = Coor4D::geo(58.75, 08.25, 0., 0.); - assert_eq!(geoid.contains(&c, 0.0), false); - assert_eq!(geoid.contains(&c, 1.0), true); - - let n = geoid.at(&c, 1.0).unwrap(); - assert!((n[0] - 58.83).abs() < 0.1); - - let c = Coor4D::geo(53.25, 8.0, 0., 0.); - assert_eq!(geoid.contains(&c, 0.0), false); - assert_eq!(geoid.contains(&c, 1.0), true); - - let n = geoid.at(&c, 1.0).unwrap(); - assert!((n[0] - (53.25 + 0.08)).abs() < 0.001); - Ok(()) } } diff --git a/src/grid/ntv2/subgrid.rs b/src/grid/ntv2/subgrid.rs index 35734db..ffe4aee 100644 --- a/src/grid/ntv2/subgrid.rs +++ b/src/grid/ntv2/subgrid.rs @@ -6,7 +6,8 @@ pub(super) fn ntv2_subgrid(parser: &NTv2Parser, head_offset: usize) -> Result Result<(), Error> { let mut ctx = Plain::default(); - let op = ctx.op("gridshift grids=../../geodesy/datum/test.datum")?; + let op = ctx.op("gridshift grids=test.datum")?; let cph = Coor4D::geo(55., 12., 0., 0.); let mut data = [cph]; @@ -189,11 +189,29 @@ mod tests { Ok(()) } + #[test] + fn ntv2() -> Result<(), Error> { + let mut ctx = Plain::default(); + let op = ctx.op("gridshift grids=100800401.gsb")?; + let bcn = Coor2D::geo(41.3874, 2.1686); + let mut data = [bcn]; + + ctx.apply(op, Fwd, &mut data)?; + let res = data[0].to_geo(); + assert!((res[0] - 41.38627500250805).abs() < 1e-8); + assert!((res[1] - 2.167450821894838).abs() < 1e-8); + + ctx.apply(op, Inv, &mut data)?; + assert!((data[0][0] - bcn[0]).abs() < 1e-10); + assert!((data[0][1] - bcn[1]).abs() < 1e-10); + + Ok(()) + } + #[test] fn multiple_grids() -> Result<(), Error> { let mut ctx = Plain::default(); - let op = ctx - .op("gridshift grids=../../geodesy/datum/test.datum,../../geodesy/datum/test.datum")?; + let op = ctx.op("gridshift grids=test.datum, test.datum")?; let cph = Coor4D::geo(55., 12., 0., 0.); let mut data = [cph]; @@ -212,7 +230,7 @@ mod tests { #[test] fn fails_without_null_grid() -> Result<(), Error> { let mut ctx = Plain::default(); - let op = ctx.op("gridshift grids=../../geodesy/datum/test.datum")?; + let op = ctx.op("gridshift grids=test.datum")?; let ldn = Coor4D::geo(51.505, -0.09, 0., 0.); let mut data = [ldn]; @@ -228,7 +246,7 @@ mod tests { #[test] fn passes_with_null_grid() -> Result<(), Error> { let mut ctx = Plain::default(); - let op = ctx.op("gridshift grids=../../geodesy/datum/test.datum, @null")?; + let op = ctx.op("gridshift grids=test.datum, @null")?; let ldn = Coor4D::geo(51.505, -0.09, 0., 0.); let mut data = [ldn]; @@ -250,7 +268,7 @@ mod tests { #[test] fn optional_grid() -> Result<(), Error> { let mut ctx = Plain::default(); - let op = ctx.op("gridshift grids=@../../geodesy/datum/test_subset.datum, @missing.gsb, ../../geodesy/datum/test.datum")?; + let op = ctx.op("gridshift grids=@test_subset.datum, @missing.gsb, test.datum")?; // Copenhagen is outside of the (optional, but present, subset grid) let cph = Coor4D::geo(55., 12., 0., 0.); @@ -271,7 +289,6 @@ mod tests { let expected_correction = Coor4D([11.331, 55.971, 0., 0.]); ctx.apply(op, Fwd, &mut data)?; let correction = ((data[0] - haby) * Coor4D([3600., 3600., 3600., 3600.])).to_degrees(); - dbg!(correction); assert!((correction - expected_correction)[0].abs() < 1e-6); assert!((correction - expected_correction)[1].abs() < 1e-6);