cmems_crs_zarr.md

In these CMEMs ZARR there's evidence of the dummy "crs" var used via "grid_mapping" in NetCDF, but there's no such var.

https://data.marine.copernicus.eu/product/SEALEVEL_GLO_PHY_L4_MY_008_047/description

curl -s "https://s3.waw3-1.cloudferro.com/mdl-arco-time-045/arco/SEALEVEL_GLO_PHY_L4_MY_008_047/cmems_obs-sl_glo_phy-ssh_my_allsat-demo-l4-duacs-0.125deg_P1D-i_202511/timeChunked.zarr/.zmetadata" | grep -i "grid_mapping\|crs\|spatial_ref"
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",
      "grid_mapping": "crs",


curl -s "https://s3.waw3-1.cloudferro.com/mdl-arco-time-045/arco/SEALEVEL_GLO_PHY_L4_MY_008_047/cmems_obs-sl_glo_phy-ssh_my_allsat-demo-l4-duacs-0.125deg_P1D-i_202511/timeChunked.zarr/.zmetadata" | grep -i "\\.zarray"
    "adt/.zarray": {
    "err_sla/.zarray": {
    "err_ugosa/.zarray": {
    "err_vgosa/.zarray": {
    "flag_ice/.zarray": {
    "lat_bnds/.zarray": {
    "latitude/.zarray": {
    "lon_bnds/.zarray": {
    "longitude/.zarray": {
    "nv/.zarray": {
    "sla/.zarray": {
    "time/.zarray": {
    "tpa_correction/.zarray": {
    "ugos/.zarray": {
    "ugosa/.zarray": {
    "vgos/.zarray": {
    "vgosa/.zarray": {

it's funny because I very clearly remember earlier products by AVISO that shipped in transposed Mercator grid but with rectilinear long lat arrays (ewk), there was a brief happy time when CF provided crs and this stuff just worked from the NetCDFs.

here is a candidate fix-up function

import zarr
import json
from pathlib import Path

def fix_zarr_crs(zarr_path, crs="EPSG:4326"):
    """Add missing CRS variable to a Zarr store (in-place).
    
    Parameters
    ----------
    zarr_path : str or Path
        Path to the Zarr store
    crs : str
        CRS identifier (e.g., "EPSG:4326", "EPSG:3857")
    """
    import pyproj
    
    zarr_path = Path(zarr_path)
    
    # Get CRS info from pyproj
    crs_obj = pyproj.CRS(crs)
    
    # Build attributes
    crs_attrs = {
        "_ARRAY_DIMENSIONS": [],
        "crs_wkt": crs_obj.to_wkt(),
        "spatial_ref": crs_obj.to_wkt(),
    }
    
    # Add grid_mapping_name and ellipsoid params for geographic CRS
    if crs_obj.is_geographic:
        crs_attrs["grid_mapping_name"] = "latitude_longitude"
        ellipsoid = crs_obj.ellipsoid
        if ellipsoid:
            crs_attrs["semi_major_axis"] = ellipsoid.semi_major_metre
            crs_attrs["inverse_flattening"] = ellipsoid.inverse_flattening
        prime_meridian = crs_obj.prime_meridian
        if prime_meridian:
            crs_attrs["longitude_of_prime_meridian"] = prime_meridian.longitude
    
    # For projected CRS, extract projection params
    elif crs_obj.is_projected:
        cf = crs_obj.to_cf()
        crs_attrs.update(cf)  # to_cf() gives CF-compatible dict
    
    # Write crs/.zarray (scalar int)
    crs_zarray = {
        "chunks": [],
        "compressor": None,
        "dtype": "<i4",
        "fill_value": 0,
        "filters": None,
        "order": "C",
        "shape": [],
        "zarr_format": 2
    }
    
    crs_dir = zarr_path / "crs"
    crs_dir.mkdir(exist_ok=True)
    
    (crs_dir / ".zarray").write_text(json.dumps(crs_zarray))
    (crs_dir / ".zattrs").write_text(json.dumps(crs_attrs))
    (crs_dir / "0").write_bytes(b'\x00\x00\x00\x00')
    
    # Update .zmetadata if it exists
    zmeta_path = zarr_path / ".zmetadata"
    if zmeta_path.exists():
        meta = json.loads(zmeta_path.read_text())
        meta["metadata"]["crs/.zarray"] = crs_zarray
        meta["metadata"]["crs/.zattrs"] = crs_attrs
        zmeta_path.write_text(json.dumps(meta, indent=2))
    
    print(f"Added CRS ({crs}) to {zarr_path}")


# Quick test
if __name__ == "__main__":
    import pyproj
    
    # Show what pyproj gives us
    for epsg in ["EPSG:4326", "EPSG:3857", "EPSG:6931"]:
        crs = pyproj.CRS(epsg)
        print(f"\n{epsg}:")
        print(f"  is_geographic: {crs.is_geographic}")
        print(f"  is_projected: {crs.is_projected}")
        if crs.ellipsoid:
            print(f"  semi_major: {crs.ellipsoid.semi_major_metre}")
            print(f"  inv_flat: {crs.ellipsoid.inverse_flattening}")
        print(f"  to_cf(): {list(crs.to_cf().keys())}")