Skip to content

Commit

Permalink
Merge pull request #171 from opendatacube/switch_to_geodatasets
Browse files Browse the repository at this point in the history
Update so tests pass with latest versions of upstream libraries.
  • Loading branch information
SpacemanPaul authored Jul 16, 2024
2 parents d020f21 + 65a16cb commit c1585d3
Show file tree
Hide file tree
Showing 5 changed files with 71 additions and 26 deletions.
16 changes: 11 additions & 5 deletions odc/geo/_blocks.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,15 @@
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:
Expand All @@ -33,7 +39,7 @@ def __init__(
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
Expand Down Expand Up @@ -126,7 +132,7 @@ def extract(
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)

Expand Down
8 changes: 6 additions & 2 deletions odc/geo/_xr_interop.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
from __future__ import annotations

import functools
import json
import math
import warnings
from dataclasses import dataclass
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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


Expand Down
49 changes: 43 additions & 6 deletions odc/geo/data/__init__.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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)
Expand Down
20 changes: 9 additions & 11 deletions tests/test_converters.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@

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
Expand All @@ -16,27 +15,26 @@


@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():
Expand Down
4 changes: 2 additions & 2 deletions tests/test_gcp.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down

0 comments on commit c1585d3

Please sign in to comment.