Odc geo follow throughs

datacube/api/core.py

@@ -2,6 +2,7 @@
 #
 # Copyright (c) 2015-2021 ODC Contributors
 # SPDX-License-Identifier: Apache-2.0
+import warnings
 import uuid
 import collections.abc
 from itertools import groupby
@@ -15,7 +16,7 @@
 from datacube.config import LocalConfig
 from datacube.storage import reproject_and_fuse, BandInfo
 from datacube.utils import ignore_exceptions_if
-from odc.geo import CRS, XY, Resolution
+from odc.geo import CRS, yx_, res_, resyx_
 from odc.geo.xr import xr_coords
 from datacube.utils.dates import normalise_dt
 from odc.geo.geom import intersects, box, bbox_union
@@ -280,7 +281,7 @@ def load(self, product=None, measurements=None, output_crs=None, resolution=None
                         y=(-35.15, -35.2),
                         time=('1990', '1991'),
                         output_crs='EPSG:3577`,
-                        resolution=(-30, 30),
+                        resolution=30,
                         resampling='cubic'
                 )
 
@@ -307,9 +308,9 @@ def load(self, product=None, measurements=None, output_crs=None, resolution=None
             This differs from the ``crs`` parameter desribed above, which is used to define the CRS
             of the coordinates in the query itself.
 
-        :param (float,float) resolution:
-            A tuple of the spatial resolution of the returned data. Units are in the coordinate
-            space of ``output_crs``.
+        :param float|(float,float) resolution:
+            A the spatial resolution of the returned data, provided as a tuple if not using square pixels
+            with inverted Y axis. Units are in the coordinate space of ``output_crs``.
 
             This includes the direction (as indicated by a positive or negative number).
             For most CRSs, the first number will be negative, e.g. ``(-30, 30)``.
@@ -430,6 +431,12 @@ def filter_jan(dataset): return dataset.time.begin.month == 1
             if extra_dims.has_empty_dim():
                 return xarray.Dataset()
 
+        # if resolution is square, raise warning to use a single number instead of a tuple
+        if type(resolution) is tuple and resolution[0] == -resolution[1]:
+            warnings.warn("When using square pixels with an inverted Y axis, "
+                          "provide resolution as a single number instead of a tuple.")
+            resolution = resolution[1]
+
         geobox = output_geobox(like=like, output_crs=output_crs, resolution=resolution, align=align,
                                grid_spec=datacube_product.grid_spec,
                                load_hints=datacube_product.load_hints(),
@@ -908,11 +915,10 @@ def output_geobox(like=None, output_crs=None, resolution=None, align=None,
 
             geopolygon = get_bounds(datasets, crs)
 
-    if resolution is not None and type(resolution) is not Resolution:
-        resolution = Resolution(*resolution)
+    resolution = resyx_(*resolution) if type(resolution) is tuple else res_(resolution)
 
     if align is not None:
-        align = XY(*align)
+        align = yx_(align)
 
     return GeoBox.from_geopolygon(geopolygon, resolution, crs, align)
 

datacube/api/grid_workflow.py

@@ -172,7 +172,7 @@ def cell_observations(self, cell_index=None, geopolygon=None, tile_buffer=None,
         :param odc.geo.geom.Geometry geopolygon:
             Only return observations with data inside polygon.
         :param (float,float) tile_buffer:
-            buffer tiles by (y, x) in CRS units
+            buffer tiles by (x, y) in CRS units
         :param (int,int) cell_index:
             The cell index. E.g. (14, -40)
         :param indexers:
@@ -188,6 +188,8 @@ def cell_observations(self, cell_index=None, geopolygon=None, tile_buffer=None,
         # pylint: disable=too-many-locals
         # TODO: split this method into 3: cell/polygon/unconstrained querying
 
+        if tile_buffer is not None:
+            _LOG.warning("tile_buffer is now expected in (x, y) order")
         if tile_buffer is not None and geopolygon is not None:
             raise ValueError('Cannot process tile_buffering and geopolygon together.')
         cells = {}

datacube/drivers/netcdf/writer.py

@@ -13,11 +13,11 @@
 import numpy
 
 from datacube.utils.masking import describe_flags_def
-from datacube.utils import data_resolution_and_offset
 from netCDF4 import Dataset
 
 from odc.geo import CRS
 from odc.geo.geom import box
+from odc.geo.math import data_resolution_and_offset
 
 from datacube import __version__
 

datacube/model/__init__.py

@@ -34,9 +34,10 @@
     "ExtraDimensions", "IngestorConfig"
 ]
 
-from odc.geo import CRS, BoundingBox, Geometry, wh_
+from odc.geo import CRS, BoundingBox, Geometry, wh_, xy_
 from odc.geo.geobox import GeoBox
 from odc.geo.geom import intersects, polygon
+import odc.geo.gridspec as gridspec
 from deprecat import deprecat
 
 _LOG = logging.getLogger(__name__)
@@ -524,7 +525,7 @@ def extra_dimensions(self) -> "ExtraDimensions":
         return ExtraDimensions(self._extra_dimensions)
 
     @cached_property
-    def grid_spec(self) -> Optional['GridSpec']:
+    def grid_spec(self) -> Optional['gridspec.GridSpec']:
         """
         Grid specification for this product
         """
@@ -549,7 +550,13 @@ def extract_point(name):
         if not complete:
             return None
 
-        return GridSpec(crs=crs, **gs_params)
+        # tile_size param has been renamed to tile_shape in odc-geo version of GridSpec
+        gs_params['tile_shape'] = gs_params.pop('tile_size')
+
+        # convert origin to XY
+        gs_params['origin'] = xy_(gs_params['origin'])
+
+        return gridspec.GridSpec(crs=crs, **gs_params)
 
     @staticmethod
     def validate_extra_dims(definition: Mapping[str, Any]):
@@ -712,7 +719,7 @@ def to_dict(self) -> Mapping[str, Any]:
             row.update({
                 'crs': str(self.grid_spec.crs),
                 'spatial_dimensions': self.grid_spec.dimensions,
-                'tile_size': self.grid_spec.tile_size,
+                'tile_size': self.grid_spec.tile_shape,
                 'resolution': self.grid_spec.resolution,
             })
         return row
@@ -750,6 +757,7 @@ class IngestorConfig:
     pass
 
 
+@deprecat(reason='GridSpec is now defined in odc-geo.', version='1.9.0')
 class GridSpec:
     """
     Definition for a regular spatial grid

datacube/storage/_read.py

@@ -7,7 +7,7 @@
 import numpy as np
 from typing import Optional, Tuple
 
-from ..utils.math import is_almost_int, valid_mask
+from ..utils.math import valid_mask
 
 from odc.geo import wh_
 from odc.geo.roi import (
@@ -26,6 +26,7 @@
     Nodata,
 )
 from odc.geo.overlap import compute_reproject_roi, is_affine_st
+from odc.geo.math import is_almost_int
 
 
 def rdr_geobox(rdr) -> GeoBox:

datacube/testutils/geom.py

@@ -10,7 +10,7 @@
 from odc.geo import CRS
 from odc.geo.geobox import GeoBox
 from odc.geo.math import apply_affine
-from datacube.model import GridSpec
+from odc.geo.gridspec import GridSpec
 
 # pylint: disable=invalid-name
 
@@ -19,8 +19,8 @@
 epsg3857 = CRS('EPSG:3857')
 
 AlbersGS = GridSpec(crs=epsg3577,
-                    tile_size=(100000.0, 100000.0),
-                    resolution=(-25, 25),
+                    tile_shape=(100000.0, 100000.0),
+                    resolution=25,
                     origin=(0.0, 0.0))
 
 SAMPLE_WKT_WITHOUT_AUTHORITY = '''PROJCS["unnamed",

datacube/utils/__init__.py

@@ -31,7 +31,6 @@
     unsqueeze_data_array,
     spatial_dims,
     iter_slices,
-    data_resolution_and_offset,
 )
 from ._misc import (
     DatacubeException,
@@ -66,7 +65,6 @@
     "unsqueeze_dataset",
     "spatial_dims",
     "iter_slices",
-    "data_resolution_and_offset",
     "DatacubeException",
     "schema_validated",
     "write_user_secret_file",

datacube/utils/geometry/gbox.py

@@ -12,7 +12,7 @@
 
 from . import Geometry, GeoBox, BoundingBox
 from .tools import align_up
-from datacube.utils.math import clamp
+from odc.geo.math import clamp
 
 # pylint: disable=invalid-name
 MaybeInt = Optional[int]

datacube/utils/geometry/tools.py

@@ -446,7 +446,7 @@ def box_overlap(src_shape, dst_shape, ST, tol):  # noqa: N803
                       direction is: Xsrc = ST*Xdst
     :param       tol: Sub-pixel translation tolerance that's scaled by resolution.
     """
-    from datacube.utils.math import maybe_int, snap_scale
+    from odc.geo.math import maybe_int, snap_scale
 
     (sx, _, tx,
      _, sy, ty,

datacube/utils/math.py

@@ -3,11 +3,13 @@
 # Copyright (c) 2015-2020 ODC Contributors
 # SPDX-License-Identifier: Apache-2.0
 from typing import Tuple, Union, Optional, Any, cast
-from math import ceil, fmod
+from math import ceil
 
 import numpy
 import xarray as xr
-from affine import Affine
+import odc.geo.math as geomath
+
+from deprecat import deprecat
 
 
 def unsqueeze_data_array(da: xr.DataArray,
@@ -70,70 +72,29 @@ def spatial_dims(xx: Union[xr.DataArray, xr.Dataset],
     return None
 
 
+@deprecat(reason='This method has been moved to odc-geo.', version='1.9.0')
 def maybe_zero(x: float, tol: float) -> float:
-    """ Turn almost zeros to actual zeros
-    """
-    if abs(x) < tol:
-        return 0
-    return x
+    return geomath.maybe_zero(x, tol)
 
 
+@deprecat(reason='This method has been moved to odc-geo.', version='1.9.0')
 def maybe_int(x: float, tol: float) -> Union[int, float]:
-    """ Turn almost ints to actual ints, pass through other values unmodified
-    """
-    def split(x):
-        x_part = fmod(x, 1.0)
-        x_whole = x - x_part
-        if x_part > 0.5:
-            x_part -= 1
-            x_whole += 1
-        elif x_part < -0.5:
-            x_part += 1
-            x_whole -= 1
-        return (x_whole, x_part)
-
-    x_whole, x_part = split(x)
-
-    if abs(x_part) < tol:  # almost int
-        return int(x_whole)
-    else:
-        return x
+    return geomath.maybe_int(x, tol)
 
 
+@deprecat(reason='This method has been moved to odc-geo.', version='1.9.0')
 def snap_scale(s, tol=1e-6):
-    """ Snap scale to the nearest integer and simple fractions in the form 1/<int>
-    """
-    if abs(s) >= 1 - tol:
-        return maybe_int(s, tol)
-    else:
-        # Check of s is 0
-        if abs(s) < tol:
-            return s
-
-        # Check for simple fractions
-        s_inv = 1 / s
-        s_inv_snapped = maybe_int(s_inv, tol)
-        if s_inv_snapped is s_inv:
-            return s
-        return 1 / s_inv_snapped
+    return geomath.snap_scale(s, tol)
 
 
+@deprecat(reason='This method has been moved to odc-geo.', version='1.9.0')
 def clamp(x, lo, up):
-    """
-    clamp x to be lo <= x <= up
-    """
-    assert lo <= up
-    return lo if x < lo else up if x > up else x
+    return geomath.clamp(x, lo, up)
 
 
+@deprecat(reason='This method has been moved to odc-geo.', version='1.9.0')
 def is_almost_int(x: float, tol: float):
-    """
-    Check if number is close enough to an integer
-    """
-    x = abs(fmod(x, 1))
-    if x > 0.5:
-        x = 1 - x
-    return x < tol
+    return geomath.is_almost_int(x, tol)
 
 
 def dtype_is_float(dtype) -> bool:
@@ -199,65 +160,14 @@ def num2numpy(x, dtype, ignore_range=None):
     return None
 
 
+@deprecat(reason='This method has been moved to odc-geo.', version='1.9.0')
 def data_resolution_and_offset(data, fallback_resolution=None):
-    """ Compute resolution and offset from x/y axis data.
-
-        Only uses first two coordinate values, assumes that data is regularly
-        sampled.
-
-        Returns
-        =======
-        (resolution: float, offset: float)
-    """
-    if data.size < 2:
-        if data.size < 1:
-            raise ValueError("Can't calculate resolution for empty data")
-        if fallback_resolution is None:
-            raise ValueError("Can't calculate resolution with data size < 2")
-        res = fallback_resolution
-    else:
-        res = (data[data.size - 1] - data[0]) / (data.size - 1.0)
-        res = res.item()
-
-    off = data[0] - 0.5 * res
-    return res, off.item()
+    return geomath.data_resolution_and_offset(data, fallback_resolution)
 
 
+@deprecat(reason='This method has been moved to odc-geo.', version='1.9.0')
 def affine_from_axis(xx, yy, fallback_resolution=None):
-    """ Compute Affine transform from pixel to real space given X,Y coordinates.
-
-        :param xx: X axis coordinates
-        :param yy: Y axis coordinates
-        :param fallback_resolution: None|float|(resx:float, resy:float) resolution to
-                                    assume for single element axis.
-
-        (0, 0) in pixel space is defined as top left corner of the top left pixel
-            \
-            `` 0   1
-             +---+---+
-           0 |   |   |
-             +---+---+
-           1 |   |   |
-             +---+---+
-
-        Only uses first two coordinate values, assumes that data is regularly
-        sampled.
-
-        raises ValueError when any axis is empty
-        raises ValueError when any axis has single value and fallback resolution was not supplied.
-    """
-    if fallback_resolution is not None:
-        if isinstance(fallback_resolution, (float, int)):
-            frx, fry = fallback_resolution, fallback_resolution
-        else:
-            frx, fry = fallback_resolution
-    else:
-        frx, fry = None, None
-
-    xres, xoff = data_resolution_and_offset(xx, frx)
-    yres, yoff = data_resolution_and_offset(yy, fry)
-
-    return Affine.translation(xoff, yoff) * Affine.scale(xres, yres)
+    return geomath.affine_from_axis(xx, yy, fallback_resolution)
 
 
 def iter_slices(shape, chunk_size):

datacube/utils/xarray_geoextensions.py

@@ -16,7 +16,7 @@
 import warnings
 import xarray
 from datacube.utils import spatial_dims
-from datacube.utils.math import affine_from_axis
+from odc.geo.math import affine_from_axis
 from odc.geo._xr_interop import _xarray_geobox as _xr_geobox