Skip to content

fix(data): planetary CRS tests match equidistant cylindrical on its ellipsoid (main CI) #1588

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
giswqs wants to merge 1 commit into from
Closed

fix(data): planetary CRS tests match equidistant cylindrical on its ellipsoid (main CI) #1588

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 20 additions & 7 deletions crates/data/tests/planetary_crs.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,14 @@ fn planetographic_systems_without_east_axes_are_refused() {
let dir = tempfile::tempdir().unwrap();
let r = 3_396_190.0f64;
let metres = LONLAT.map(|(lon, lat)| (r * lon.to_radians(), r * lat.to_radians()));
// The same points projected on the planetographic ellipsoid, as RustGIS computes an
// equidistant cylindrical projection on an ellipsoid (EPSG method 1028).
let f = 1.0 / 169.894447223612;
let on_ellipsoid = rustgis_geom::EquidistantCylindrical::new(r, f * (2.0 - f), 0.0, 0.0, 0.0, 0.0, 0.0).unwrap();

@github-actions github-actions Bot Nov 8, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Quality / test coverage (confidence: medium). The new fixture metres come from rustgis_geom::EquidistantCylindrical, the same code the reader uses, so the round trip cannot catch an error in that projection. The PR body says whether reading eqc on an ellipsoid is right is still open (#1586), and these assertions now pin that behaviour. Consider one hard-coded reference value per datum, worked out independently (for example from the GN7-2 formula, or a published PROJ result), next to the computed fixtures. Left for the release quality pass if you prefer; it doesn't block the CI fix.

let ellipsoid_metres = LONLAT.map(|(lon, lat)| {
let m = on_ellipsoid.forward(Coord::new(lon, lat)).unwrap();
(m.x, m.y)
});
let with_axes = |wkt: &str, axes: &str| format!("{},{axes}]", wkt.strip_suffix(']').unwrap());
let east_north = r#"AXIS["Lon",EAST],AXIS["Lat",NORTH]"#;
let unnamed = MARS_OGRAPHIC_ESRI.replace("GCS_Mars_2015_Ographic", "GCS_Mars");
Expand Down Expand Up @@ -194,14 +202,14 @@ fn planetographic_systems_without_east_axes_are_refused() {
("esri_east", with_axes(MARS_OGRAPHIC_ESRI, east_north), LONLAT),
("wkt2_east", MARS_WEST.replace(",west,", ",east,"), LONLAT),
("ocentric_code", with_axes(&unnamed, r#"AUTHORITY["IAU_2015","49900"]"#), LONLAT),
("eqc_easting", eqc_on(&with_axes(MARS_OGRAPHIC_ESRI, east_north), "Mars_Equirectangular"), metres),
("eqc_easting", eqc_on(&with_axes(MARS_OGRAPHIC_ESRI, east_north), "Mars_Equirectangular"), ellipsoid_metres),
(
"eqc_base_easting",
with_axes(&eqc_on(MARS_OGRAPHIC_ESRI, "Mars_Equirectangular"), r#"AXIS["Easting",EAST],AXIS["Northing",NORTH]"#),
metres,
ellipsoid_metres,
),
("eqc_ocentric_code", with_axes(MARS_EQC, r#"AUTHORITY["IAU_2015","49910"]"#), metres),
("wkt2_eqc_easting", wkt2_eqc("Mars (2015) / Ographic", "east"), metres),
("wkt2_eqc_easting", wkt2_eqc("Mars (2015) / Ographic", "east"), ellipsoid_metres),
];
for (name, prj, native) in &read {
let path = shapefile(dir.path(), name, *native, prj);
Expand Down Expand Up @@ -239,7 +247,12 @@ fn planetographic_systems_without_east_axes_are_refused() {
let src = source::open_path(&path, None).unwrap();
assert_eq!(said(&src.info().notes), 1, "{:?}", src.info().notes);
let a = 6_378_137.0f64;
let earth_metres = LONLAT.map(|(lon, lat)| (a * lon.to_radians(), a * lat.to_radians()));
let f = 1.0 / 298.257223563;

@github-actions github-actions Bot Nov 8, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit (confidence: high, left for the release quality pass). This let f shadows the earlier let f in the same function, and the two blocks that turn lon/lat into ellipsoidal metres are near-copies. A small helper, such as fn ellipsoid_metres(a: f64, inv_f: f64, lonlat: [(f64, f64); 2]) -> [(f64, f64); 2], would remove both.

let on_wgs84 = rustgis_geom::EquidistantCylindrical::new(a, f * (2.0 - f), 0.0, 0.0, 0.0, 0.0, 0.0).unwrap();
let earth_metres = LONLAT.map(|(lon, lat)| {
let m = on_wgs84.forward(Coord::new(lon, lat)).unwrap();
(m.x, m.y)
});
let earth_eqc = eqc_on(&west, "World_Equidistant_Cylindrical");
assert!(earth_eqc.contains("WGS_1984") && earth_eqc.contains("WEST"), "{earth_eqc}");
let path = shapefile(dir.path(), "earth_eqc_west", earth_metres, &earth_eqc);
Expand Down Expand Up @@ -744,13 +757,13 @@ fn a_mars_projection_in_a_geotiff_is_said_or_refused() {
// On the Earth, such a citation changes nothing.
let earth = cited_on(6_378_137.0, 3073, "Ographic test grid|").unwrap();
assert!(!earth.crs.is_planetary() && earth.lost.is_empty(), "{:?}", earth.lost);
// A transformation RustGIS doesn't read (23, the New Zealand Map Grid's): refused on Mars, even
// A transformation RustGIS doesn't read (26, the New Zealand Map Grid's): refused on Mars, even
// where the numbers look like degrees, as they are metres.
let mut unread = eqc.clone();
let at = unread.chunks(4).position(|k| k[0] == 3075).unwrap() * 4 + 3;
unread[at] = 23;
unread[at] = 26;
let e = raster::read_bytes("unread.tif", &geotiff(-180.0, &unread, &[3_396_190.0], "")).unwrap_err().to_string();
assert!(e.contains(ANOTHER_BODY) && e.contains("transformation 23"), "{e}");
assert!(e.contains(ANOTHER_BODY) && e.contains("transformation 26"), "{e}");
// On the Earth it is read as longitude/latitude where the numbers fit, as before.
let earth = raster::read_bytes("earth.tif", &geotiff(-180.0, &unread, &[6_378_137.0], "")).unwrap();
assert_eq!(earth.crs, Crs::LonLat);
Expand Down
Loading