Merge pull request #43 from stac-utils/guhidalgo/inferepsg4326fromcf

Infer epsg:4326 CRS if latitude and longitude coordinates detected

Author: ghidalgo3

.github/workflows/continuous-integration.yml

@@ -22,7 +22,7 @@ jobs:
           python-version: ${{ matrix.python-version }}
 
       - name: Install dependencies
-        run: python -m pip install .[dev]
+        run: python -m pip install -r requirements.dev.txt && python -m pip install .[dev] --no-deps
 
       - name: Run tests
         run: pytest tests -v

pyproject.toml

@@ -28,7 +28,7 @@ requires-python = ">=3.10"
 dynamic = ["version", "description"]
 
 [project.optional-dependencies]
-dev = ["pytest", "pre-commit", "shapely", "cftime", "kerchunk", "h5netcdf", "fsspec"]
+dev = ["pytest", "pre-commit", "shapely", "cftime", "kerchunk", "h5netcdf", "fsspec", "pooch"]
 
 [project.scripts]
 xstac = "xstac._generate:main"

requirements.dev.txt

@@ -0,0 +1,191 @@
+# This file was autogenerated by uv via the following command:
+#    uv pip compile -o ./requirements.dev.txt ./pyproject.toml --extra dev
+aiohappyeyeballs==2.4.6
+    # via aiohttp
+aiohttp==3.11.12
+    # via fsspec
+aiosignal==1.3.2
+    # via aiohttp
+annotated-types==0.7.0
+    # via pydantic
+asciitree==0.3.3
+    # via zarr
+async-timeout==5.0.1
+    # via aiohttp
+attrs==25.1.0
+    # via
+    #   aiohttp
+    #   jsonschema
+    #   referencing
+certifi==2025.1.31
+    # via
+    #   pyproj
+    #   requests
+cf-xarray==0.10.0
+    # via xstac (./pyproject.toml)
+cfgv==3.4.0
+    # via pre-commit
+cftime==1.6.4.post1
+    # via xstac (./pyproject.toml)
+charset-normalizer==3.4.1
+    # via requests
+click==8.1.8
+    # via dask
+cloudpickle==3.1.1
+    # via dask
+dask==2025.2.0
+    # via xstac (./pyproject.toml)
+distlib==0.3.9
+    # via virtualenv
+exceptiongroup==1.2.2
+    # via pytest
+fasteners==0.19
+    # via zarr
+filelock==3.17.0
+    # via virtualenv
+frozenlist==1.5.0
+    # via
+    #   aiohttp
+    #   aiosignal
+fsspec==2025.2.0
+    # via
+    #   xstac (./pyproject.toml)
+    #   dask
+    #   kerchunk
+h5netcdf==1.5.0
+    # via xstac (./pyproject.toml)
+h5py==3.13.0
+    # via h5netcdf
+identify==2.6.7
+    # via pre-commit
+idna==3.10
+    # via
+    #   requests
+    #   yarl
+importlib-metadata==8.6.1
+    # via dask
+iniconfig==2.0.0
+    # via pytest
+jsonschema==4.23.0
+    # via xstac (./pyproject.toml)
+jsonschema-specifications==2024.10.1
+    # via jsonschema
+kerchunk==0.2.7
+    # via xstac (./pyproject.toml)
+locket==1.0.0
+    # via partd
+multidict==6.1.0
+    # via
+    #   aiohttp
+    #   yarl
+nodeenv==1.9.1
+    # via pre-commit
+numcodecs==0.13.1
+    # via
+    #   kerchunk
+    #   zarr
+numpy==2.2.3
+    # via
+    #   xstac (./pyproject.toml)
+    #   cftime
+    #   dask
+    #   h5py
+    #   kerchunk
+    #   numcodecs
+    #   pandas
+    #   shapely
+    #   xarray
+    #   zarr
+packaging==24.2
+    # via
+    #   dask
+    #   h5netcdf
+    #   pooch
+    #   pytest
+    #   xarray
+pandas==2.2.3
+    # via
+    #   xstac (./pyproject.toml)
+    #   xarray
+partd==1.4.2
+    # via dask
+platformdirs==4.3.6
+    # via
+    #   pooch
+    #   virtualenv
+pluggy==1.5.0
+    # via pytest
+pooch==1.8.2
+    # via xstac (./pyproject.toml)
+pre-commit==4.1.0
+    # via xstac (./pyproject.toml)
+propcache==0.3.0
+    # via
+    #   aiohttp
+    #   yarl
+pydantic==2.10.6
+    # via xstac (./pyproject.toml)
+pydantic-core==2.27.2
+    # via pydantic
+pyproj==3.7.1
+    # via xstac (./pyproject.toml)
+pystac==1.12.1
+    # via xstac (./pyproject.toml)
+pytest==8.3.4
+    # via xstac (./pyproject.toml)
+python-dateutil==2.9.0.post0
+    # via
+    #   pandas
+    #   pystac
+pytz==2025.1
+    # via pandas
+pyyaml==6.0.2
+    # via
+    #   dask
+    #   pre-commit
+referencing==0.36.2
+    # via
+    #   jsonschema
+    #   jsonschema-specifications
+requests==2.32.3
+    # via pooch
+rpds-py==0.23.0
+    # via
+    #   jsonschema
+    #   referencing
+shapely==2.0.7
+    # via xstac (./pyproject.toml)
+six==1.17.0
+    # via python-dateutil
+tomli==2.2.1
+    # via pytest
+toolz==1.0.0
+    # via
+    #   dask
+    #   partd
+typing-extensions==4.12.2
+    # via
+    #   multidict
+    #   pydantic
+    #   pydantic-core
+    #   referencing
+tzdata==2025.1
+    # via pandas
+ujson==5.10.0
+    # via kerchunk
+urllib3==2.3.0
+    # via requests
+virtualenv==20.29.2
+    # via pre-commit
+xarray==2025.1.2
+    # via
+    #   xstac (./pyproject.toml)
+    #   cf-xarray
+yarl==1.18.3
+    # via aiohttp
+zarr==2.18.3
+    # via
+    #   xstac (./pyproject.toml)
+    #   kerchunk
+zipp==3.21.0
+    # via importlib-metadata

tests/test_xstac.py

@@ -165,16 +165,28 @@ def test_from_pystac_object(ds_without_spatial_dims):
 
 
 @pytest.mark.parametrize(
-    "ds",
+    ("ds", "expected"),
     [
-        xr.Dataset(coords={"epsg": xr.DataArray(32633, name="epsg")}),
-        xr.Dataset(coords={"proj:epsg": xr.DataArray(32633, name="proj:epsg")}),
-        xr.Dataset(attrs={"crs": "epsg:32633"}),
+        (xr.Dataset(coords={"epsg": xr.DataArray(32633, name="epsg")}), 32633),
+        (
+            xr.Dataset(coords={"proj:epsg": xr.DataArray(32633, name="proj:epsg")}),
+            32633,
+        ),
+        (xr.Dataset(attrs={"crs": "epsg:32633"}), 32633),
     ],
 )
-def test_maybe_infer_reference_system(ds):
+def test_maybe_infer_reference_system(ds, expected):
     result = maybe_infer_reference_system(ds, reference_system=None)
-    expected = pyproj.crs.CRS.from_epsg(32633).to_json_dict()
+    expected = pyproj.crs.CRS.from_epsg(expected).to_json_dict()
+    assert result == expected
+
+
+def test_maybe_infer_reference_system_from_cf_coordinates(ds):
+    for n, variable in ds.variables.items():
+        if "grid_mapping_name" in variable.attrs:
+            ds[n].attrs.pop("grid_mapping_name")
+    result = maybe_infer_reference_system(ds, reference_system=None)
+    expected = pyproj.crs.CRS.from_epsg(4326).to_json_dict()
     assert result == expected
 
 

xstac/_xstac.py

@@ -204,6 +204,7 @@ def maybe_infer_reference_system(ds, reference_system) -> dict:
     * proj:epsg from the coords
     * crs from the attrs
     * variables with grid_mapping_name in their attrs.
+    * CF convention detection through cf_xarray
 
     """
     if reference_system is None:
@@ -220,17 +221,21 @@ def maybe_infer_reference_system(ds, reference_system) -> dict:
                 ds.attrs["crs"]
             ).to_json_dict()
         else:
-            # try to infer it
+            # try to infer it from variables
             names = [
                 k for k, v in ds.variables.items() if "grid_mapping_name" in v.attrs
             ]
-            if not names:
-                raise ValueError("Couldn't find a reference system")
-            elif len(names) > 1:
+            if len(names) > 1:
                 raise ValueError("Too many reference systems: %s", names)
-            (name,) = names
-            crs = CRS.from_cf(ds[name].attrs)
-            reference_system = crs.to_json_dict()
+            elif len(names) == 1:
+                (name,) = names
+                crs = CRS.from_cf(ds[name].attrs)
+                reference_system = crs.to_json_dict()
+            elif {"latitude", "longitude"} <= ds.cf.coordinates.keys():
+                reference_system = pyproj.CRS.from_epsg(4326).to_json_dict()
+
+            if not reference_system:
+                raise ValueError("Couldn't find a reference system")
 
     elif reference_system is False:
         reference_system = None