You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by pa...@apache.org on 2026/06/19 19:28:53 UTC
(sedona-db) branch main updated: feat(rust/sedona-raster-zarr): resolve CRS from CF/rioxarray conventions (#985)
This is an automated email from the ASF dual-hosted git repository.
paleolimbot pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/sedona-db.git
The following commit(s) were added to refs/heads/main by this push:
new e25f081e7 feat(rust/sedona-raster-zarr): resolve CRS from CF/rioxarray conventions (#985)
e25f081e7 is described below
commit e25f081e739bacd42aabf8b9205cdc6e070359ff
Author: James Willis <ja...@users.noreply.github.com>
AuthorDate: Fri Jun 19 12:28:48 2026 -0700
feat(rust/sedona-raster-zarr): resolve CRS from CF/rioxarray conventions (#985)
---
rust/sedona-raster-zarr/src/geozarr.rs | 88 +++++++++++++++-
rust/sedona-raster-zarr/src/loader.rs | 133 +++++++++++++++++++++++-
rust/sedona-raster-zarr/tests/zarr_roundtrip.rs | 125 ++++++++++++++++++++++
3 files changed, 343 insertions(+), 3 deletions(-)
diff --git a/rust/sedona-raster-zarr/src/geozarr.rs b/rust/sedona-raster-zarr/src/geozarr.rs
index 1bc3f13d5..b3acd15d8 100644
--- a/rust/sedona-raster-zarr/src/geozarr.rs
+++ b/rust/sedona-raster-zarr/src/geozarr.rs
@@ -31,8 +31,12 @@ use arrow_schema::ArrowError;
#[derive(Debug, Clone, PartialEq)]
pub struct GroupGeoMetadata {
/// CRS string in PROJ or WKT format (whichever the group declared).
- /// `None` if no `proj:wkt2` / `proj:projjson` / `proj:code` (or the
- /// legacy `proj:epsg`) attribute is present on the group.
+ /// `None` if the group declares no CRS through either the GeoZarr
+ /// `proj:*` convention (`proj:wkt2` / `proj:projjson` / `proj:code`, or
+ /// the legacy `proj:epsg`) or the CF / rioxarray convention (`crs_wkt` /
+ /// `spatial_ref`). The CF convention also commonly puts these attributes
+ /// on a separate grid-mapping variable rather than the group; that case
+ /// is resolved in the loader (see `crs_from_cf_attributes`).
pub crs: Option<String>,
/// Affine transform, stored in GDAL GeoTransform order:
/// `[origin_x, scale_x, skew_x, origin_y, skew_y, scale_y]`. The source
@@ -97,6 +101,34 @@ fn parse_crs(
})?;
return Ok(Some(format!("EPSG:{code}")));
}
+ // No GeoZarr `proj:*` key; fall back to the CF / rioxarray convention,
+ // which much of the public Zarr corpus (xarray + rioxarray writers) uses.
+ crs_from_cf_attributes(attrs)
+}
+
+/// Read a CRS from the CF / rioxarray attribute convention.
+///
+/// rioxarray writes the CRS as `crs_wkt` (WKT2) and, redundantly, as a
+/// `spatial_ref` attribute holding the same WKT; some writers (and group
+/// roots) put an `authority:code` string in `spatial_ref` instead. `crs_wkt`
+/// wins over `spatial_ref` when both are present.
+///
+/// These attributes live either on the group (handled by `parse_crs`) or on
+/// a separate grid-mapping variable (handled by the loader), so this is
+/// factored out for both callers.
+///
+/// PROJ.4 strings (`proj4` / `proj4_params`) are intentionally not read: the
+/// CRS layer accepts authority codes, WKT and PROJJSON, but not PROJ.4, so
+/// surfacing one would only produce a downstream parse error.
+pub fn crs_from_cf_attributes(
+ attrs: &serde_json::Map<String, serde_json::Value>,
+) -> Result<Option<String>, ArrowError> {
+ if let Some(v) = attrs.get("crs_wkt") {
+ return Ok(Some(json_value_to_string(v, "crs_wkt")?));
+ }
+ if let Some(v) = attrs.get("spatial_ref") {
+ return Ok(Some(json_value_to_string(v, "spatial_ref")?));
+ }
Ok(None)
}
@@ -235,6 +267,58 @@ mod tests {
assert!(crs.contains("GeographicCRS"));
}
+ #[test]
+ fn cf_crs_wkt_parses_to_wkt_string() {
+ let g = GroupGeoMetadata::from_attributes(&map(json!({
+ "crs_wkt": "PROJCS[\"unknown\",GEOGCS[\"unknown\", ...]]"
+ })))
+ .unwrap();
+ assert!(g.crs.as_deref().unwrap().starts_with("PROJCS"));
+ }
+
+ #[test]
+ fn cf_spatial_ref_authority_code_parses() {
+ // rioxarray / CarbonPlan group roots put an authority code here.
+ let g = GroupGeoMetadata::from_attributes(&map(json!({
+ "spatial_ref": "EPSG:3031",
+ "proj4_params": "+proj=stere +lat_0=-90"
+ })))
+ .unwrap();
+ assert_eq!(g.crs.as_deref(), Some("EPSG:3031"));
+ }
+
+ #[test]
+ fn cf_crs_wkt_takes_precedence_over_spatial_ref() {
+ let g = GroupGeoMetadata::from_attributes(&map(json!({
+ "crs_wkt": "PROJCS[\"a\"]",
+ "spatial_ref": "EPSG:3031"
+ })))
+ .unwrap();
+ assert_eq!(g.crs.as_deref(), Some("PROJCS[\"a\"]"));
+ }
+
+ #[test]
+ fn geozarr_proj_code_takes_precedence_over_cf() {
+ // A group declaring both conventions resolves via GeoZarr `proj:*`.
+ let g = GroupGeoMetadata::from_attributes(&map(json!({
+ "proj:code": "EPSG:4326",
+ "crs_wkt": "PROJCS[\"other\"]"
+ })))
+ .unwrap();
+ assert_eq!(g.crs.as_deref(), Some("EPSG:4326"));
+ }
+
+ #[test]
+ fn cf_proj4_alone_is_not_read() {
+ // PROJ.4 is not accepted by the CRS layer, so a group that declares
+ // only `proj4` stays CRS-less rather than surfacing an unusable string.
+ let g = GroupGeoMetadata::from_attributes(&map(json!({
+ "proj4": "+proj=stere +lat_0=-90 +lat_ts=-71"
+ })))
+ .unwrap();
+ assert!(g.crs.is_none());
+ }
+
#[test]
fn transform_affine_reorders_to_gdal() {
// Affine [a, b, c, d, e, f] = [1, 0, 100, 0, -1, 200]: north-up,
diff --git a/rust/sedona-raster-zarr/src/loader.rs b/rust/sedona-raster-zarr/src/loader.rs
index f2514b882..622f26a23 100644
--- a/rust/sedona-raster-zarr/src/loader.rs
+++ b/rust/sedona-raster-zarr/src/loader.rs
@@ -42,7 +42,7 @@ use zarrs::storage::{AsyncReadableListableStorage, AsyncReadableListableStorageT
use crate::coords;
use crate::dtype::zarr_to_band_data_type;
-use crate::geozarr::GroupGeoMetadata;
+use crate::geozarr::{crs_from_cf_attributes, GroupGeoMetadata};
use crate::source_uri::build_chunk_anchor;
/// Streaming reader over the chunk grid of a Zarr group.
@@ -335,6 +335,16 @@ async fn open_and_validate(
}
};
+ // CF / rioxarray groups declare their CRS on a separate scalar
+ // grid-mapping variable (commonly `spatial_ref`) rather than in the group
+ // attributes. When the group attributes didn't already yield a CRS,
+ // consult that variable before we drop the open array handles.
+ if geo.crs.is_none() {
+ if let Some(crs) = crs_from_grid_mapping_variable(&storage, &arrays).await {
+ geo.crs = Some(crs);
+ }
+ }
+
let array_infos = collect_array_infos(arrays)?;
validate_group_constraints(&array_infos)?;
@@ -480,6 +490,83 @@ async fn open_named_arrays(
Ok(out)
}
+/// Resolve a CRS from a CF / rioxarray grid-mapping variable.
+///
+/// rioxarray writes the CRS to a scalar variable (conventionally
+/// `spatial_ref`, sometimes `crs`) carrying a `crs_wkt` / `spatial_ref`
+/// attribute, linked from each data variable via the CF `grid_mapping`
+/// attribute — or, as rioxarray does, listed in the data variable's
+/// `coordinates` attribute. That variable is a 0-D scalar, so it never
+/// appears among the rank ≥ 2 data arrays and must be opened on its own.
+///
+/// Candidate names are gathered in priority order — `grid_mapping` targets
+/// first, then `coordinates` tokens, then the conventional `spatial_ref` /
+/// `crs` — and the first that opens and yields a CRS wins. A candidate that
+/// is missing, unreadable, or carries no CRS attribute is skipped, so a group
+/// without a grid-mapping variable simply stays CRS-less.
+async fn crs_from_grid_mapping_variable(
+ storage: &AsyncReadableListableStorage,
+ data_arrays: &[Array<dyn AsyncReadableListableStorageTraits>],
+) -> Option<String> {
+ let candidates = cf_crs_candidate_names(data_arrays.iter().map(|a| a.attributes()));
+ for name in candidates {
+ let path = format!("/{}", name.trim_start_matches('/'));
+ let array = match Array::async_open(storage.clone(), &path).await {
+ Ok(array) => array,
+ // An absent candidate and a transient/permission failure look
+ // alike here; CRS is optional so we keep trying, but leave a
+ // breadcrumb so a genuinely broken store stays traceable.
+ Err(e) => {
+ log::debug!("Zarr grid-mapping candidate {name:?}: open failed: {e}");
+ continue;
+ }
+ };
+ match crs_from_cf_attributes(array.attributes()) {
+ Ok(Some(crs)) => return Some(crs),
+ Ok(None) => continue,
+ // A malformed CRS attribute on a candidate variable shouldn't take
+ // the whole read offline; log and keep looking.
+ Err(e) => log::warn!("Zarr grid-mapping variable {name:?}: ignoring CRS: {e}"),
+ }
+ }
+ None
+}
+
+/// Gather candidate grid-mapping variable names from the data arrays'
+/// attributes, in the order they should be tried: every `grid_mapping`
+/// target, then every `coordinates` token, then the conventional
+/// `spatial_ref` / `crs`. De-duplicated, first occurrence wins.
+///
+/// CF 1.7+ allows an extended `grid_mapping` form, `"<crs var>: <coord>
+/// <coord> ..."`, where the grid-mapping variable name carries a trailing
+/// colon; we strip it so that name resolves. The colon-less coordinate names
+/// in that form are tried too — harmless, they carry no CRS attribute.
+fn cf_crs_candidate_names<'a>(
+ array_attrs: impl Iterator<Item = &'a serde_json::Map<String, serde_json::Value>>,
+) -> Vec<String> {
+ let mut names: Vec<String> = Vec::new();
+ let mut push = |s: &str| {
+ for tok in s.split_whitespace() {
+ let tok = tok.trim_end_matches(':');
+ if !tok.is_empty() && !names.iter().any(|n| n == tok) {
+ names.push(tok.to_string());
+ }
+ }
+ };
+ let attrs: Vec<_> = array_attrs.collect();
+ for key in ["grid_mapping", "coordinates"] {
+ for a in &attrs {
+ if let Some(s) = a.get(key).and_then(|v| v.as_str()) {
+ push(s);
+ }
+ }
+ }
+ for conventional in ["spatial_ref", "crs"] {
+ push(conventional);
+ }
+ names
+}
+
/// Discover the group's direct child arrays.
///
/// Discovery goes through zarrs's [`Group::async_children`], which uses
@@ -949,6 +1036,50 @@ mod tests {
assert_eq!(idx, vec![0, 1]);
}
+ fn attrs(json: serde_json::Value) -> serde_json::Map<String, serde_json::Value> {
+ json.as_object().unwrap().clone()
+ }
+
+ #[test]
+ fn cf_crs_candidate_names_prioritizes_grid_mapping_then_coordinates() {
+ // grid_mapping target wins, coordinates tokens follow (real coord
+ // vars like `time` get tried too — harmless, they carry no CRS), and
+ // the conventional names are appended last.
+ let a = attrs(serde_json::json!({
+ "grid_mapping": "crs",
+ "coordinates": "spatial_ref time"
+ }));
+ let got = cf_crs_candidate_names(std::iter::once(&a));
+ assert_eq!(got, vec!["crs", "spatial_ref", "time"]);
+ }
+
+ #[test]
+ fn cf_crs_candidate_names_falls_back_to_conventional() {
+ // No grid_mapping / coordinates anywhere: try the conventional names.
+ let a = attrs(serde_json::json!({"units": "K"}));
+ let got = cf_crs_candidate_names(std::iter::once(&a));
+ assert_eq!(got, vec!["spatial_ref", "crs"]);
+ }
+
+ #[test]
+ fn cf_crs_candidate_names_strips_extended_grid_mapping_colon() {
+ // CF 1.7+ extended form: the colon-suffixed token is the grid-mapping
+ // variable; the rest are the coordinates it applies to.
+ let a = attrs(serde_json::json!({"grid_mapping": "crs: lon lat"}));
+ let got = cf_crs_candidate_names(std::iter::once(&a));
+ // `crs` (colon stripped) first, then the coord tokens, then the
+ // conventional `spatial_ref` (`crs` already present, deduped).
+ assert_eq!(got, vec!["crs", "lon", "lat", "spatial_ref"]);
+ }
+
+ #[test]
+ fn cf_crs_candidate_names_dedupes_across_arrays() {
+ let a = attrs(serde_json::json!({"coordinates": "spatial_ref time"}));
+ let b = attrs(serde_json::json!({"coordinates": "spatial_ref lat lon"}));
+ let got = cf_crs_candidate_names([&a, &b].into_iter());
+ assert_eq!(got, vec!["spatial_ref", "time", "lat", "lon", "crs"]);
+ }
+
#[test]
fn resolve_spatial_dim_indices_explicit_missing_errors() {
let names = vec!["a".into(), "b".into()];
diff --git a/rust/sedona-raster-zarr/tests/zarr_roundtrip.rs b/rust/sedona-raster-zarr/tests/zarr_roundtrip.rs
index fcd0b6fe4..217628324 100644
--- a/rust/sedona-raster-zarr/tests/zarr_roundtrip.rs
+++ b/rust/sedona-raster-zarr/tests/zarr_roundtrip.rs
@@ -541,6 +541,131 @@ async fn irregular_coordinate_arrays_fall_back_to_identity() {
assert_eq!(raster.crs(), None);
}
+/// A CF / rioxarray group whose CRS lives on a separate scalar grid-mapping
+/// variable rather than in the group attributes. `temperature` has generic
+/// `y`/`x` dims (so no CRS is inferred from the names) and links to the
+/// `spatial_ref` variable via `link_attr` (`"coordinates"` as rioxarray
+/// writes it, or the CF `"grid_mapping"` attribute). The scalar `spatial_ref`
+/// variable carries the WKT in `crs_wkt` when `crs_wkt` is `Some` — `None`
+/// builds the variable with no CRS attribute, exercising the
+/// stays-CRS-less path. Regularly spaced coordinate arrays supply the
+/// transform either way.
+fn build_grid_mapping_fixture(crs_wkt: Option<&str>, link_attr: &str) -> TempDir {
+ let tmp = TempDir::new().unwrap();
+ let store = Arc::new(FilesystemStore::new(tmp.path()).unwrap());
+
+ GroupBuilder::new()
+ .build(store.clone(), "/")
+ .unwrap()
+ .store_metadata()
+ .unwrap();
+
+ // 1-D coordinate arrays (regular spacing -> derivable transform).
+ for (name, vals) in [("y", &[20.0, 19.0][..]), ("x", &[10.0, 11.0, 12.0][..])] {
+ let n = vals.len() as u64;
+ let arr = ArrayBuilder::new(vec![n], vec![n], data_type::float64(), 0.0f64)
+ .dimension_names(Some([name]))
+ .build(store.clone(), &format!("/{name}"))
+ .unwrap();
+ arr.store_metadata().unwrap();
+ arr.store_chunk(&[0], vals).unwrap();
+ }
+
+ // Scalar grid-mapping variable, with the CRS attribute only when asked.
+ let mut crs_attrs = serde_json::Map::new();
+ if let Some(wkt) = crs_wkt {
+ crs_attrs.insert("crs_wkt".into(), serde_json::json!(wkt));
+ }
+ let crs_var = ArrayBuilder::new(
+ Vec::<u64>::new(),
+ Vec::<u64>::new(),
+ data_type::int64(),
+ 0i64,
+ )
+ .attributes(crs_attrs)
+ .build(store.clone(), "/spatial_ref")
+ .unwrap();
+ crs_var.store_metadata().unwrap();
+
+ // 2-D data array linking to the CRS variable via `link_attr`.
+ let mut data_attrs = serde_json::Map::new();
+ data_attrs.insert(link_attr.into(), serde_json::json!("spatial_ref"));
+ let data = ArrayBuilder::new(vec![2u64, 3], vec![2u64, 3], data_type::uint8(), 0u8)
+ .dimension_names(Some(["y", "x"]))
+ .attributes(data_attrs)
+ .build(store.clone(), "/temperature")
+ .unwrap();
+ data.store_metadata().unwrap();
+ data.store_chunk(&[0, 0], vec![0u8; 6]).unwrap();
+
+ tmp
+}
+
+const POLAR_STEREO_WKT: &str = "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"WGS_1984\",\
+ SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],\
+ UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],\
+ PARAMETER[\"latitude_of_origin\",-71],UNIT[\"metre\",1]]";
+
+#[tokio::test]
+async fn resolves_crs_from_cf_grid_mapping_variable() {
+ // CRS is on the `spatial_ref` variable, not the group; the loader must
+ // open it and lift the WKT even though `y`/`x` imply no CRS themselves.
+ // rioxarray links it via the `coordinates` attribute.
+ let tmp = build_grid_mapping_fixture(Some(POLAR_STEREO_WKT), "coordinates");
+ let uri = format!("file://{}", tmp.path().display());
+ let arr = read_all(&uri, None).await.unwrap();
+
+ let raster = RasterStructArray::try_new(&arr).unwrap().get(0).unwrap();
+ // Transform still derived from the regular coordinate arrays.
+ assert_eq!(
+ raster.transform().to_vec(),
+ vec![9.5, 1.0, 0.0, 20.5, 0.0, -1.0]
+ );
+ assert_eq!(raster.crs(), Some(POLAR_STEREO_WKT));
+}
+
+#[tokio::test]
+async fn resolves_crs_via_cf_grid_mapping_attribute() {
+ // The CF-standard linkage: the data variable's `grid_mapping` attribute
+ // names the CRS variable directly (vs. rioxarray's `coordinates` list).
+ let tmp = build_grid_mapping_fixture(Some(POLAR_STEREO_WKT), "grid_mapping");
+ let uri = format!("file://{}", tmp.path().display());
+ let arr = read_all(&uri, None).await.unwrap();
+
+ let raster = RasterStructArray::try_new(&arr).unwrap().get(0).unwrap();
+ assert_eq!(raster.crs(), Some(POLAR_STEREO_WKT));
+}
+
+#[tokio::test]
+async fn grid_mapping_variable_without_crs_attribute_stays_crs_less() {
+ // The `spatial_ref` variable exists and opens, but carries no CRS
+ // attribute — the group must read back CRS-less rather than erroring.
+ let tmp = build_grid_mapping_fixture(None, "coordinates");
+ let uri = format!("file://{}", tmp.path().display());
+ let arr = read_all(&uri, None).await.unwrap();
+
+ let raster = RasterStructArray::try_new(&arr).unwrap().get(0).unwrap();
+ // Transform still derived from the coordinate arrays; just no CRS.
+ assert_eq!(
+ raster.transform().to_vec(),
+ vec![9.5, 1.0, 0.0, 20.5, 0.0, -1.0]
+ );
+ assert_eq!(raster.crs(), None);
+}
+
+#[tokio::test]
+async fn resolves_crs_from_grid_mapping_variable_with_arrays_filter() {
+ // The grid-mapping variable isn't in the filter, but the loader opens it
+ // by name — mirroring how coordinate arrays are resolved.
+ let tmp = build_grid_mapping_fixture(Some(POLAR_STEREO_WKT), "coordinates");
+ let uri = format!("file://{}", tmp.path().display());
+ let arrays = ["temperature".to_string()];
+ let arr = read_all(&uri, Some(&arrays)).await.unwrap();
+
+ let raster = RasterStructArray::try_new(&arr).unwrap().get(0).unwrap();
+ assert_eq!(raster.crs(), Some(POLAR_STEREO_WKT));
+}
+
#[tokio::test]
async fn derives_transform_with_explicit_arrays_filter() {
// The coordinate arrays are not in the `arrays` filter, but the reader