Skip to content

Commit

Permalink
xy/set_xy for CoordinateSet. Tests via tmerc/merc
Browse files Browse the repository at this point in the history
  • Loading branch information
@busstoptaktik
busstoptaktik committed Apr 5, 2024
1 parent 94e6fea commit ef087cc
Show file tree
Hide file tree
Showing 4 changed files with 120 additions and 40 deletions.
21 changes: 21 additions & 0 deletions src/coordinate/mod.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,27 @@ pub trait CoordinateSet: CoordinateMetadata {
self.len() == 0
}

/// Replace the two first elements of the `index`th `CoordinateTuple` with `x` and `y`.
/// If the dimension in less than 2, fill the coordinate with `f64::NAN`.
/// For efficiency, consider implemented a specific implementation, when implementing
/// the CoordinateSet trait for a concrete data type: The provided version is not
/// very efficient.
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
let mut coord = self.get_coord(index);
coord.set_nth_unchecked(0, x);
coord.set_nth_unchecked(1, y);
self.set_coord(index, &coord);
}

/// Replace the two first elements of the `index`th `CoordinateTuple` with `x` and `y`.
/// If the dimension in less than 2, fill the coordinate with `f64::NAN`.
/// For efficiency, consider implemented a specific implementation, when implementing
/// the CoordinateSet trait for a concrete data type: The provided version is not
/// very efficient.
fn xy(&self, index: usize) -> (f64, f64) {
self.get_coord(index).xy()
}

/// Set all coordinate tuples in the set to NaN
fn stomp(&mut self) {
let nanny = Coor4D::nan();
Expand Down
60 changes: 60 additions & 0 deletions src/coordinate/set.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,12 @@ impl CoordinateSet for &mut [Coor4D] {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = *value;
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

impl<const N: usize> CoordinateSet for [Coor4D; N] {
Expand All @@ -30,6 +36,12 @@ impl<const N: usize> CoordinateSet for [Coor4D; N] {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = *value;
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

impl CoordinateSet for Vec<Coor4D> {
Expand All @@ -45,6 +57,12 @@ impl CoordinateSet for Vec<Coor4D> {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = *value;
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

// ----- CoordinateSet implementations for some Coor3D containers ------------
Expand All @@ -62,6 +80,12 @@ impl<const N: usize> CoordinateSet for [Coor3D; N] {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = Coor3D([value[0], value[1], value[2]]);
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

impl CoordinateSet for &mut [Coor3D] {
Expand All @@ -77,6 +101,12 @@ impl CoordinateSet for &mut [Coor3D] {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = Coor3D([value[0], value[1], value[2]]);
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

impl CoordinateSet for Vec<Coor3D> {
Expand All @@ -92,6 +122,12 @@ impl CoordinateSet for Vec<Coor3D> {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = Coor3D([value[0], value[1], value[2]]);
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

// ----- CoordinateSet implementations for some Coor2D containers ------------
Expand Down Expand Up @@ -127,6 +163,12 @@ impl<const N: usize> CoordinateSet for [Coor2D; N] {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = Coor2D([value[0], value[1]]);
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

impl CoordinateSet for &mut [Coor2D] {
Expand All @@ -142,6 +184,12 @@ impl CoordinateSet for &mut [Coor2D] {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = Coor2D([value[0], value[1]]);
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

impl CoordinateSet for Vec<Coor2D> {
Expand All @@ -157,6 +205,12 @@ impl CoordinateSet for Vec<Coor2D> {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = Coor2D([value[0], value[1]]);
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

/// User defined values for third and fourth coordinate dimension.
Expand Down Expand Up @@ -218,6 +272,12 @@ impl CoordinateSet for &mut [Coor32] {
fn set_coord(&mut self, index: usize, value: &Coor4D) {
self[index] = Coor32::raw(value[0], value[1]);
}
fn xy(&self, index: usize) -> (f64, f64) {
self[index].xy()
}
fn set_xy(&mut self, index: usize, x: f64, y: f64) {
self[index].set_xy(x, y);
}
}

// ----- Implementations: Coordinate Metadata ---------------------------------
Expand Down
26 changes: 12 additions & 14 deletions src/inner_op/merc.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,13 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
let mut successes = 0_usize;
let length = operands.len();
for i in 0..length {
let mut coord = operands.get_coord(i);
// Easting
coord[0] = (coord[0] - lon_0) * k_0 * a - x_0;
let (lon, lat) = operands.xy(i);

// Northing
let lat = coord[1] + lat_0;
coord[1] = a * k_0 * ellps.latitude_geographic_to_isometric(lat) - y_0;
let easting = (lon - lon_0) * k_0 * a - x_0;
let isometric = ellps.latitude_geographic_to_isometric(lat + lat_0);
let northing = a * k_0 * isometric - y_0;

operands.set_coord(i, &coord);
operands.set_xy(i, easting, northing);
successes += 1;
}

Expand All @@ -44,17 +42,17 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
let mut successes = 0_usize;
let length = operands.len();
for i in 0..length {
let mut coord = operands.get_coord(i);
let (mut x, mut y) = operands.xy(i);

// Easting -> Longitude
let x = coord[0] + x_0;
coord[0] = x / (a * k_0) - lon_0;
x += x_0;
let lon = x / (a * k_0) - lon_0;

// Northing -> Latitude
let y = coord[1] + y_0;
y += y_0;
let psi = y / (a * k_0);
coord[1] = ellps.latitude_isometric_to_geographic(psi) - lat_0;
operands.set_coord(i, &coord);
let lat = ellps.latitude_isometric_to_geographic(psi) - lat_0;
operands.set_xy(i, lon, lat);
successes += 1;
}

Expand Down Expand Up @@ -160,9 +158,9 @@ mod tests {
let definition = "merc lat_ts=56";
let op = ctx.op(definition)?;

// Validation value from PROJ: echo 12 55 0 0 | cct -d18 +proj=merc +lat_ts=56
let geo = [Coor4D::geo(55., 12., 0., 0.)];

// Validation value from PROJ: echo 12 55 0 0 | cct -d18 +proj=merc +lat_ts=56
let projected = [Coor4D::raw(
748_713.257_925_886_8,
4_106_573.862_841_270_4,
Expand Down
53 changes: 27 additions & 26 deletions src/inner_op/tmerc.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,14 +29,15 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
let range = 0..operands.len();
let mut successes = 0_usize;
for i in range {
let mut coord = operands.get_coord(i);
//let mut coord = operands.get_coord(i);
let (lon, lat) = operands.xy(i);

// --- 1. Geographical -> Conformal latitude, rotated longitude

// The conformal latitude
let lat = ellps.latitude_geographic_to_conformal(coord[1], conformal);
let lat = ellps.latitude_geographic_to_conformal(lat, conformal);
// The longitude as reckoned from the central meridian
let lon = coord[0] - lon_0;
let lon = lon - lon_0;

// --- 2. Conformal LAT, LNG -> complex spherical LAT

Expand All @@ -48,7 +49,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
// --- 3. Complex spherical N, E -> ellipsoidal normalized N, E

// Some numerical optimizations from PROJ modifications by Even Rouault,
let inv_denom_tan_lon = 1. / sin_lat.hypot(cos_lat_lon);
let inv_denom_tan_lon = sin_lat.hypot(cos_lat_lon).recip();
let tan_lon = sin_lon * cos_lat * inv_denom_tan_lon;
// Inverse Gudermannian, using the precomputed tan(lon)
let mut lon = tan_lon.asinh();
Expand All @@ -74,17 +75,18 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {

// Don't wanna play if we're too far from the center meridian
if lon.abs() > 2.623395162778 {
coord[0] = f64::NAN;
coord[1] = f64::NAN;
operands.set_xy(i, f64::NAN, f64::NAN);
continue;
}

// --- 4. ellipsoidal normalized N, E -> metric N, E

coord[0] = qs * lon + x_0; // Easting
coord[1] = qs * lat + zb; // Northing
let easting = qs * lon + x_0; // Easting
let northing = qs * lat + zb; // Northing

// Done!
operands.set_xy(i, easting, northing);
successes += 1;
operands.set_coord(i, &coord);
}

successes
Expand Down Expand Up @@ -118,17 +120,16 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
let range = 0..operands.len();
let mut successes = 0_usize;
for i in range {
let mut coord = operands.get_coord(i);
let (x, y) = operands.xy(i);

// --- 1. Normalize N, E

let mut lon = (coord[0] - x_0) / qs;
let mut lat = (coord[1] - zb) / qs;
let mut lon = (x - x_0) / qs;
let mut lat = (y - zb) / qs;

// Don't wanna play if we're too far from the center meridian
if lon.abs() > 2.623395162778 {
coord[0] = f64::NAN;
coord[1] = f64::NAN;
operands.set_xy(i, f64::NAN, f64::NAN);
continue;
}

Expand All @@ -151,10 +152,10 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {

let lon = angular::normalize_symmetric(lon + lon_0);
let lat = ellps.latitude_conformal_to_geographic(lat, conformal);
(coord[0], coord[1]) = (lon, lat);

// Done!
operands.set_xy(i, lon, lat);
successes += 1;
operands.set_coord(i, &coord);
}

successes
Expand Down Expand Up @@ -374,8 +375,8 @@ mod tests {
let geo = [
Coor2D::geo( 55., 12.),
Coor2D::geo(-55., 12.),
Coor2D::geo( 55., -6.),
Coor2D::geo(-55., -6.)
Coor2D::geo( 55., -6.),
Coor2D::geo(-55., -6.)
];

#[rustfmt::skip]
Expand Down Expand Up @@ -407,18 +408,18 @@ mod tests {

#[rustfmt::skip]
let geo = [
Coor4D::geo( 55., 12., 0., 0.),
Coor4D::geo(-55., 12., 0., 0.),
Coor4D::geo( 55., -6., 0., 0.),
Coor4D::geo(-55., -6., 0., 0.)
Coor2D::geo( 55., 12.),
Coor2D::geo(-55., 12.),
Coor2D::geo( 55., -6.,),
Coor2D::geo(-55., -6.,)
];

#[rustfmt::skip]
let projected = [
Coor4D::raw( 691_875.632_139_661, 1e7+6_098_907.825_005_012, 0., 0.),
Coor4D::raw( 691_875.632_139_661, 1e7-6_098_907.825_005_012, 0., 0.),
Coor4D::raw(-455_673.814_189_040, 1e7+6_198_246.671_090_279, 0., 0.),
Coor4D::raw(-455_673.814_189_040, 1e7-6_198_246.671_090_279, 0., 0.)
Coor2D::raw( 691_875.632_139_661, 1e7+6_098_907.825_005_012),
Coor2D::raw( 691_875.632_139_661, 1e7-6_098_907.825_005_012),
Coor2D::raw(-455_673.814_189_040, 1e7+6_198_246.671_090_279),
Coor2D::raw(-455_673.814_189_040, 1e7-6_198_246.671_090_279)
];

let mut operands = geo;
Expand Down

0 comments on commit ef087cc

Please sign in to comment.