use super::Coord; use geo_macros::Prj; use num_traits::{AsPrimitive, FromPrimitive, Num}; use proj::{Proj, ProjError}; use thiserror::Error; type Lat = f64; type Lon = f64; #[derive(Clone, Copy)] pub struct Range(pub f64, pub f64); impl Range { pub fn key_points(&self, max_points: usize) -> Vec { if max_points == 0 { return vec![]; } let range = (self.0.min(self.1) as f64, self.1.max(self.0) as f64); assert!(!(range.0.is_nan() || range.1.is_nan())); if (range.0 - range.1).abs() < std::f64::EPSILON { return vec![range.0 as f64]; } let mut scale = (10f64).powf((range.1 - range.0).log(10.0).floor()); // The value granularity controls how we round the values. // To avoid generating key points like 1.00000000001, we round to the nearest multiple of the // value granularity. // By default, we make the granularity as the 1/10 of the scale. let mut value_granularity = scale / 10.0; fn rem_euclid(a: f64, b: f64) -> f64 { let ret = if b > 0.0 { a - (a / b).floor() * b } else { a - (a / b).ceil() * b }; if (ret - b).abs() < std::f64::EPSILON { 0.0 } else { ret } } // At this point we need to make sure that the loop invariant: // The scale must yield number of points than requested if 1 + ((range.1 - range.0) / scale).floor() as usize > max_points { scale *= 10.0; value_granularity *= 10.0; } 'outer: loop { let old_scale = scale; for nxt in [2.0, 5.0, 10.0].iter() { let mut new_left = range.0 - rem_euclid(range.0, old_scale / nxt); if new_left < range.0 { new_left += old_scale / nxt; } let new_right = range.1 - rem_euclid(range.1, old_scale / nxt); let npoints = 1.0 + ((new_right - new_left) / old_scale * nxt); if npoints.round() as usize > max_points { break 'outer; } scale = old_scale / nxt; } scale = old_scale / 10.0; value_granularity /= 10.0; } let mut ret = vec![]; // In some extreme cases, left might be too big, so that (left + scale) - left == 0 due to // floating point error. // In this case, we may loop forever. To avoid this, we need to use two variables to store // the current left value. So we need keep a left_base and a left_relative. let left = { let mut value = range.0 - rem_euclid(range.0, scale); if value < range.0 { value += scale; } value }; let left_base = (left / value_granularity).floor() * value_granularity; let mut left_relative = left - left_base; let right = range.1 - rem_euclid(range.1, scale); while (right - left_relative - left_base) >= -std::f64::EPSILON { let new_left_relative = (left_relative / value_granularity).round() * value_granularity; if new_left_relative < 0.0 { left_relative += value_granularity; } ret.push((left_relative + left_base) as f64); left_relative += scale; } return ret; } } #[derive(Error, Debug)] pub enum CoordError { #[error("")] ProjError { #[from] source: ProjError, }, } impl + Num> From<(T, T)> for Range { fn from(value: (T, T)) -> Self { let value = (value.0.as_(), value.1.as_()); let (_min, _max) = (value.0.min(value.1), value.0.max(value.1)); Self(_min, _max) } } pub struct LatLonCoord { // actual: ((f64, f64), (f64, f64)), actual: (Range, Range), logical: (Range, Range), pcs: PCS, } struct PCS { pub lon_range: Range, pub lat_range: Range, pub proj_param: T, transformer: Proj, } unsafe impl Sync for PCS {} unsafe impl Send for PCS {} impl LatLonCoord { /// Creates a new `LatLonCoord` instance. /// /// # Arguments /// /// * `lon` - An optional longitude range. /// * `lat` - An optional latitude range. /// * `actual` - A tuple containing the actual ranges. /// * `project` - A projection. /// /// # Returns /// /// A new `LatLonCoord` instance. pub fn new( lon: Option, lat: Option, actual: ((i32, i32), (i32, i32)), // actual: (Range, Range), project: T, ) -> Self { let pcs = PCS::new(project, lon, lat); let _box = pcs.bbox().unwrap(); Self { actual: (Range::from(actual.0), Range::from(actual.1)), pcs: pcs, logical: _box, } } pub fn set_actual(&mut self, actual: ((i32, i32), (i32, i32))) { self.actual = (Range::from(actual.0), Range::from(actual.1)); } pub fn lon_range(&self) -> Range { self.pcs.lon_range } pub fn lat_range(&self) -> Range { self.pcs.lat_range } } unsafe impl Sync for LatLonCoord {} unsafe impl Send for LatLonCoord {} impl Coord for LatLonCoord { fn map(&self, axis_1: f64, axis_2: f64) -> super::ScreenCoord { let point = self.pcs.map((axis_1, axis_2)); let logical_dim1_span = self.logical.0 .1 - self.logical.0 .0; let dim1_rate = (point.0 - self.logical.0 .0) / logical_dim1_span; let logical_dim2_span = self.logical.1 .1 - self.logical.1 .0; let dim2_rate = (point.1 - self.logical.1 .0) / logical_dim2_span; ( (dim1_rate * (self.actual.0 .1 - self.actual.0 .0) as f64) + self.actual.0 .0, (dim2_rate * (self.actual.1 .1 - self.actual.1 .0) as f64) + self.actual.1 .0, ) } fn dim1_range(&self) -> (f64, f64) { let v = self.lon_range(); (v.0, v.1) } fn dim2_range(&self) -> (f64, f64) { let v = self.lat_range(); (v.0, v.1) } } #[derive(Clone, Debug)] pub enum Projection { PlateCarree, LambertConformal, LambertCylindrical, Mercator, } /// A trait for defining a projection. pub trait ProjectionS: Sync + Send { /// Returns a proj-string of the projection. fn build(&self) -> String; /// Returns the logical range of the projection. /// In common, different projections have different logical ranges. /// # Arguments /// /// * `lon_range` - An optional longitude range. /// * `lat_range` - An optional latitude range. /// /// # Returns /// /// A tuple containing the longitude and latitude ranges. fn logic_range(&self, lon_range: Option, lat_range: Option) -> (Range, Range); } impl PCS { fn new(proj_param: T, lon_range: Option, lat_range: Option) -> Self { let (lon_range, lat_range) = proj_param.logic_range(lon_range, lat_range); Self { lon_range, lat_range, transformer: Proj::new(proj_param.build().as_str()).unwrap(), proj_param: proj_param, } } fn bbox(&self) -> Result<(Range, Range), CoordError> { let _proj_transformer = &self.transformer; let lb = (self.lon_range.0.to_radians(), self.lat_range.0.to_radians()); let rt = (self.lon_range.1.to_radians(), self.lat_range.1.to_radians()); let bl = _proj_transformer.convert(lb)?; let rt = _proj_transformer.convert(rt)?; Ok((Range::from((bl.0, rt.0)), Range::from((bl.1, rt.1)))) } fn map(&self, lon_lat: (Lat, Lon)) -> (f64, f64) { let _proj_transformer = &self.transformer; let _lon_lat = _proj_transformer .convert((lon_lat.0.to_radians(), lon_lat.1.to_radians())) .unwrap(); _lon_lat } } #[derive(Prj)] /// A struct representing the Mercator projection. pub struct Mercator { /// The central longitude of the projection. pub central_lon: f64, /// The minimum latitude of the projection. pub min_latitude: f64, /// The maximum latitude of the projection. pub max_latitude: f64, /// The false easting of the projection. pub false_easting: f64, /// The false northing of the projection. pub false_northing: f64, /// The latitude of true scale of the projection. pub latitude_true_scale: f64, } fn proj_string<'a>(vs: Vec<(&'a str, &'a str)>) -> String { vs.into_iter() .map(|(option, value)| format!("+{}={}", option, value)) .collect::>() .join(" ") } impl Mercator { /// Creates a new instance of the `Mercator` projection with default values. pub fn new() -> Self { Self { central_lon: 0.0, min_latitude: -82.0, max_latitude: 82.0, false_easting: 0.0, false_northing: 0.0, latitude_true_scale: 0.0, } } } impl ProjectionS for Mercator { fn build(&self) -> String { let _central_lon = format!("{:.1}", &self.central_lon); let _false_easting = format!("{:.1}", &self.false_easting); let _false_northing = format!("{:.1}", &self.false_northing); let _ts = format!("{:.1}", &self.latitude_true_scale); let input = vec![ ("proj", "merc"), ("lon_0", _central_lon.as_str()), ("x_0", _false_easting.as_str()), ("y_0", _false_northing.as_str()), ("lat_ts", _ts.as_str()), ("units", "m"), ("ellps", "GRS80"), ]; let _proj_string = proj_string(input); _proj_string } fn logic_range(&self, lon_range: Option, lat_range: Option) -> (Range, Range) { let lon_range = lon_range.unwrap_or(Range { 0: -180f64, 1: 180f64, }); let lat_range = lat_range.unwrap_or(Range { 0: self.min_latitude, 1: self.max_latitude, }); let lat_range = Range { 0: lat_range.0.max(self.min_latitude), 1: lat_range.1.min(self.max_latitude), }; (lon_range, lat_range) } }