diff --git a/crates/engine/src/cmd/qgis.rs b/crates/engine/src/cmd/qgis.rs index 7039437f9..4ab560396 100644 --- a/crates/engine/src/cmd/qgis.rs +++ b/crates/engine/src/cmd/qgis.rs @@ -84,13 +84,21 @@ fn import_project(s: &mut Session, p: &Value) -> CmdResult { let m = s.map_mut(); m.basemap = rustgis_map::Basemap::Blank; m.graticule = false; + let degrees = project.crs.as_deref().and_then(Geographic::of); if let Some(crs) = &project.crs { let code = match crs.as_str() { "EPSG:3857" | "EPSG:900913" => "mercator", + // QGIS draws a project in degrees as plate carrée: longitude and latitude as x and y. + _ if degrees.is_some() => "equirectangular", other => other, }; - if let Err(e) = s.run_part("map.projection", &json!({"type": code})) { - report.approximated("", format!("the project's coordinate system {crs} ({e}), drawn in Web Mercator")); + match s.run_part("map.projection", &json!({"type": code})) { + Err(e) => report.approximated("", format!("the project's coordinate system {crs} ({e}), drawn in Web Mercator")), + Ok(_) if degrees.is_some_and(|g| !g.wgs84 || !g.greenwich) => { + let meridian = if degrees.is_some_and(|g| !g.greenwich) { ", its longitudes taken as from Greenwich" } else { "" }; + report.approximated("", format!("the project's coordinate system {crs}, drawn in Equirectangular on WGS 84{meridian}")); + } + Ok(_) => {} } } // Bottom first, so each is added on top of the ones below it. @@ -102,7 +110,7 @@ fn import_project(s: &mut Session, p: &Value) -> CmdResult { Err(e) => report.skipped(&layer.name, format!("the layer: {e}")), } } - go_to_view(s, &project); + go_to_view(s, &project, degrees); let notes: Vec = report.notes.iter().map(note_json).collect(); if !report.notes.is_empty() { let left = report.notes.iter().filter(|n| n.kind == NoteKind::Skipped).count(); @@ -302,12 +310,21 @@ pub(crate) fn web_settings(s: &mut Session, added: &Value, before: Option, } /// Goes to the project's saved view, or to its layers when that view shows none of them. -fn go_to_view(s: &mut Session, project: &qgis::Project) { +fn go_to_view(s: &mut Session, project: &qgis::Project, degrees: Option) { let data = s.map().full_extent(); - let saved = project.extent.and_then(|[x0, y0, x1, y1]| { - let m = s.map(); - let corners = [(x0, y0), (x1, y0), (x0, y1), (x1, y1)].map(|(x, y)| m.to_lonlat(Coord::new(x, y))); - corners.iter().flatten().fold(None, |e: Option, c| Some(e.unwrap_or(Extent::EMPTY).union(&Extent::new(c.x, c.y, c.x, c.y)))) + let saved = project.extent.and_then(|[x0, y0, x1, y1]| match degrees { + // In degrees already (from Greenwich: the view of a project on another prime meridian + // is left to its data). + Some(g) if g.greenwich => Some(Extent::new(x0.max(-180.0), y0.max(-90.0), x1.min(180.0), y1.min(90.0))).filter(|e| !e.is_empty()), + Some(_) => None, + None => { + let m = s.map(); + let corners = [(x0, y0), (x1, y0), (x0, y1), (x1, y1)].map(|(x, y)| m.to_lonlat(Coord::new(x, y))); + corners + .iter() + .flatten() + .fold(None, |e: Option, c| Some(e.unwrap_or(Extent::EMPTY).union(&Extent::new(c.x, c.y, c.x, c.y)))) + } }); let target = match saved { Some(e) if !e.is_empty() && (data.is_empty() || e.intersects(&data)) => e, @@ -318,6 +335,40 @@ fn go_to_view(s: &mut Session, project: &qgis::Project) { } } +/// A coordinate system in degrees (EPSG:4326, ETRS89, NAD83, …), which QGIS draws as plate carrée. +#[derive(Clone, Copy, Debug)] +struct Geographic { + /// On WGS 84, or a datum PROJ takes as the same (no shift to WGS 84). + wgs84: bool, + /// Longitudes from Greenwich (not from another prime meridian, such as Paris's). + greenwich: bool, +} + +impl Geographic { + /// `None` for a coordinate system that isn't in degrees, or isn't known: an `EPSG:` code, + /// `OGC:CRS84`, or a PROJ string (a project QGIS saved without an authority code). + fn of(crs: &str) -> Option { + if crs.eq_ignore_ascii_case("OGC:CRS84") { + return Some(Geographic { wgs84: true, greenwich: true }); + } + let def = if crs.trim_start().starts_with('+') { + crs + } else { + let code = crs.strip_prefix("EPSG:")?.trim().parse::().ok()?; + rustgis_data::raster::crs::epsg_proj4(code)? + }; + let has = |p: &str| def.split_whitespace().any(|t| t == p); + ["+proj=longlat", "+proj=latlong", "+proj=lonlat", "+proj=latlon"].iter().any(|p| has(p)).then(|| Geographic { + wgs84: has("+datum=WGS84") || has("+datum=NAD83") || has("+towgs84=0,0,0,0,0,0,0") || has("+towgs84=0,0,0"), + // `+pm=` names a meridian or gives its offset in degrees (`0`, `0.0`: Greenwich). + greenwich: def + .split_whitespace() + .filter_map(|t| t.strip_prefix("+pm=")) + .all(|pm| pm.eq_ignore_ascii_case("greenwich") || pm.parse::().is_ok_and(|d| d == 0.0)), + }) + } +} + fn note_json(n: &qgis::Note) -> Value { json!({"layer": n.layer, "kind": n.kind, "text": n.text, "count": n.count}) } diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index e3a15e5da..609131a5b 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -623,6 +623,73 @@ fn rivers_are_labelled_along_their_lines_in_italics() { assert!(!text.contains("italic"), "{text}"); } +#[test] +fn a_qgis_project_in_degrees_opens_in_equirectangular() { + // QGIS draws a project in degrees as plate carrée, longitude and latitude as x and y. + // A PROJ string goes where QGIS keeps it for a system without an authority code. + let qgs = |crs: &str| { + let tag = if crs.starts_with('+') { "proj4" } else { "authid" }; + format!( + r#"World<{tag}>{crs} + 5451550"# + ) + }; + let mut s = world(); + for (crs, datum_note) in [ + ("EPSG:4326", false), + ("EPSG:4258", false), + ("OGC:CRS84", false), + ("EPSG:4314", true), + // Saved without an authority code: its PROJ string. + ("+proj=longlat +datum=WGS84 +no_defs", false), + ("+proj=longlat +datum=WGS84 +pm=0.0 +no_defs", false), + ("+proj=longlat +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +no_defs", true), + ] { + let out = s.run("project.importQgis", &json!({"text": qgs(crs)})).unwrap(); + assert_eq!(s.map().projection.code(), "equirectangular", "{crs}"); + let text = serde_json::to_string(&out["notes"]).unwrap(); + assert!(!text.contains("Web Mercator"), "{crs}: {text}"); + assert_eq!(text.contains("drawn in Equirectangular on WGS 84"), datum_note, "{crs}: {text}"); + // The saved view, read in degrees. + let c = s.map().to_lonlat(s.map().camera.center).unwrap(); + assert!((c.x - 10.0).abs() < 0.5 && (c.y - 47.5).abs() < 0.5, "{crs}: {c:?}"); + } + // Longitudes from Paris (in grads): drawn the same way, the view left to the data. + let out = s.run("project.importQgis", &json!({"text": qgs("EPSG:4807")})).unwrap(); + assert_eq!(s.map().projection.code(), "equirectangular"); + let text = serde_json::to_string(&out["notes"]).unwrap(); + assert!(text.contains("on WGS 84") && text.contains("from Greenwich"), "{text}"); + // A projected one keeps its own projection. + s.run("project.importQgis", &json!({"text": qgs("EPSG:2056").replace("5", "2600000")})).unwrap(); + assert_eq!(s.map().projection.code(), "EPSG:2056"); +} + +#[test] +fn a_qgis_project_with_a_geographic_proj_string_keeps_its_saved_view() { + for (crs, datum_note, saved_view) in [ + ("+proj=longlat +datum=WGS84", false, true), + ("+proj=latlong +datum=WGS84", false, true), + ("+proj=lonlat +datum=WGS84", false, true), + ("+proj=latlon +datum=WGS84", false, true), + (" +datum=WGS84 +proj=longlat +no_defs", false, true), + ("+proj=longlat +datum=potsdam", true, true), + ("+proj=longlat +datum=WGS84 +pm=paris", true, false), + ] { + let text = format!( + r#"{crs} +5451550"# + ); + let mut s = world(); + let out = s.run("project.importQgis", &json!({"text": text})).unwrap(); + assert_eq!(s.map().projection.code(), "equirectangular", "{crs}"); + let notes = serde_json::to_string(&out["notes"]).unwrap(); + assert!(!notes.contains("Web Mercator"), "{crs}: {notes}"); + assert_eq!(notes.contains("drawn in Equirectangular on WGS 84"), datum_note, "{crs}: {notes}"); + let c = s.map().to_lonlat(s.map().camera.center).unwrap(); + assert_eq!((c.x - 10.0).abs() < 0.5 && (c.y - 47.5).abs() < 0.5, saved_view, "{crs}: {c:?}"); + } +} + #[test] fn a_qgis_project_opens_as_a_map_and_a_qgis_style_styles_a_layer() { let dir = std::env::temp_dir().join(format!("rustgis-qgis-{}", std::process::id())); diff --git a/docs/user-guide/workspace.md b/docs/user-guide/workspace.md index 146d20596..054c0d295 100644 --- a/docs/user-guide/workspace.md +++ b/docs/user-guide/workspace.md @@ -27,7 +27,9 @@ Its layers come in the layer tree's order, each from its file (GeoPackage, Shape GeoTIFF and the other formats RustGIS reads, relative paths from the project's folder), a CSV of points, XYZ tiles or a WMS layer, unchecked layers and groups hidden. The map takes the project's coordinate system and its saved view (or its layers' extent, when that view shows none of them), -with no basemap and no graticule, as QGIS draws it. The import is one undo step. +with no basemap and no graticule, as QGIS draws it. A project in degrees (EPSG:4326 and other supported +longitude/latitude systems) opens in the Equirectangular projection, the plate carrée QGIS draws +it in. The import is one undo step. Each layer keeps its look: