Map projections

use projections.qweruse proj_point, proj_points from projections      // Quillo 0.6: only these

proj4rs = { version = "=0.3.0", features = ["crs-definitions"] }crs-definitions = "=0.6.0"serde_json = { version = "1", features = ["preserve_order"] }

// A planning site: parcels from the national grid shown on a web map, GPS points put on the grid, distances in metres.////   quillo run projections.qwe      then open http://localhost:8080/pointsuse "projections.qwer"// Big Ben, in GPS degrees (longitude first).let BIG_BEN = [-0.124625, 51.500729]route GET /points {    let zone = proj_utm_zone(BIG_BEN[1], BIG_BEN[0])?    {        web_mercator: proj_point(BIG_BEN[0], BIG_BEN[1], "EPSG:4326", "EPSG:3857")?,        british_grid: proj_point(BIG_BEN[0], BIG_BEN[1], "EPSG:4326", "EPSG:27700")?,        back: proj_point(530268.0, 179640.0, "EPSG:27700", "EPSG:4326")?,        zone: zone,        utm: proj_point(BIG_BEN[0], BIG_BEN[1], "EPSG:4326", zone)?,        dutch: proj_points([[4.8897, 52.3731], [5.1214, 52.0907]], "EPSG:4326", "EPSG:28992")?,        crs: proj_crs("EPSG:27700")?,        unknown: proj_point(0.0, 0.0, "EPSG:1", "EPSG:4326").error,    }}// A parcel drawn on the British grid, in GPS degrees for Leaflet.route GET /parcel = proj_geojson({    type: "Feature",    properties: { parcel: "TQ 3079" },    geometry: { type: "Polygon", coordinates: [[[530000.0, 179000.0], [531000.0, 179000.0], [531000.0, 180000.0], [530000.0, 180000.0], [530000.0, 179000.0]]] },}, "EPSG:27700", "EPSG:4326")?

//! Coordinates moved between map projections by EPSG code or PROJ string, offline: GPS degrees (EPSG:4326) to Web Mercator (3857), UTM zones, the British National Grid (27700), Dutch RD (28992), Swiss LV95 (2056) and some 6,000 others, points and whole GeoJSON.//! version: 1.0.0//!//! [dependencies]//! proj4rs = { version = "=0.3.0", features = ["crs-definitions"] }//! crs-definitions = "=0.6.0"//! serde_json = { version = "1", features = ["preserve_order"] }use proj4rs::Proj;use serde_json::{json, Value};fn code_of(crs: &str) -> Option<u16> {    let t = crs.trim();    let up = t.to_ascii_uppercase();    let num = up.strip_prefix("EPSG:").or_else(|| up.strip_prefix("EPSG::")).unwrap_or(&up);    num.trim().parse::<u16>().ok()}fn proj(crs: &str) -> Result<Proj, String> {    let t = crs.trim();    if t.starts_with('+') {        return Proj::from_proj_string(t).map_err(|e| format!("`{}`: {}", t, e));    }    if t.eq_ignore_ascii_case("WGS84") {        return Proj::from_proj_string("+proj=longlat +datum=WGS84 +no_defs").map_err(|e| e.to_string());    }    let code = code_of(t).ok_or_else(|| format!("`{}` is neither an EPSG code (\"EPSG:3857\") nor a PROJ string (\"+proj=...\")", t))?;    let def = crs_definitions::from_code(code).ok_or_else(|| format!("EPSG:{} is not one this knows", code))?;    Proj::from_proj_string(def.proj4).map_err(|e| format!("EPSG:{}: {}", code, e))}struct Pair { from: Proj, to: Proj }impl Pair {    fn new(from: &str, to: &str) -> Result<Pair, String> {        Ok(Pair { from: proj(from)?, to: proj(to)? })    }    fn point(&self, x: f64, y: f64, z: f64) -> Result<(f64, f64, f64), String> {        let mut p = if self.from.is_latlong() { (x.to_radians(), y.to_radians(), z) } else { (x, y, z) };        proj4rs::transform::transform(&self.from, &self.to, &mut p).map_err(|e| format!("({}, {}): {}", x, y, e))?;        if !p.0.is_finite() || !p.1.is_finite() {            return Err(format!("({}, {}) has no place in the target projection", x, y));        }        Ok(if self.to.is_latlong() {            (round(p.0.to_degrees(), 1e9), round(p.1.to_degrees(), 1e9), p.2)        } else {            (round(p.0, 1e4), round(p.1, 1e4), p.2)        })    }}fn round(v: f64, k: f64) -> f64 {    (v * k).round() / k}fn walk(v: &mut Value, pair: &Pair) -> Result<(), String> {    match v {        Value::Array(a) if a.len() >= 2 && a[0].is_number() => {            let x = a[0].as_f64().unwrap_or(0.0);            let y = a[1].as_f64().ok_or("a coordinate is a number")?;            let z = a.get(2).and_then(|z| z.as_f64()).unwrap_or(0.0);            let (px, py, pz) = pair.point(x, y, z)?;            a[0] = json!(px);            a[1] = json!(py);            if a.len() >= 3 { a[2] = json!(round(pz, 1e4)) }        }        Value::Array(a) => for x in a { walk(x, pair)? },        Value::Object(o) => {            for key in ["coordinates", "geometry", "geometries", "features"] {                if let Some(x) = o.get_mut(key) { walk(x, pair)? }            }            o.remove("bbox");        }        _ => {}    }    Ok(())}/// One point moved from projection `from` to `to` (`"EPSG:4326"`, `"3857"`, `"+proj=utm +zone=31 +datum=WGS84"`): `[x, y]`. Degrees go longitude first, as GeoJSON has them; metres come rounded to 0.1 mm.pub fn proj_point(x: f64, y: f64, from: &str, to: &str) -> Result<Vec<f64>, String> {    let (px, py, _) = Pair::new(from, to)?.point(x, y, 0.0)?;    Ok(vec![px, py])}/// Many points (`[[x, y], ...]`, an optional height third) moved from `from` to `to` at once.pub fn proj_points(points: Vec<Vec<f64>>, from: &str, to: &str) -> Result<Vec<Vec<f64>>, String> {    let pair = Pair::new(from, to)?;    points.iter().map(|p| {        if p.len() < 2 { return Err("each point is [x, y]".to_string()) }        let (x, y, z) = pair.point(p[0], p[1], p.get(2).copied().unwrap_or(0.0))?;        Ok(if p.len() >= 3 { vec![x, y, round(z, 1e4)] } else { vec![x, y] })    }).collect()}/// A GeoJSON geometry, feature or FeatureCollection with every position moved from `from` to `to`: data from a national grid to WGS 84 for a web map, or back. Any `bbox` is dropped, being no longer right.pub fn proj_geojson(geojson: serde_json::Value, from: &str, to: &str) -> Result<Value, String> {    let pair = Pair::new(from, to)?;    let mut g = geojson;    walk(&mut g, &pair)?;    Ok(g)}/// What an EPSG code is: `{ code, name, proj4, geographic, unit }` (`"degrees"` or `"m"`...).pub fn proj_crs(crs: &str) -> Result<Value, String> {    let code = code_of(crs).ok_or_else(|| format!("`{}` is not an EPSG code", crs))?;    let def = crs_definitions::from_code(code).ok_or_else(|| format!("EPSG:{} is not one this knows", code))?;    let p = Proj::from_proj_string(def.proj4).map_err(|e| format!("EPSG:{}: {}", code, e))?;    let name = def.wkt.split('"').nth(1).unwrap_or("").to_string();    Ok(json!({        "code": format!("EPSG:{}", code),        "name": name,        "proj4": def.proj4,        "geographic": p.is_latlong(),        "unit": if p.is_latlong() { "degrees" } else { p.units() },    }))}/// The WGS 84 UTM zone of a point, as an EPSG code: `"EPSG:32631"` in the north, `"EPSG:327xx"` in the south, Norway's and Svalbard's exceptions counted. Metres for measuring a small area.pub fn proj_utm_zone(lat: f64, lon: f64) -> Result<String, String> {    if !(-80.0..=84.0).contains(&lat) || !(-180.0..=180.0).contains(&lon) {        return Err(format!("({}, {}) is outside the UTM zones (80°S to 84°N)", lat, lon));    }    let mut zone = (((lon + 180.0) / 6.0).floor() as i64 % 60) + 1;    if (56.0..64.0).contains(&lat) && (3.0..12.0).contains(&lon) { zone = 32 }    if (72.0..=84.0).contains(&lat) {        zone = if (0.0..9.0).contains(&lon) { 31 } else if (9.0..21.0).contains(&lon) { 33 } else if (21.0..33.0).contains(&lon) { 35 } else if (33.0..42.0).contains(&lon) { 37 } else { zone };    }    Ok(format!("EPSG:{}", if lat >= 0.0 { 32600 + zone } else { 32700 + zone }))}