Repository: opendatacube/odc-geo
Pull request #179 (2024-09-24 04:26:26)

Expose odc.{crs_coord,grid_mapping}
- access CRS coordinate variable/name

These seem useful



Patch:
diff --git a/README.rst b/README.rst
index 4d5637c..ac1730a 100644
--- a/README.rst
+++ b/README.rst
@@ -1,7 +1,7 @@
 odc.geo
 #######
 
-|Documentation Status| |Test Status| |Test Coverage| |Conda Version| |Discord|
+|Documentation Status| |Test Status| |Test Coverage| |Conda Version|
 
 This library combines geometry shape classes from shapely_ with CRS from pyproj_ to provide
 projection aware ``Geometry``. It exposes all the functionality provided by ``shapely`` modules, but
@@ -49,10 +49,6 @@ For details and motivation see `ODC-EP-06`_ enhancement proposal.
    :target: https://anaconda.org/conda-forge/odc-geo
    :alt:  Conda Version
 
-.. |Discord| image:: https://img.shields.io/discord/1212501566326571070?label=Discord&logo=discord&logoColor=white&color=7289DA
-   :target: https://discord.gg/4hhBQVas5U
-   :alt: Join Discord for support
-
 .. |Intro Image| image:: docs/_static/intro.svg
    :alt: Map with GeoBoxes
    :width: 256
diff --git a/docs/api-geobox.rst b/docs/api-geobox.rst
index e7ac554..446376f 100644
--- a/docs/api-geobox.rst
+++ b/docs/api-geobox.rst
@@ -8,7 +8,6 @@ GeoBox
    GeoBox
    GeoBox.from_geopolygon
    GeoBox.from_bbox
-   GeoBox.from_rio
    GeoBox.affine
    GeoBox.anchor
    GeoBox.boundary
diff --git a/docs/rtd-requirements.txt b/docs/rtd-requirements.txt
index 2ffa258..4077512 100644
--- a/docs/rtd-requirements.txt
+++ b/docs/rtd-requirements.txt
@@ -19,5 +19,5 @@ sphinx
 sphinx_rtd_theme
 sphinx-autodoc-typehints
 jupyter_sphinx
-geopandas<1
+geopandas
 ipykernel
diff --git a/odc/geo/_blocks.py b/odc/geo/_blocks.py
index 84b9580..cc78b0b 100644
--- a/odc/geo/_blocks.py
+++ b/odc/geo/_blocks.py
@@ -13,15 +13,9 @@ from .roi import VariableSizedTiles, roi_intersect3, roi_normalise, roi_shape
 from .types import Chunks2d
 
 
-def _find_common_type(array_types, scalar_type=None):
-    if scalar_type is None:
-        return np.result_type(*array_types)
-
-    if np.issubdtype(scalar_type, np.floating):
-        array_types = array_types + [0.0]
-    elif np.issubdtype(scalar_type, np.complexfloating):
-        array_types = array_types + [0j]
-    return np.result_type(*array_types)
+def _find_common_type(array_types, scalar_types):
+    # TODO: don't use find_common_type as it's being removed from numpy
+    return np.find_common_type(array_types, scalar_types)
 
 
 class BlockAssembler:
@@ -39,7 +33,7 @@ class BlockAssembler:
         self._dtype = (
             np.dtype("float32")
             if len(blocks) == 0
-            else _find_common_type([b.dtype for b in blocks.values()])
+            else _find_common_type([b.dtype for b in blocks.values()], [])
         )
         self._axis = axis
         self._blocks = blocks
@@ -132,7 +126,7 @@ class BlockAssembler:
             dtype = self._dtype
             if fill_value is not None:
                 # possibly upgrade to float based on fill_value
-                dtype = _find_common_type([dtype], np.min_scalar_type(fill_value))
+                dtype = _find_common_type([dtype], [np.min_scalar_type(fill_value)])
         else:
             dtype = np.dtype(dtype)
 
diff --git a/odc/geo/_version.py b/odc/geo/_version.py
index 84cf5b1..56038e9 100644
--- a/odc/geo/_version.py
+++ b/odc/geo/_version.py
@@ -1,3 +1,3 @@
 """version information only."""
 
-__version__ = "0.4.8"
+__version__ = "0.4.7rc2"
diff --git a/odc/geo/_xr_interop.py b/odc/geo/_xr_interop.py
index 2ca3cb0..6404d7b 100644
--- a/odc/geo/_xr_interop.py
+++ b/odc/geo/_xr_interop.py
@@ -5,11 +5,9 @@
 """
 Add ``.odc.`` extension to :py:class:`xarray.Dataset` and :class:`xarray.DataArray`.
 """
-
 from __future__ import annotations
 
 import functools
-import json
 import math
 import warnings
 from dataclasses import dataclass
@@ -211,8 +209,7 @@ def _mk_crs_coord(
     crs_wkt = cf.get("crs_wkt", None) or crs.wkt
 
     if gcps is not None:
-        # Store as string
-        cf["gcps"] = json.dumps(_gcps_to_json(gcps))
+        cf["gcps"] = _gcps_to_json(gcps)
 
     if transform is not None:
         cf["GeoTransform"] = _render_geo_transform(transform, precision=24)
@@ -526,24 +523,22 @@ def _extract_gcps(crs_coord: xarray.DataArray) -> Optional[GCPMapping]:
         return None
     crs = _extract_crs(crs_coord)
     try:
-        if isinstance(gcps, str):
-            gcps = json.loads(gcps)
         wld = Geometry(gcps, crs=crs)
         pix = [
             xy_(f["properties"]["col"], f["properties"]["row"])
             for f in gcps["features"]
         ]
         return GCPMapping(pix, wld)
-    except (IndexError, KeyError, ValueError, json.JSONDecodeError):
+    except (IndexError, KeyError, ValueError):
         return None
 
 
 def _extract_geo_transform(crs_coord: xarray.DataArray) -> Optional[Affine]:
-    geo_transform_parts = crs_coord.attrs.get("GeoTransform", "").split(" ")
-    if len(geo_transform_parts) != 6:
+    geo_transfrom_parts = crs_coord.attrs.get("GeoTransform", "").split(" ")
+    if len(geo_transfrom_parts) != 6:
         return None
     try:
-        c, a, b, f, d, e = map(float, geo_transform_parts)
+        c, a, b, f, d, e = map(float, geo_transfrom_parts)
     except ValueError:
         return None
 
diff --git a/odc/geo/crs.py b/odc/geo/crs.py
index 8d70a76..d037ffb 100644
--- a/odc/geo/crs.py
+++ b/odc/geo/crs.py
@@ -26,6 +26,8 @@ from pyproj.transformer import Transformer
 
 from .types import XY, Unset
 
+SomeCRS = Union[str, int, "CRS", _CRS, Dict[str, Any]]
+MaybeCRS = Union[SomeCRS, Unset, None]
 EPSG_UNSET = 0
 
 
@@ -393,9 +395,6 @@ class CRSMismatchError(ValueError):
     """
 
 
-SomeCRS = Union[str, int, CRS, _CRS, Dict[str, Any]]
-MaybeCRS = Union[SomeCRS, Unset, None]
-
 # fmt: off
 @overload
 def norm_crs(crs: SomeCRS) -> CRS: ...
diff --git a/odc/geo/data/__init__.py b/odc/geo/data/__init__.py
index c9433b4..2d91b7d 100644
--- a/odc/geo/data/__init__.py
+++ b/odc/geo/data/__init__.py
@@ -1,7 +1,5 @@
 import json
 import lzma
-import threading
-
 from functools import lru_cache
 from pathlib import Path
 from typing import Any, Dict, Optional, Tuple
@@ -48,54 +46,19 @@ def gbox_css() -> str:
         return src.read()
 
 
-class _CachedGeoDataFrame:
-    _instance = None
-
-    # Override in sub-classes
-    _lock = threading.Lock()
-    _data_url = ""
-
-    def __init__(self):
-        # Thread safe class-cached dataload
-        if self._instance is None:
-            with self._lock:
-                if self._instance is None:
-                    self.__class__._instance = self._load_from_url()
-
-    def _load_from_url(self):
-        # pylint: disable=import-outside-toplevel
-        import geopandas as gpd
-
-        with catch_warnings():
-            filterwarnings("ignore", category=FutureWarning)
-            df = gpd.read_file(self._data_url)
-        return df
-
-
-class Countries(_CachedGeoDataFrame):
-    """
-    Cache-wrapper around the Natural Earth low-res countries geodataset.
-    """
-
-    _lock = threading.Lock()
-    _data_url = (
-        "https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip"
-    )
-
-    def frame_by_iso3(self, iso3):
-        df = self._instance
-        return df[df.ISO_A3 == iso3]
-
-
 def country_geom(iso3: str, crs: MaybeCRS = None) -> Geometry:
     """
     Extract geometry for a country from geopandas sample data.
     """
     # pylint: disable=import-outside-toplevel
+    import geopandas as gpd
+
     from ..converters import from_geopandas
 
-    countries = Countries()
-    (gg,) = from_geopandas(countries.frame_by_iso3(iso3))
+    with catch_warnings():
+        filterwarnings("ignore", category=FutureWarning)
+        df = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
+    (gg,) = from_geopandas(df[df.iso_a3 == iso3])
     crs = norm_crs(crs)
     if crs is not None:
         gg = gg.to_crs(crs)
diff --git a/odc/geo/math.py b/odc/geo/math.py
index a809272..ac70992 100644
--- a/odc/geo/math.py
+++ b/odc/geo/math.py
@@ -26,7 +26,6 @@ import numpy as np
 from affine import Affine
 from numpy.polynomial.polynomial import polygrid2d, polyval2d
 
-from ._interop import is_dask_collection
 from .types import (
     XY,
     AnchorEnum,
@@ -287,7 +286,7 @@ def data_resolution_and_offset(
 
     :returns: ``(resolution, offset)``
     """
-    if is_dask_collection(data):
+    if not isinstance(data, np.ndarray):
         data = data.values
 
     if data.size < 2:


Test patch:
diff --git a/tests/test_converters.py b/tests/test_converters.py
index dc45ba7..33a73b9 100644
--- a/tests/test_converters.py
+++ b/tests/test_converters.py
@@ -7,6 +7,7 @@ import pytest
 
 rasterio = pytest.importorskip("rasterio")
 gpd = pytest.importorskip("geopandas")
+gpd_datasets = pytest.importorskip("geopandas.datasets")
 
 from odc.geo._interop import have
 from odc.geo.converters import extract_gcps, from_geopandas, map_crs, rio_geobox
@@ -15,26 +16,27 @@ from odc.geo.geobox import GeoBox, GeoBoxBase
 
 
 @pytest.fixture
-def countries_geodataframe():
-    from odc.geo.data import Countries
+def ne_lowres_path():
+    with catch_warnings():
+        filterwarnings("ignore")
+        path = gpd_datasets.get_path("naturalearth_lowres")
+    yield path
 
-    yield Countries()._instance
 
-
-def test_from_geopandas(countries_geodataframe):
-    df = countries_geodataframe
+def test_from_geopandas(ne_lowres_path):
+    df = gpd.read_file(ne_lowres_path)
     gg = from_geopandas(df)
     assert isinstance(gg, list)
     assert len(gg) == len(df)
     assert gg[0].crs == "epsg:4326"
 
-    (au,) = from_geopandas(df[df.ISO_A3 == "AUS"].to_crs(epsg=3577))
+    (au,) = from_geopandas(df[df.iso_a3 == "AUS"].to_crs(epsg=3577))
     assert au.crs.epsg == 3577
 
-    (au,) = from_geopandas(df[df.ISO_A3 == "AUS"].to_crs(epsg=3857).geometry)
+    (au,) = from_geopandas(df[df.iso_a3 == "AUS"].to_crs(epsg=3857).geometry)
     assert au.crs.epsg == 3857
 
-    assert from_geopandas(df.CONTINENT) == []
+    assert from_geopandas(df.continent) == []
 
 
 def test_have():
diff --git a/tests/test_gcp.py b/tests/test_gcp.py
index 1f9a2cf..3d263cd 100644
--- a/tests/test_gcp.py
+++ b/tests/test_gcp.py
@@ -125,8 +125,8 @@ def test_gcp_geobox_xr(au_gcp_geobox: GCPGeoBox):
     assert _gbox.crs == gbox.crs
     assert (_gbox.extent ^ gbox.extent).is_empty
 
-    # corrupt gcps
-    yy.spatial_ref.attrs["gcps"] = "not even geojson"
+    # corrupt some gcps
+    yy.spatial_ref.attrs["gcps"]["features"][0].pop("properties")
     # should not throw, just return None
     assert yy.odc.uncached.geobox is None