Repository: opendatacube/odc-geo
Pull request #171 (2024-07-11 07:50:32)

Update so tests pass with latest versions of upstream libraries.
1. Directly download naturalearth countries dataset.  (See #165)

     Noting that the country dataset from naturalearth used previously is NOT packaged with geodatasets - so we may as well just hardcode the url and download it directly.

2. Refactor use of `np.find_common_type` (which is not available in Numpy 2.x) to `np.result_type` (which is).

3. Check if the `gcps` metadata field on the xarray is of string type, and if so run `json.loads` over it.  (Required for `rioxarray>=0.15.6` as it now stores the gcps definition as a string)



Patch:
diff --git a/odc/geo/_blocks.py b/odc/geo/_blocks.py
index cc78b0b..84b9580 100644
--- a/odc/geo/_blocks.py
+++ b/odc/geo/_blocks.py
@@ -13,9 +13,15 @@ from .roi import VariableSizedTiles, roi_intersect3, roi_normalise, roi_shape
 from .types import Chunks2d
 
 
-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)
+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)
 
 
 class BlockAssembler:
@@ -33,7 +39,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
@@ -126,7 +132,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/_xr_interop.py b/odc/geo/_xr_interop.py
index 6404d7b..53a62f9 100644
--- a/odc/geo/_xr_interop.py
+++ b/odc/geo/_xr_interop.py
@@ -8,6 +8,7 @@ Add ``.odc.`` extension to :py:class:`xarray.Dataset` and :class:`xarray.DataArr
 from __future__ import annotations
 
 import functools
+import json
 import math
 import warnings
 from dataclasses import dataclass
@@ -209,7 +210,8 @@ def _mk_crs_coord(
     crs_wkt = cf.get("crs_wkt", None) or crs.wkt
 
     if gcps is not None:
-        cf["gcps"] = _gcps_to_json(gcps)
+        # Store as string
+        cf["gcps"] = json.dumps(_gcps_to_json(gcps))
 
     if transform is not None:
         cf["GeoTransform"] = _render_geo_transform(transform, precision=24)
@@ -523,13 +525,15 @@ 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):
+    except (IndexError, KeyError, ValueError, json.JSONDecodeError):
         return None
 
 
diff --git a/odc/geo/data/__init__.py b/odc/geo/data/__init__.py
index 2d91b7d..c9433b4 100644
--- a/odc/geo/data/__init__.py
+++ b/odc/geo/data/__init__.py
@@ -1,5 +1,7 @@
 import json
 import lzma
+import threading
+
 from functools import lru_cache
 from pathlib import Path
 from typing import Any, Dict, Optional, Tuple
@@ -46,19 +48,54 @@ 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
 
-    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])
+    countries = Countries()
+    (gg,) = from_geopandas(countries.frame_by_iso3(iso3))
     crs = norm_crs(crs)
     if crs is not None:
         gg = gg.to_crs(crs)


Test patch:
diff --git a/tests/test_converters.py b/tests/test_converters.py
index 33a73b9..dc45ba7 100644
--- a/tests/test_converters.py
+++ b/tests/test_converters.py
@@ -7,7 +7,6 @@ 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
@@ -16,27 +15,26 @@ from odc.geo.geobox import GeoBox, GeoBoxBase
 
 
 @pytest.fixture
-def ne_lowres_path():
-    with catch_warnings():
-        filterwarnings("ignore")
-        path = gpd_datasets.get_path("naturalearth_lowres")
-    yield path
+def countries_geodataframe():
+    from odc.geo.data import Countries
 
+    yield Countries()._instance
 
-def test_from_geopandas(ne_lowres_path):
-    df = gpd.read_file(ne_lowres_path)
+
+def test_from_geopandas(countries_geodataframe):
+    df = countries_geodataframe
     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 3d263cd..1f9a2cf 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 some gcps
-    yy.spatial_ref.attrs["gcps"]["features"][0].pop("properties")
+    # corrupt gcps
+    yy.spatial_ref.attrs["gcps"] = "not even geojson"
     # should not throw, just return None
     assert yy.odc.uncached.geobox is None