1use crate::geo::Coord;
2use crate::{Result, TrailgenError};
3use serde::{Deserialize, Serialize};
4
5const WEB_MERCATOR_R_M: f64 = 6_378_137.0;
6
7#[derive(Clone, Copy, Debug, Eq, PartialEq)]
8pub enum CrsVerdict {
9 AssumedGeographic,
10 Geographic(GeodeticDatum),
11 WebMercator,
12 Utm(UtmCrs),
13}
14
15#[derive(Clone, Copy, Debug, Eq, PartialEq)]
16pub enum CoordProjector {
17 Identity,
18 WebMercator,
19 Utm(UtmCrs),
20}
21
22#[derive(Clone, Copy, Debug, Eq, PartialEq, Serialize, Deserialize)]
23#[serde(rename_all = "kebab-case")]
24pub enum GeodeticDatum {
25 Wgs84,
26 Nad83,
27}
28
29impl GeodeticDatum {
30 #[must_use]
31 pub fn from_normalized_geographic(normalized: &str) -> Option<Self> {
32 if normalized.contains("EPSG4326")
33 || normalized.contains("OGC13CRS84")
34 || normalized.contains("OGC14CRS84")
35 || normalized.contains("CRS84")
36 || normalized.contains("WGS84")
37 || normalized.contains("WGS1984")
38 {
39 Some(Self::Wgs84)
40 } else if normalized.contains("EPSG4269")
41 || normalized.contains("NAD83")
42 || normalized.contains("NORTHAMERICANDATUM1983")
43 {
44 Some(Self::Nad83)
45 } else {
46 None
47 }
48 }
49
50 const fn ellipsoid(self) -> Ellipsoid {
51 match self {
52 Self::Wgs84 => Ellipsoid {
53 a: 6_378_137.0,
54 inv_f: 298.257_223_563,
55 },
56 Self::Nad83 => Ellipsoid {
57 a: 6_378_137.0,
58 inv_f: 298.257_222_101,
59 },
60 }
61 }
62}
63
64#[derive(Clone, Copy, Debug, Eq, PartialEq, Serialize, Deserialize)]
65#[serde(rename_all = "kebab-case")]
66pub enum UtmHemisphere {
67 North,
68 South,
69}
70
71#[derive(Clone, Copy, Debug, Eq, PartialEq, Serialize, Deserialize)]
72pub struct UtmCrs {
73 pub datum: GeodeticDatum,
74 pub zone: u8,
75 pub hemisphere: UtmHemisphere,
76}
77
78impl UtmCrs {
79 #[must_use]
80 pub const fn from_parts(
81 datum: GeodeticDatum,
82 zone: u8,
83 hemisphere: UtmHemisphere,
84 ) -> Option<Self> {
85 if zone >= 1 && zone <= 60 {
86 Some(Self {
87 datum,
88 zone,
89 hemisphere,
90 })
91 } else {
92 None
93 }
94 }
95
96 #[must_use]
97 pub fn from_epsg(epsg: u16) -> Option<Self> {
98 match epsg {
99 32601..=32660 => Self::from_parts(
100 GeodeticDatum::Wgs84,
101 u8::try_from(epsg - 32600).ok()?,
102 UtmHemisphere::North,
103 ),
104 32701..=32760 => Self::from_parts(
105 GeodeticDatum::Wgs84,
106 u8::try_from(epsg - 32700).ok()?,
107 UtmHemisphere::South,
108 ),
109 26901..=26923 => Self::from_parts(
110 GeodeticDatum::Nad83,
111 u8::try_from(epsg - 26900).ok()?,
112 UtmHemisphere::North,
113 ),
114 _ => None,
115 }
116 }
117
118 #[must_use]
119 pub fn from_normalized_srs(normalized: &str) -> Option<Self> {
120 (1_u16..=60).find_map(|zone| {
121 let zone_u8 = u8::try_from(zone).ok()?;
122 [
123 (
124 format!("EPSG326{zone:02}"),
125 Self::from_parts(GeodeticDatum::Wgs84, zone_u8, UtmHemisphere::North)?,
126 ),
127 (
128 format!("EPSG327{zone:02}"),
129 Self::from_parts(GeodeticDatum::Wgs84, zone_u8, UtmHemisphere::South)?,
130 ),
131 (
132 format!("EPSG269{zone:02}"),
133 Self::from_parts(GeodeticDatum::Nad83, zone_u8, UtmHemisphere::North)?,
134 ),
135 (
136 format!("WGS84UTMZONE{zone}N"),
137 Self::from_parts(GeodeticDatum::Wgs84, zone_u8, UtmHemisphere::North)?,
138 ),
139 (
140 format!("WGS84UTMZONE{zone}S"),
141 Self::from_parts(GeodeticDatum::Wgs84, zone_u8, UtmHemisphere::South)?,
142 ),
143 (
144 format!("NAD83UTMZONE{zone}N"),
145 Self::from_parts(GeodeticDatum::Nad83, zone_u8, UtmHemisphere::North)?,
146 ),
147 (
148 format!("NAD1983UTMZONE{zone}N"),
149 Self::from_parts(GeodeticDatum::Nad83, zone_u8, UtmHemisphere::North)?,
150 ),
151 ]
152 .into_iter()
153 .find_map(|(needle, crs)| normalized.contains(&needle).then_some(crs))
154 })
155 }
156
157 const fn false_northing_m(self) -> f64 {
158 match self.hemisphere {
159 UtmHemisphere::North => 0.0,
160 UtmHemisphere::South => 10_000_000.0,
161 }
162 }
163
164 fn λ0(self) -> f64 {
165 ((f64::from(self.zone) - 1.0).mul_add(6.0, -177.0)).to_radians()
166 }
167}
168
169#[derive(Clone, Copy)]
170struct Ellipsoid {
171 a: f64,
172 inv_f: f64,
173}
174
175impl Ellipsoid {
176 fn f(self) -> f64 {
177 1.0 / self.inv_f
178 }
179}
180
181impl CoordProjector {
182 #[must_use]
183 pub fn project(self, x: f64, y: f64, ele: Option<f64>) -> Coord {
184 match self {
185 Self::Identity => Coord {
186 lon: x,
187 lat: y,
188 ele,
189 },
190 Self::WebMercator => web_mercator_to_wgs84(x, y, ele),
191 Self::Utm(crs) => utm_to_geographic(x, y, ele, crs),
192 }
193 }
194}
195
196#[must_use]
197pub(crate) fn wgs84_to_web_mercator(coord: Coord) -> (f64, f64) {
198 let lat = coord.lat.clamp(-85.051_128_78, 85.051_128_78).to_radians();
199 (
200 WEB_MERCATOR_R_M * coord.lon.to_radians(),
201 WEB_MERCATOR_R_M * (std::f64::consts::FRAC_PI_4 + lat / 2.0).tan().ln(),
202 )
203}
204
205#[must_use]
206#[allow(clippy::many_single_char_names, clippy::suboptimal_flops)]
207pub fn wgs84_to_utm(coord: Coord, zone: u8, north: bool) -> Option<(f64, f64)> {
208 let hemisphere = if north {
209 UtmHemisphere::North
210 } else {
211 UtmHemisphere::South
212 };
213 geographic_to_utm(
214 coord,
215 UtmCrs::from_parts(GeodeticDatum::Wgs84, zone, hemisphere)?,
216 )
217}
218
219#[must_use]
220#[allow(clippy::many_single_char_names, clippy::suboptimal_flops)]
221pub fn geographic_to_utm(coord: Coord, crs: UtmCrs) -> Option<(f64, f64)> {
222 if !coord.lon.is_finite() || !coord.lat.is_finite() {
223 return None;
224 }
225 let φ = coord.lat.to_radians();
226 let λ = coord.lon.to_radians();
227 let λ0 = crs.λ0();
228 let ellipsoid = crs.datum.ellipsoid();
229 let a = ellipsoid.a;
230 let f = ellipsoid.f();
231 let k0 = 0.9996;
232 let e2 = f * (2.0 - f);
233 let e4 = e2 * e2;
234 let e6 = e4 * e2;
235 let ep2 = e2 / (1.0 - e2);
236 let sinφ = φ.sin();
237 let cosφ = φ.cos();
238 let tanφ = φ.tan();
239 let n = a / (1.0 - e2 * sinφ * sinφ).sqrt();
240 let t = tanφ * tanφ;
241 let c = ep2 * cosφ * cosφ;
242 let aa = cosφ * (λ - λ0);
243 let m = a
244 * ((1.0 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0) * φ
245 - (3.0 * e2 / 8.0 + 3.0 * e4 / 32.0 + 45.0 * e6 / 1024.0) * (2.0 * φ).sin()
246 + (15.0 * e4 / 256.0 + 45.0 * e6 / 1024.0) * (4.0 * φ).sin()
247 - (35.0 * e6 / 3072.0) * (6.0 * φ).sin());
248 let easting = 500_000.0
249 + k0 * n
250 * (aa
251 + (1.0 - t + c) * aa.powi(3) / 6.0
252 + (5.0 - 18.0 * t + t * t + 72.0 * c - 58.0 * ep2) * aa.powi(5) / 120.0);
253 let mut northing = k0
254 * (m + n
255 * tanφ
256 * (aa.powi(2) / 2.0
257 + (5.0 - t + 9.0 * c + 4.0 * c * c) * aa.powi(4) / 24.0
258 + (61.0 - 58.0 * t + t * t + 600.0 * c - 330.0 * ep2) * aa.powi(6) / 720.0));
259 northing += crs.false_northing_m();
260 (easting.is_finite() && northing.is_finite()).then_some((easting, northing))
261}
262
263#[derive(Clone, Copy, Debug, Eq, PartialEq)]
264pub enum VectorCrsKind {
265 GeoJson,
266 ShapefilePrj,
267}
268
269impl VectorCrsKind {
270 const fn label(self) -> &'static str {
271 match self {
272 Self::GeoJson => "GeoJSON CRS",
273 Self::ShapefilePrj => "shapefile .prj CRS",
274 }
275 }
276}
277
278pub fn validate_crs_name(kind: VectorCrsKind, name: &str) -> Result<CrsVerdict> {
279 let normalized = normalize(name);
280 match UtmCrs::from_normalized_srs(&normalized) {
281 Some(crs) => Ok(CrsVerdict::Utm(crs)),
282 None if let Some(datum) = GeodeticDatum::from_normalized_geographic(&normalized) => {
283 Ok(CrsVerdict::Geographic(datum))
284 }
285 None if is_web_mercator(&normalized) => Ok(CrsVerdict::WebMercator),
286 None => Err(TrailgenError::InvalidData(format!(
287 "{} {name:?} is not supported; reproject input to geographic lon/lat WGS84/NAD83 (EPSG:4326/4269/CRS84), EPSG:3857 Web Mercator, or WGS84/NAD83 UTM (EPSG:326xx/327xx/269xx), before ingestion",
288 kind.label()
289 ))),
290 }
291}
292
293pub fn validate_prj_wkt(wkt: &str) -> Result<CrsVerdict> {
294 let normalized = normalize(wkt);
295 if is_web_mercator(&normalized) {
296 return Ok(CrsVerdict::WebMercator);
297 }
298 if let Some(crs) = UtmCrs::from_normalized_srs(&normalized) {
299 return Ok(CrsVerdict::Utm(crs));
300 }
301 if is_projected(&normalized) {
302 return Err(TrailgenError::InvalidData(
303 "shapefile .prj advertises an unsupported projected CRS; reproject input to geographic lon/lat WGS84/NAD83 (EPSG:4326/4269/CRS84), EPSG:3857 Web Mercator, or WGS84/NAD83 UTM (EPSG:326xx/327xx/269xx), before ingestion"
304 .to_owned(),
305 ));
306 }
307 validate_crs_name(VectorCrsKind::ShapefilePrj, wkt)
308}
309
310#[must_use]
311pub const fn projector(verdict: CrsVerdict) -> CoordProjector {
312 match verdict {
313 CrsVerdict::AssumedGeographic | CrsVerdict::Geographic(_) => CoordProjector::Identity,
314 CrsVerdict::WebMercator => CoordProjector::WebMercator,
315 CrsVerdict::Utm(crs) => CoordProjector::Utm(crs),
316 }
317}
318
319fn normalize(raw: &str) -> String {
320 raw.chars()
321 .filter(char::is_ascii_alphanumeric)
322 .flat_map(char::to_uppercase)
323 .collect()
324}
325
326fn is_web_mercator(normalized: &str) -> bool {
327 normalized.contains("EPSG3857")
328 || normalized.contains("EPSG900913")
329 || normalized.contains("WEBMERCATOR")
330 || normalized.contains("PSEUDOMERCATOR")
331 || normalized.contains("WGS84PSEUDOMERCATOR")
332}
333
334fn is_projected(normalized: &str) -> bool {
335 normalized.contains("PROJCS")
336 || normalized.contains("PROJCRS")
337 || normalized.contains("PROJECTION")
338}
339
340fn web_mercator_to_wgs84(x_m: f64, y_m: f64, ele: Option<f64>) -> Coord {
341 let lon = (x_m / WEB_MERCATOR_R_M).to_degrees();
342 let lat = 2.0f64
343 .mul_add(
344 (y_m / WEB_MERCATOR_R_M).exp().atan(),
345 -std::f64::consts::FRAC_PI_2,
346 )
347 .to_degrees();
348 Coord { lon, lat, ele }
349}
350
351#[allow(clippy::many_single_char_names, clippy::suboptimal_flops)]
352fn utm_to_geographic(easting_m: f64, northing_m: f64, ele: Option<f64>, crs: UtmCrs) -> Coord {
353 let x = easting_m - 500_000.0;
354 let y = northing_m - crs.false_northing_m();
355 let ellipsoid = crs.datum.ellipsoid();
356 let a = ellipsoid.a;
357 let f = ellipsoid.f();
358 let k0 = 0.9996;
359 let e2 = f * (2.0 - f);
360 let e4 = e2 * e2;
361 let e6 = e4 * e2;
362 let ep2 = e2 / (1.0 - e2);
363 let e1 = (1.0 - (1.0 - e2).sqrt()) / (1.0 + (1.0 - e2).sqrt());
364 let e1_2 = e1 * e1;
365 let e1_3 = e1_2 * e1;
366 let e1_4 = e1_2 * e1_2;
367 let μ = y / (k0 * a * (1.0 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0));
368 let φ1 = μ
369 + (3.0 * e1 / 2.0 - 27.0 * e1_3 / 32.0) * (2.0 * μ).sin()
370 + (21.0 * e1_2 / 16.0 - 55.0 * e1_4 / 32.0) * (4.0 * μ).sin()
371 + (151.0 * e1_3 / 96.0) * (6.0 * μ).sin()
372 + (1097.0 * e1_4 / 512.0) * (8.0 * μ).sin();
373 let sinφ1 = φ1.sin();
374 let cosφ1 = φ1.cos();
375 let tanφ1 = φ1.tan();
376 let n1 = a / (1.0 - e2 * sinφ1 * sinφ1).sqrt();
377 let r1 = a * (1.0 - e2) / (1.0 - e2 * sinφ1 * sinφ1).powf(1.5);
378 let t1 = tanφ1 * tanφ1;
379 let c1 = ep2 * cosφ1 * cosφ1;
380 let d = x / (n1 * k0);
381 let lat = φ1
382 - (n1 * tanφ1 / r1)
383 * (d.powi(2) / 2.0
384 - (5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1 * c1 - 9.0 * ep2) * d.powi(4) / 24.0
385 + (61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1 * t1 - 252.0 * ep2 - 3.0 * c1 * c1)
386 * d.powi(6)
387 / 720.0);
388 let lon = crs.λ0()
389 + (d - (1.0 + 2.0 * t1 + c1) * d.powi(3) / 6.0
390 + (5.0 - 2.0 * c1 + 28.0 * t1 - 3.0 * c1 * c1 + 8.0 * ep2 + 24.0 * t1 * t1)
391 * d.powi(5)
392 / 120.0)
393 / cosφ1;
394 Coord {
395 lon: lon.to_degrees(),
396 lat: lat.to_degrees(),
397 ele,
398 }
399}