Repository: opendatacube/odc-geo
Pull request #2 (2022-01-12 05:28:02)

Refactor: split overly large `_base.py` file
What was a single `_base.py` is now

- `_crs.py` (pyproj CRS wrapper)
- `_geom.py` ("Shapely Geometry with CRS" class and helpers)
- `_geobox.py` (`GeoBox` class and related utilities)

This had to be done carefully to preserve git history across the split, hence two commits (step 1|2)

Checking history is preserved on all 3 new files:

- https://github.com/opendatacube/odc-geo/blame/refactor/odc/geo/_geobox.py
- https://github.com/opendatacube/odc-geo/blame/refactor/odc/geo/_geom.py
- https://github.com/opendatacube/odc-geo/blame/refactor/odc/geo/_crs.py

Note that there is a weak intentional cyclic dependency between `_crs.py` and
`_geom.py`. Mostly `_geom` needs CRS, that's a fundamental dependency. But CRS
also uses Geometry classes for some convenience methods, namely to report "valid
area" for a given CRS object.



codecov[bot]: # [Codecov](https://codecov.io/gh/opendatacube/odc-geo/pull/2?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) Report
> Merging [#2](https://codecov.io/gh/opendatacube/odc-geo/pull/2?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) (4b8cb7a) into [develop](https://codecov.io/gh/opendatacube/odc-geo/commit/c51b8e9dd935af13d453e3640739dc529858e9bd?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) (c51b8e9) will **decrease** coverage by `0.05%`.
> The diff coverage is `91.24%`.

[![Impacted file tree graph](https://codecov.io/gh/opendatacube/odc-geo/pull/2/graphs/tree.svg?width=650&height=150&src=pr&token=Qc7isqXNGF&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube)](https://codecov.io/gh/opendatacube/odc-geo/pull/2?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube)

```diff
@@             Coverage Diff             @@
##           develop       #2      +/-   ##
===========================================
- Coverage    84.60%   84.55%   -0.06%     
===========================================
  Files           10       12       +2     
  Lines         1507     1528      +21     
===========================================
+ Hits          1275     1292      +17     
- Misses         232      236       +4     
```


| [Impacted Files](https://codecov.io/gh/opendatacube/odc-geo/pull/2?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) | Coverage Δ | |
|---|---|---|
| [odc/geo/xarray\_geoextensions.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby94YXJyYXlfZ2VvZXh0ZW5zaW9ucy5weQ==) | `0.00% <0.00%> (ø)` | |
| [odc/geo/\_geobox.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9fZ2VvYm94LnB5) | `91.93% <91.93%> (ø)` | |
| [odc/geo/\_crs.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9fY3JzLnB5) | `97.46% <97.46%> (ø)` | |
| [odc/geo/\_\_init\_\_.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |
| [odc/geo/\_geom.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9fZ2VvbS5weQ==) | `98.30% <100.00%> (ø)` | |
| [odc/geo/\_gridspec.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9fZ3JpZHNwZWMucHk=) | `95.89% <100.00%> (+0.11%)` | :arrow_up: |
| [odc/geo/\_warp.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9fd2FycC5weQ==) | `29.78% <100.00%> (ø)` | |
| [odc/geo/gbox.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9nYm94LnB5) | `99.09% <100.00%> (+<0.01%)` | :arrow_up: |
| [odc/geo/math.py](https://codecov.io/gh/opendatacube/odc-geo/pull/2/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-b2RjL2dlby9tYXRoLnB5) | `37.89% <100.00%> (+2.04%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/opendatacube/odc-geo/pull/2?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/opendatacube/odc-geo/pull/2?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube). Last update [c51b8e9...4b8cb7a](https://codecov.io/gh/opendatacube/odc-geo/pull/2?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).

Kirill888: Right now almost everything is exposed from one flat namespace, hence the `_`. But we should think about having better hierarchy than just flat namespace.

Patch:
diff --git a/odc/geo/__init__.py b/odc/geo/__init__.py
index e75f40a..cd8a9f0 100644
--- a/odc/geo/__init__.py
+++ b/odc/geo/__init__.py
@@ -4,30 +4,33 @@
 # SPDX-License-Identifier: Apache-2.0
 """ Geometric shapes and operations on them
 """
+# import order is important here
+#  _crs <-- _geom <-- _geobox <- other
+# isort: skip_file
 
 from ._version import __version__
 
-from ._base import (
+from ._crs import (
     CRS,
-    BoundingBox,
-    Coordinate,
-    CoordList,
     CRSError,
     CRSMismatchError,
-    GeoBox,
-    Geometry,
     MaybeCRS,
     SomeCRS,
-    assign_crs,
+    crs_units_per_degree,
+    norm_crs,
+    norm_crs_or_error,
+)
+
+from ._geom import (
+    BoundingBox,
+    CoordList,
+    Geometry,
     bbox_intersection,
     bbox_union,
     box,
     chop_along_antimeridian,
     clip_lon180,
     common_crs,
-    crs_units_per_degree,
-    geobox_intersection_conservative,
-    geobox_union_conservative,
     intersects,
     line,
     lonlat_bounds,
@@ -37,15 +40,25 @@ from ._base import (
     multipoint,
     multipolygon,
     point,
-    polygon,
     polygon_from_transform,
+    polygon,
     projected_lon,
-    scaled_down_geobox,
     sides,
     unary_intersection,
     unary_union,
 )
-from ._warp import rio_reproject, warp_affine
+
+from ._geobox import (
+    Coordinate,
+    GeoBox,
+    assign_crs,
+    geobox_intersection_conservative,
+    geobox_union_conservative,
+    scaled_down_geobox,
+)
+
+from ._gridspec import GridSpec
+
 from .tools import (
     affine_from_pts,
     apply_affine,
@@ -70,70 +83,72 @@ from .tools import (
     split_translation,
     w_,
 )
+from ._warp import rio_reproject, warp_affine
 
-from ._gridspec import GridSpec
 
 __all__ = [
-    "Coordinate",
+    "__version__",
+    "affine_from_pts",
+    "apply_affine",
+    "assign_crs",
+    "bbox_intersection",
+    "bbox_union",
     "BoundingBox",
+    "box",
+    "chop_along_antimeridian",
+    "clip_lon180",
+    "common_crs",
+    "compute_axis_overlap",
+    "compute_reproject_roi",
+    "Coordinate",
+    "CoordList",
+    "crs_units_per_degree",
+    "CRS",
     "CRSError",
     "CRSMismatchError",
-    "CRS",
-    "MaybeCRS",
-    "SomeCRS",
-    "CoordList",
-    "Geometry",
+    "decompose_rws",
+    "geobox_intersection_conservative",
+    "geobox_union_conservative",
     "GeoBox",
+    "Geometry",
+    "get_scale_at_point",
     "GridSpec",
-    "assign_crs",
-    "common_crs",
-    "bbox_union",
-    "bbox_intersection",
-    "crs_units_per_degree",
-    "geobox_union_conservative",
-    "geobox_intersection_conservative",
     "intersects",
-    "point",
-    "multipoint",
+    "is_affine_st",
     "line",
+    "lonlat_bounds",
+    "MaybeCRS",
+    "mid_longitude",
+    "multigeom",
     "multiline",
-    "polygon",
+    "multipoint",
     "multipolygon",
-    "multigeom",
-    "box",
-    "sides",
+    "native_pix_transform",
+    "norm_crs",
+    "norm_crs_or_error",
+    "point",
     "polygon_from_transform",
-    "unary_union",
-    "unary_intersection",
-    "lonlat_bounds",
+    "polygon",
     "projected_lon",
-    "clip_lon180",
-    "chop_along_antimeridian",
-    "mid_longitude",
-    "is_affine_st",
-    "apply_affine",
-    "compute_axis_overlap",
+    "rio_reproject",
     "roi_boundary",
+    "roi_center",
+    "roi_from_points",
+    "roi_intersect",
     "roi_is_empty",
     "roi_is_full",
-    "roi_intersect",
-    "roi_shape",
     "roi_normalise",
-    "roi_from_points",
-    "roi_center",
     "roi_pad",
+    "roi_shape",
     "scaled_down_geobox",
-    "scaled_down_shape",
     "scaled_down_roi",
+    "scaled_down_shape",
     "scaled_up_roi",
-    "decompose_rws",
-    "affine_from_pts",
-    "get_scale_at_point",
-    "native_pix_transform",
-    "compute_reproject_roi",
+    "sides",
+    "SomeCRS",
     "split_translation",
-    "warp_affine",
-    "rio_reproject",
+    "unary_intersection",
+    "unary_union",
     "w_",
-    "__version__",
+    "warp_affine",
 ]
diff --git a/odc/geo/_crs.py b/odc/geo/_crs.py
new file mode 100644
index 0000000..536b14d
--- /dev/null
+++ b/odc/geo/_crs.py
@@ -0,0 +1,303 @@
+# This file is part of the Open Data Cube, see https://opendatacube.org for more information
+#
+# Copyright (c) 2015-2020 ODC Contributors
+# SPDX-License-Identifier: Apache-2.0
+import math
+import warnings
+from distutils.version import LooseVersion
+from typing import TYPE_CHECKING, Any, Callable, Dict, Optional, Tuple, Union
+
+import cachetools
+import numpy
+import rasterio
+from pyproj.crs import CRS as _CRS
+from pyproj.enums import WktVersion
+from pyproj.exceptions import CRSError
+from pyproj.transformer import Transformer
+
+SomeCRS = Union[str, "CRS", _CRS, Dict[str, Any]]
+MaybeCRS = Optional[SomeCRS]
+
+
+@cachetools.cached({})
+def _make_crs(crs_str: str) -> Tuple[_CRS, Optional[int]]:
+    crs = _CRS.from_user_input(crs_str)
+    return (crs, crs.to_epsg())
+
+
+def _make_crs_transform_key(from_crs, to_crs, always_xy):
+    return (id(from_crs), id(to_crs), always_xy)
+
+
+@cachetools.cached({}, key=_make_crs_transform_key)
+def _make_crs_transform(from_crs, to_crs, always_xy):
+    return Transformer.from_crs(from_crs, to_crs, always_xy=always_xy).transform
+
+
+def _guess_crs_str(crs_spec: Any) -> Optional[str]:
+    """
+    Returns a string representation of the crs spec.
+    Returns `None` if it does not understand the spec.
+    """
+    if isinstance(crs_spec, str):
+        return crs_spec
+    if isinstance(crs_spec, dict):
+        crs_spec = _CRS.from_dict(crs_spec)
+
+    if hasattr(crs_spec, "to_epsg"):
+        epsg = crs_spec.to_epsg()
+        if epsg is not None:
+            return f"EPSG:{crs_spec.to_epsg()}"
+    if hasattr(crs_spec, "to_wkt"):
+        return crs_spec.to_wkt()
+    return None
+
+
+class CRS:
+    """
+    Wrapper around `pyproj.CRS` for backwards compatibility.
+    """
+
+    DEFAULT_WKT_VERSION = (
+        WktVersion.WKT1_GDAL
+        if LooseVersion(rasterio.__gdal_version__) < LooseVersion("3.0.0")
+        else WktVersion.WKT2_2019
+    )
+
+    __slots__ = ("_crs", "_epsg", "_str")
+
+    def __init__(self, crs_str: Any):
+        """
+        :param crs_str: string representation of a CRS, often an EPSG code like 'EPSG:4326'
+        :raises: `pyproj.exceptions.CRSError`
+        """
+        crs_str = _guess_crs_str(crs_str)
+        if crs_str is None:
+            raise CRSError(
+                "Expect string or any object with `.to_epsg()` or `.to_wkt()` method"
+            )
+
+        _crs, _epsg = _make_crs(crs_str)
+
+        self._crs = _crs
+        self._epsg = _epsg
+        self._str = crs_str
+
+    def __getstate__(self):
+        return {"crs_str": self._str}
+
+    def __setstate__(self, state):
+        self.__init__(state["crs_str"])
+
+    def to_wkt(self, pretty: bool = False, version: Optional[WktVersion] = None) -> str:
+        """
+        WKT representation of the CRS
+        """
+        if version is None:
+            version = self.DEFAULT_WKT_VERSION
+
+        return self._crs.to_wkt(pretty=pretty, version=version)
+
+    @property
+    def wkt(self) -> str:
+        return self.to_wkt(version=WktVersion.WKT1_GDAL)
+
+    def to_epsg(self) -> Optional[int]:
+        """
+        EPSG Code of the CRS or None
+        """
+        return self._epsg
+
+    @property
+    def epsg(self) -> Optional[int]:
+        return self._epsg
+
+    @property
+    def semi_major_axis(self):
+        return self._crs.ellipsoid.semi_major_metre
+
+    @property
+    def semi_minor_axis(self):
+        return self._crs.ellipsoid.semi_minor_metre
+
+    @property
+    def inverse_flattening(self):
+        return self._crs.ellipsoid.inverse_flattening
+
+    @property
+    def geographic(self) -> bool:
+        return self._crs.is_geographic
+
+    @property
+    def projected(self) -> bool:
+        return self._crs.is_projected
+
+    @property
+    def dimensions(self) -> Tuple[str, str]:
+        """
+        List of dimension names of the CRS.
+        The ordering of the names is intended to reflect the `numpy` array axis order of the loaded raster.
+        """
+        if self.geographic:
+            return "latitude", "longitude"
+
+        if self.projected:
+            return "y", "x"
+
+        raise ValueError("Neither projected nor geographic")  # pragma: no cover
+
+    @property
+    def units(self) -> Tuple[str, str]:
+        """
+        List of dimension units of the CRS.
+        The ordering of the units is intended to reflect the `numpy` array axis order of the loaded raster.
+        """
+        if self.geographic:
+            return "degrees_north", "degrees_east"
+
+        if self.projected:
+            x, y = self._crs.axis_info
+            return x.unit_name, y.unit_name
+
+        raise ValueError("Neither projected nor geographic")  # pragma: no cover
+
+    def __str__(self) -> str:
+        return self._str
+
+    def __hash__(self) -> int:
+        return hash(self.to_wkt())
+
+    def __repr__(self) -> str:
+        return f"CRS('{self._str}')"
+
+    def __eq__(self, other) -> bool:
+        if not isinstance(other, CRS):
+            try:
+                other = CRS(other)
+            except Exception:  # pylint: disable=broad-except
+                return False
+
+        if self._crs is other._crs:
+            return True
+
+        if self.epsg is not None and other.epsg is not None:
+            return self.epsg == other.epsg
+
+        return self._crs == other._crs
+
+    def __ne__(self, other) -> bool:
+        return not self == other
+
+    @property
+    def proj(self) -> _CRS:
+        """Access proj.CRS object that this wraps"""
+        return self._crs
+
+    @property
+    def valid_region(self) -> Optional["_geom.Geometry"]:
+        """Return valid region of this CRS.
+
+        Bounding box in Lon/Lat as a 4 point Polygon in EPSG:4326.
+        None if not defined
+        """
+
+        from . import _geom  # pylint: disable=import-outside-toplevel
+
+        aou = self._crs.area_of_use
+        if aou is None:
+            return None
+        return _geom.box(aou.west, aou.south, aou.east, aou.north, "EPSG:4326")
+
+    @property
+    def crs_str(self) -> str:
+        """DEPRECATED"""
+        warnings.warn(
+            "Please use `str(crs)` instead of `crs.crs_str`",
+            category=DeprecationWarning,
+        )
+        return self._str
+
+    def transformer_to_crs(
+        self, other: "CRS", always_xy=True
+    ) -> Callable[[Any, Any], Tuple[Any, Any]]:
+        """
+        Returns a function that maps x, y -> x', y' where x, y are coordinates in
+        this stored either as scalars or ndarray objects and x', y' are the same
+        points in the `other` CRS.
+        """
+
+        # pylint: disable=protected-access
+        transform = _make_crs_transform(self._crs, other._crs, always_xy=always_xy)
+
+        def result(x, y):
+            rx, ry = transform(x, y)  # pylint: disable=unpacking-non-sequence
+
+            if not isinstance(rx, numpy.ndarray) or not isinstance(ry, numpy.ndarray):
+                return (rx, ry)
+
+            missing = numpy.isnan(rx) | numpy.isnan(ry)
+            rx[missing] = numpy.nan
+            ry[missing] = numpy.nan
+            return (rx, ry)
+
+        return result
+
+
+class CRSMismatchError(ValueError):
+    """
+    Raised when geometry operation is attempted on geometries in different
+    coordinate references.
+    """
+
+
+def norm_crs(crs: MaybeCRS) -> Optional[CRS]:
+    """Normalise CRS representation."""
+    if isinstance(crs, CRS):
+        return crs
+    if crs is None:
+        return None
+    return CRS(crs)
+
+
+def norm_crs_or_error(crs: MaybeCRS) -> CRS:
+    """Normalise CRS representation, raise error if input is ``None``."""
+    if isinstance(crs, CRS):
+        return crs
+    if crs is None:
+        raise ValueError("Expect valid CRS")
+    return CRS(crs)
+
+
+def crs_units_per_degree(
+    crs: SomeCRS,
+    lon: Union[float, Tuple[float, float]],
+    lat: float = 0,
+    step: float = 0.1,
+) -> float:
+    """Compute number of CRS units per degree for a projected CRS at a given location
+    in lon/lat.
+
+    Location can be supplied as a tuple or as two arguments.
+
+    Returns
+    -------
+    A floating number S such that `S*degrees -> meters`
+    """
+
+    from . import _geom  # pylint: disable=import-outside-toplevel
+
+    if isinstance(lon, tuple):
+        lon, lat = lon
+
+    lon2 = lon + step
+    if lon2 > 180:
+        lon2 = lon - step
+
+    ll = _geom.line([(lon, lat), (lon2, lat)], "EPSG:4326")
+    xy = ll.to_crs(crs, resolution=math.inf)
+
+    return xy.length / step
+
+
+if TYPE_CHECKING:
+    from . import _geom
diff --git a/odc/geo/_geobox.py b/odc/geo/_geobox.py
new file mode 100644
index 0000000..d917dcd
--- /dev/null
+++ b/odc/geo/_geobox.py
@@ -0,0 +1,445 @@
+# This file is part of the Open Data Cube, see https://opendatacube.org for more information
+#
+# Copyright (c) 2015-2020 ODC Contributors
+# SPDX-License-Identifier: Apache-2.0
+import math
+from collections import OrderedDict, namedtuple
+from typing import Dict, Hashable, List, Optional, Tuple, Union
+
+import numpy
+import xarray as xr
+from affine import Affine
+
+from ._crs import CRS, MaybeCRS, norm_crs_or_error
+from ._geom import (
+    BoundingBox,
+    Geometry,
+    bbox_intersection,
+    bbox_union,
+    polygon_from_transform,
+)
+from .math import is_almost_int
+from .tools import is_affine_st, roi_normalise, roi_shape
+
+Coordinate = namedtuple("Coordinate", ("values", "units", "resolution"))
+
+
+def _align_pix(left: float, right: float, res: float, off: float) -> Tuple[float, int]:
+    if res < 0:
+        res = -res
+        val = math.ceil((right - off) / res) * res + off
+        width = max(1, int(math.ceil((val - left - 0.1 * res) / res)))
+    else:
+        val = math.floor((left - off) / res) * res + off
+        width = max(1, int(math.ceil((right - val - 0.1 * res) / res)))
+    return val, width
+
+
+class GeoBox:
+    """
+    Defines the location and resolution of a rectangular grid of data,
+    including it's :py:class:`CRS`.
+
+    :param crs: Coordinate Reference System
+    :param affine: Affine transformation defining the location of the geobox
+    """
+
+    def __init__(self, width: int, height: int, affine: Affine, crs: MaybeCRS):
+        assert is_affine_st(
+            affine
+        ), "Only axis-aligned geoboxes are currently supported"
+        self.width = width
+        self.height = height
+        self.affine = affine
+        self.extent = polygon_from_transform(width, height, affine, crs=crs)
+
+    @classmethod
+    def from_geopolygon(
+        cls,
+        geopolygon: Geometry,
+        resolution: Tuple[float, float],
+        crs: MaybeCRS = None,
+        align: Optional[Tuple[float, float]] = None,
+    ) -> "GeoBox":
+        """
+        :param resolution: (y_resolution, x_resolution)
+        :param crs: CRS to use, if different from the geopolygon
+        :param align: Align geobox such that point 'align' lies on the pixel boundary.
+        """
+        align = align or (0.0, 0.0)
+        assert (
+            0.0 <= align[1] <= abs(resolution[1])
+        ), "X align must be in [0, abs(x_resolution)] range"
+        assert (
+            0.0 <= align[0] <= abs(resolution[0])
+        ), "Y align must be in [0, abs(y_resolution)] range"
+
+        if crs is None:
+            crs = geopolygon.crs
+        else:
+            geopolygon = geopolygon.to_crs(crs)
+
+        bounding_box = geopolygon.boundingbox
+        offx, width = _align_pix(
+            bounding_box.left, bounding_box.right, resolution[1], align[1]
+        )
+        offy, height = _align_pix(
+            bounding_box.bottom, bounding_box.top, resolution[0], align[0]
+        )
+        affine = Affine.translation(offx, offy) * Affine.scale(
+            resolution[1], resolution[0]
+        )
+        return GeoBox(crs=crs, affine=affine, width=width, height=height)
+
+    def buffered(self, ybuff, xbuff) -> "GeoBox":
+        """
+        Produce a tile buffered by ybuff, xbuff (in CRS units)
+        """
+        by, bx = (
+            _round_to_res(buf, res) for buf, res in zip((ybuff, xbuff), self.resolution)
+        )
+        affine = self.affine * Affine.translation(-bx, -by)
+
+        return GeoBox(
+            width=self.width + 2 * bx,
+            height=self.height + 2 * by,
+            affine=affine,
+            crs=self.crs,
+        )
+
+    def __getitem__(self, roi) -> "GeoBox":
+        if isinstance(roi, int):
+            roi = (slice(roi, roi + 1), slice(None, None))
+
+        if isinstance(roi, slice):
+            roi = (roi, slice(None, None))
+
+        if len(roi) > 2:
+            raise ValueError("Expect 2d slice")
+
+        if not all(s.step is None or s.step == 1 for s in roi):
+            raise NotImplementedError("scaling not implemented, yet")
+
+        roi = roi_normalise(roi, self.shape)
+        ty, tx = (s.start for s in roi)
+        h, w = roi_shape(roi)
+
+        affine = self.affine * Affine.translation(tx, ty)
+
+        return GeoBox(width=w, height=h, affine=affine, crs=self.crs)
+
+    def __or__(self, other) -> "GeoBox":
+        """A geobox that encompasses both self and other."""
+        return geobox_union_conservative([self, other])
+
+    def __and__(self, other) -> "GeoBox":
+        """A geobox that is contained in both self and other."""
+        return geobox_intersection_conservative([self, other])
+
+    def is_empty(self) -> bool:
+        return self.width == 0 or self.height == 0
+
+    def __bool__(self) -> bool:
+        return not self.is_empty()
+
+    def __hash__(self):
+        return hash((*self.shape, self.crs, self.affine))
+
+    @property
+    def transform(self) -> Affine:
+        return self.affine
+
+    @property
+    def shape(self) -> Tuple[int, int]:
+        return self.height, self.width
+
+    @property
+    def crs(self) -> Optional[CRS]:
+        return self.extent.crs
+
+    @property
+    def dimensions(self) -> Tuple[str, str]:
+        """
+        List of dimension names of the GeoBox
+        """
+        crs = self.crs
+        if crs is None:
+            return ("y", "x")
+        return crs.dimensions
+
+    @property
+    def resolution(self) -> Tuple[float, float]:
+        """
+        Resolution in Y,X dimensions
+        """
+        return self.affine.e, self.affine.a
+
+    @property
+    def alignment(self) -> Tuple[float, float]:
+        """
+        Alignment of pixel boundaries in Y,X dimensions
+        """
+        return self.affine.yoff % abs(self.affine.e), self.affine.xoff % abs(
+            self.affine.a
+        )
+
+    @property
+    def coordinates(self) -> Dict[str, Coordinate]:
+        """
+        dict of coordinate labels
+        """
+        yres, xres = self.resolution
+        yoff, xoff = self.affine.yoff, self.affine.xoff
+
+        xs = numpy.arange(self.width) * xres + (xoff + xres / 2)
+        ys = numpy.arange(self.height) * yres + (yoff + yres / 2)
+
+        units = self.crs.units if self.crs is not None else ("1", "1")
+
+        return OrderedDict(
+            (dim, Coordinate(labels, units, res))
+            for dim, labels, units, res in zip(
+                self.dimensions, (ys, xs), units, (yres, xres)
+            )
+        )
+
+    def xr_coords(
+        self, with_crs: Union[bool, str] = False
+    ) -> Dict[Hashable, xr.DataArray]:
+        """Dictionary of Coordinates in xarray format
+
+        :param with_crs: If True include netcdf/cf style CRS Coordinate
+        with default name 'spatial_ref', if with_crs is a string then treat
+        the string as a name of the coordinate.
+
+        Returns
+        =======
+
+        OrderedDict name:str -> xr.DataArray
+
+        where names are either `y,x` for projected or `latitude, longitude` for geographic.
+
+        """
+        spatial_ref = "spatial_ref"
+        if isinstance(with_crs, str):
+            spatial_ref = with_crs
+            with_crs = True
+
+        attrs = {}
+        coords = self.coordinates
+        crs = self.crs
+        if crs is not None:
+            attrs["crs"] = str(crs)
+
+        coords: Dict[Hashable, xr.DataArray] = {
+            n: _coord_to_xr(n, c, **attrs) for n, c in coords.items()
+        }
+
+        if with_crs and crs is not None:
+            coords[spatial_ref] = _mk_crs_coord(crs, spatial_ref)
+
+        return coords
+
+    @property
+    def geographic_extent(self) -> Geometry:
+        """GeoBox extent in EPSG:4326"""
+        if self.crs is None or self.crs.geographic:
+            return self.extent
+        return self.extent.to_crs(CRS("EPSG:4326"))
+
+    coords = coordinates
+    dims = dimensions
+
+    def __str__(self):
+        return f"GeoBox({self.geographic_extent})"
+
+    def __repr__(self):
+        return f"GeoBox({self.width}, {self.height}, {self.affine!r}, {self.crs})"
+
+    def __eq__(self, other):
+        if not isinstance(other, GeoBox):
+            return False
+
+        return (
+            self.shape == other.shape
+            and self.transform == other.transform
+            and self.crs == other.crs
+        )
+
+
+def bounding_box_in_pixel_domain(geobox: GeoBox, reference: GeoBox) -> BoundingBox:
+    """
+    Returns the bounding box of `geobox` with respect to the pixel grid
+    defined by `reference` when their coordinate grids are compatible,
+    that is, have the same CRS, same pixel size and orientation, and
+    are related by whole pixel translation,
+    otherwise raises `ValueError`.
+    """
+    tol = 1.0e-8
+
+    if reference.crs != geobox.crs:
+        raise ValueError("Cannot combine geoboxes in different CRSs")
+
+    a, b, c, d, e, f, *_ = ~reference.affine * geobox.affine
+
+    if not (
+        numpy.isclose(a, 1)
+        and numpy.isclose(b, 0)
+        and is_almost_int(c, tol)
+        and numpy.isclose(d, 0)
+        and numpy.isclose(e, 1)
+        and is_almost_int(f, tol)
+    ):
+        raise ValueError("Incompatible grids")
+
+    tx, ty = round(c), round(f)
+    return BoundingBox(tx, ty, tx + geobox.width, ty + geobox.height)
+
+
+def geobox_union_conservative(geoboxes: List[GeoBox]) -> GeoBox:
+    """Union of geoboxes. Fails whenever incompatible grids are encountered."""
+    if len(geoboxes) == 0:
+        raise ValueError("No geoboxes supplied")
+
+    reference, *_ = geoboxes
+
+    bbox = bbox_union(
+        bounding_box_in_pixel_domain(geobox, reference=reference) for geobox in geoboxes
+    )
+
+    affine = reference.affine * Affine.translation(*bbox[:2])
+
+    return GeoBox(
+        width=bbox.width, height=bbox.height, affine=affine, crs=reference.crs
+    )
+
+
+def geobox_intersection_conservative(geoboxes: List[GeoBox]) -> GeoBox:
+    """
+    Intersection of geoboxes. Fails whenever incompatible grids are encountered.
+    """
+    if len(geoboxes) == 0:
+        raise ValueError("No geoboxes supplied")
+
+    reference, *_ = geoboxes
+
+    bbox = bbox_intersection(
+        bounding_box_in_pixel_domain(geobox, reference=reference) for geobox in geoboxes
+    )
+
+    # standardise empty geobox representation
+    if bbox.left > bbox.right:
+        bbox = BoundingBox(
+            left=bbox.left, bottom=bbox.bottom, right=bbox.left, top=bbox.top
+        )
+    if bbox.bottom > bbox.top:
+        bbox = BoundingBox(
+            left=bbox.left, bottom=bbox.bottom, right=bbox.right, top=bbox.bottom
+        )
+
+    affine = reference.affine * Affine.translation(*bbox[:2])
+
+    return GeoBox(
+        width=bbox.width, height=bbox.height, affine=affine, crs=reference.crs
+    )
+
+
+def scaled_down_geobox(src_geobox: GeoBox, scaler: int) -> GeoBox:
+    """Given a source geobox and integer scaler compute geobox of a scaled down image.
+
+    Output geobox will be padded when shape is not a multiple of scaler.
+    Example: 5x4, scaler=2 -> 3x2
+
+    NOTE: here we assume that pixel coordinates are 0,0 at the top-left
+          corner of a top-left pixel.
+
+    """
+    assert scaler > 1
+
+    H, W = (X // scaler + (1 if X % scaler else 0) for X in src_geobox.shape)
+
+    # Since 0,0 is at the corner of a pixel, not center, there is no
+    # translation between pixel plane coords due to scaling
+    A = src_geobox.transform * Affine.scale(scaler, scaler)
+
+    return GeoBox(W, H, A, src_geobox.crs)
+
+
+def _round_to_res(value: float, res: float) -> int:
+    res = abs(res)
+    return int(math.ceil((value - 0.1 * res) / res))
+
+
+def _mk_crs_coord(crs: CRS, name: str = "spatial_ref") -> xr.DataArray:
+    # pylint: disable=protected-access
+
+    if crs.projected:
+        grid_mapping_name = crs._crs.to_cf().get("grid_mapping_name")
+        if grid_mapping_name is None:
+            grid_mapping_name = "??"
+        grid_mapping_name = grid_mapping_name.lower()
+    else:
+        grid_mapping_name = "latitude_longitude"
+
+    epsg = 0 if crs.epsg is None else crs.epsg
+
+    return xr.DataArray(
+        numpy.asarray(epsg, "int32"),
+        name=name,
+        dims=(),
+        attrs={"spatial_ref": crs.wkt, "grid_mapping_name": grid_mapping_name},
+    )
+
+
+def _coord_to_xr(name: str, c: Coordinate, **attrs) -> xr.DataArray:
+    """Construct xr.DataArray from named Coordinate object, this can then be used
+    to define coordinates for xr.Dataset|xr.DataArray
+    """
+    attrs = dict(units=c.units, resolution=c.resolution, **attrs)
+    return xr.DataArray(c.values, coords={name: c.values}, dims=(name,), attrs=attrs)
+
+
+def assign_crs(
+    xx: Union[xr.DataArray, xr.Dataset],
+    crs: MaybeCRS = None,
+    crs_coord_name: str = "spatial_ref",
+) -> Union[xr.Dataset, xr.DataArray]:
+    """
+    Assign CRS for a non-georegistered array or dataset.
+
+    Returns a new object with CRS information populated.
+
+    Can also be called without ``crs`` argument on data that already has CRS
+    information but not in the format used by datacube, in this case CRS
+    metadata will be restructured into a shape used by datacube. This format
+    allows for better propagation of CRS information through various
+    computations.
+
+    .. code-block:: python
+
+        xx = odc.geo.assign_crs(xr.open_rasterio("some-file.tif"))
+        print(xx.geobox)
+        print(xx.astype("float32").geobox)
+
+
+    :param xx: :py:class:`xarray.Dataset` or :py:class:`~xarray.DataArray`
+    :param crs: CRS to assign, if omitted try to guess from attributes
+    :param crs_coord_name: how to name crs corodinate (defaults to ``spatial_ref``)
+    """
+    if crs is None:
+        geobox = getattr(xx, "geobox", None)
+        if geobox is None:
+            raise ValueError("Failed to guess CRS for this object")
+        crs = geobox.crs
+
+    crs = norm_crs_or_error(crs)
+    crs_coord = _mk_crs_coord(crs, name=crs_coord_name)
+    xx = xx.assign_coords({crs_coord.name: crs_coord})
+
+    xx.attrs.update(grid_mapping=crs_coord_name)
+
+    if isinstance(xx, xr.Dataset):
+        for band in xx.data_vars.values():
+            band.attrs.update(grid_mapping=crs_coord_name)
+
+    return xx
diff --git a/odc/geo/_base.py b/odc/geo/_geom.py
similarity index 53%
rename from odc/geo/_base.py
rename to odc/geo/_geom.py
index 2f22a80..4bb7429 100644
--- a/odc/geo/_base.py
+++ b/odc/geo/_geom.py
@@ -6,42 +6,25 @@ import array
 import functools
 import itertools
 import math
-import warnings
-from collections import OrderedDict, namedtuple
+from collections import namedtuple
 from collections.abc import Sequence
-from distutils.version import LooseVersion
-from typing import (
-    Any,
-    Callable,
-    Dict,
-    Hashable,
-    Iterable,
-    Iterator,
-    List,
-    Optional,
-    Tuple,
-    Union,
-)
+from typing import Any, Dict, Iterable, Iterator, List, Optional, Tuple, Union
 
-import cachetools
 import numpy
-import rasterio
-import xarray as xr
 from affine import Affine
-from pyproj.crs import CRS as _CRS
-from pyproj.enums import WktVersion
-from pyproj.exceptions import CRSError
-from pyproj.transformer import Transformer
 from shapely import geometry, ops
 from shapely.geometry import base
 
-from .math import is_almost_int
-from .tools import is_affine_st, roi_normalise, roi_shape
+from ._crs import (
+    CRS,
+    CRSMismatchError,
+    MaybeCRS,
+    SomeCRS,
+    norm_crs,
+    norm_crs_or_error,
+)
 
-Coordinate = namedtuple("Coordinate", ("values", "units", "resolution"))
 _BoundingBox = namedtuple("_BoundingBox", ("left", "bottom", "right", "top"))
-SomeCRS = Union[str, "CRS", _CRS, Dict[str, Any]]
-MaybeCRS = Optional[SomeCRS]
 CoordList = List[Tuple[float, float]]
 
 # pylint: disable=too-many-lines
@@ -126,250 +109,6 @@ class BoundingBox(_BoundingBox):
         return BoundingBox.from_xy((p1[0], p2[0]), (p1[1], p2[1]))
 
 
-@cachetools.cached({})
-def _make_crs(crs_str: str) -> Tuple[_CRS, Optional[int]]:
-    crs = _CRS.from_user_input(crs_str)
-    return (crs, crs.to_epsg())
-
-
-def _make_crs_transform_key(from_crs, to_crs, always_xy):
-    return (id(from_crs), id(to_crs), always_xy)
-
-
-@cachetools.cached({}, key=_make_crs_transform_key)
-def _make_crs_transform(from_crs, to_crs, always_xy):
-    return Transformer.from_crs(from_crs, to_crs, always_xy=always_xy).transform
-
-
-def _guess_crs_str(crs_spec: Any) -> Optional[str]:
-    """
-    Returns a string representation of the crs spec.
-    Returns `None` if it does not understand the spec.
-    """
-    if isinstance(crs_spec, str):
-        return crs_spec
-    if isinstance(crs_spec, dict):
-        crs_spec = _CRS.from_dict(crs_spec)
-
-    if hasattr(crs_spec, "to_epsg"):
-        epsg = crs_spec.to_epsg()
-        if epsg is not None:
-            return f"EPSG:{crs_spec.to_epsg()}"
-    if hasattr(crs_spec, "to_wkt"):
-        return crs_spec.to_wkt()
-    return None
-
-
-class CRS:
-    """
-    Wrapper around `pyproj.CRS` for backwards compatibility.
-    """
-
-    DEFAULT_WKT_VERSION = (
-        WktVersion.WKT1_GDAL
-        if LooseVersion(rasterio.__gdal_version__) < LooseVersion("3.0.0")
-        else WktVersion.WKT2_2019
-    )
-
-    __slots__ = ("_crs", "_epsg", "_str")
-
-    def __init__(self, crs_str: Any):
-        """
-        :param crs_str: string representation of a CRS, often an EPSG code like 'EPSG:4326'
-        :raises: `pyproj.exceptions.CRSError`
-        """
-        crs_str = _guess_crs_str(crs_str)
-        if crs_str is None:
-            raise CRSError(
-                "Expect string or any object with `.to_epsg()` or `.to_wkt()` method"
-            )
-
-        _crs, _epsg = _make_crs(crs_str)
-
-        self._crs = _crs
-        self._epsg = _epsg
-        self._str = crs_str
-
-    def __getstate__(self):
-        return {"crs_str": self._str}
-
-    def __setstate__(self, state):
-        self.__init__(state["crs_str"])
-
-    def to_wkt(self, pretty: bool = False, version: Optional[WktVersion] = None) -> str:
-        """
-        WKT representation of the CRS
-        """
-        if version is None:
-            version = self.DEFAULT_WKT_VERSION
-
-        return self._crs.to_wkt(pretty=pretty, version=version)
-
-    @property
-    def wkt(self) -> str:
-        return self.to_wkt(version=WktVersion.WKT1_GDAL)
-
-    def to_epsg(self) -> Optional[int]:
-        """
-        EPSG Code of the CRS or None
-        """
-        return self._epsg
-
-    @property
-    def epsg(self) -> Optional[int]:
-        return self._epsg
-
-    @property
-    def semi_major_axis(self):
-        return self._crs.ellipsoid.semi_major_metre
-
-    @property
-    def semi_minor_axis(self):
-        return self._crs.ellipsoid.semi_minor_metre
-
-    @property
-    def inverse_flattening(self):
-        return self._crs.ellipsoid.inverse_flattening
-
-    @property
-    def geographic(self) -> bool:
-        return self._crs.is_geographic
-
-    @property
-    def projected(self) -> bool:
-        return self._crs.is_projected
-
-    @property
-    def dimensions(self) -> Tuple[str, str]:
-        """
-        List of dimension names of the CRS.
-        The ordering of the names is intended to reflect the `numpy` array axis order of the loaded raster.
-        """
-        if self.geographic:
-            return "latitude", "longitude"
-
-        if self.projected:
-            return "y", "x"
-
-        raise ValueError("Neither projected nor geographic")  # pragma: no cover
-
-    @property
-    def units(self) -> Tuple[str, str]:
-        """
-        List of dimension units of the CRS.
-        The ordering of the units is intended to reflect the `numpy` array axis order of the loaded raster.
-        """
-        if self.geographic:
-            return "degrees_north", "degrees_east"
-
-        if self.projected:
-            x, y = self._crs.axis_info
-            return x.unit_name, y.unit_name
-
-        raise ValueError("Neither projected nor geographic")  # pragma: no cover
-
-    def __str__(self) -> str:
-        return self._str
-
-    def __hash__(self) -> int:
-        return hash(self.to_wkt())
-
-    def __repr__(self) -> str:
-        return f"CRS('{self._str}')"
-
-    def __eq__(self, other) -> bool:
-        if not isinstance(other, CRS):
-            try:
-                other = CRS(other)
-            except Exception:  # pylint: disable=broad-except
-                return False
-
-        if self._crs is other._crs:
-            return True
-
-        if self.epsg is not None and other.epsg is not None:
-            return self.epsg == other.epsg
-
-        return self._crs == other._crs
-
-    def __ne__(self, other) -> bool:
-        return not self == other
-
-    @property
-    def proj(self) -> _CRS:
-        """Access proj.CRS object that this wraps"""
-        return self._crs
-
-    @property
-    def valid_region(self) -> Optional["Geometry"]:
-        """Return valid region of this CRS.
-
-        Bounding box in Lon/Lat as a 4 point Polygon in EPSG:4326.
-        None if not defined
-        """
-        aou = self._crs.area_of_use
-        if aou is None:
-            return None
-        return box(aou.west, aou.south, aou.east, aou.north, "EPSG:4326")
-
-    @property
-    def crs_str(self) -> str:
-        """DEPRECATED"""
-        warnings.warn(
-            "Please use `str(crs)` instead of `crs.crs_str`",
-            category=DeprecationWarning,
-        )
-        return self._str
-
-    def transformer_to_crs(
-        self, other: "CRS", always_xy=True
-    ) -> Callable[[Any, Any], Tuple[Any, Any]]:
-        """
-        Returns a function that maps x, y -> x', y' where x, y are coordinates in
-        this stored either as scalars or ndarray objects and x', y' are the same
-        points in the `other` CRS.
-        """
-
-        # pylint: disable=protected-access
-        transform = _make_crs_transform(self._crs, other._crs, always_xy=always_xy)
-
-        def result(x, y):
-            rx, ry = transform(x, y)  # pylint: disable=unpacking-non-sequence
-
-            if not isinstance(rx, numpy.ndarray) or not isinstance(ry, numpy.ndarray):
-                return (rx, ry)
-
-            missing = numpy.isnan(rx) | numpy.isnan(ry)
-            rx[missing] = numpy.nan
-            ry[missing] = numpy.nan
-            return (rx, ry)
-
-        return result
-
-
-class CRSMismatchError(ValueError):
-    """
-    Raised when geometry operation is attempted on geometries in different
-    coordinate references.
-    """
-
-
-def _norm_crs(crs: MaybeCRS) -> Optional[CRS]:
-    if isinstance(crs, CRS):
-        return crs
-    if crs is None:
-        return None
-    return CRS(crs)
-
-
-def _norm_crs_or_error(crs: MaybeCRS) -> CRS:
-    if isinstance(crs, CRS):
-        return crs
-    if crs is None:
-        raise ValueError("Expect valid CRS")
-    return CRS(crs)
-
-
 def wrap_shapely(method):
     """
     Takes a method that expects shapely geometry arguments
@@ -464,7 +203,7 @@ class Geometry:
             self.geom: base.BaseGeometry = _clone_shapely_geom(geom.geom)
             return
 
-        crs = _norm_crs(crs)
+        crs = norm_crs(crs)
         self.crs = crs
         if isinstance(geom, base.BaseGeometry):
             self.geom = geom
@@ -693,7 +432,7 @@ class Geometry:
                                   when converting to geographic projections.
                                   Currently only works in few specific cases (source CRS is smooth over the dateline).
         """
-        crs = _norm_crs_or_error(crs)
+        crs = norm_crs_or_error(crs)
         if self.crs == crs:
             return self
 
@@ -776,7 +515,7 @@ def projected_lon(
     step: float = 1.0,
 ) -> Geometry:
     """Project vertical line along some longitude into given CRS."""
-    crs = _norm_crs_or_error(crs)
+    crs = norm_crs_or_error(crs)
     yy = numpy.arange(lat[0], lat[1], step, dtype="float32")
     xx = numpy.full_like(yy, lon)
     tr = CRS("EPSG:4326").transformer_to_crs(crs)
@@ -1008,352 +747,6 @@ def unary_intersection(geoms: Iterable[Geometry]) -> Geometry:
     return functools.reduce(Geometry.intersection, geoms)
 
 
-def _align_pix(left: float, right: float, res: float, off: float) -> Tuple[float, int]:
-    if res < 0:
-        res = -res
-        val = math.ceil((right - off) / res) * res + off
-        width = max(1, int(math.ceil((val - left - 0.1 * res) / res)))
-    else:
-        val = math.floor((left - off) / res) * res + off
-        width = max(1, int(math.ceil((right - val - 0.1 * res) / res)))
-    return val, width
-
-
-class GeoBox:
-    """
-    Defines the location and resolution of a rectangular grid of data,
-    including it's :py:class:`CRS`.
-
-    :param crs: Coordinate Reference System
-    :param affine: Affine transformation defining the location of the geobox
-    """
-
-    def __init__(self, width: int, height: int, affine: Affine, crs: MaybeCRS):
-        assert is_affine_st(
-            affine
-        ), "Only axis-aligned geoboxes are currently supported"
-        self.width = width
-        self.height = height
-        self.affine = affine
-        self.extent = polygon_from_transform(width, height, affine, crs=crs)
-
-    @classmethod
-    def from_geopolygon(
-        cls,
-        geopolygon: Geometry,
-        resolution: Tuple[float, float],
-        crs: MaybeCRS = None,
-        align: Optional[Tuple[float, float]] = None,
-    ) -> "GeoBox":
-        """
-        :param resolution: (y_resolution, x_resolution)
-        :param crs: CRS to use, if different from the geopolygon
-        :param align: Align geobox such that point 'align' lies on the pixel boundary.
-        """
-        align = align or (0.0, 0.0)
-        assert (
-            0.0 <= align[1] <= abs(resolution[1])
-        ), "X align must be in [0, abs(x_resolution)] range"
-        assert (
-            0.0 <= align[0] <= abs(resolution[0])
-        ), "Y align must be in [0, abs(y_resolution)] range"
-
-        if crs is None:
-            crs = geopolygon.crs
-        else:
-            geopolygon = geopolygon.to_crs(crs)
-
-        bounding_box = geopolygon.boundingbox
-        offx, width = _align_pix(
-            bounding_box.left, bounding_box.right, resolution[1], align[1]
-        )
-        offy, height = _align_pix(
-            bounding_box.bottom, bounding_box.top, resolution[0], align[0]
-        )
-        affine = Affine.translation(offx, offy) * Affine.scale(
-            resolution[1], resolution[0]
-        )
-        return GeoBox(crs=crs, affine=affine, width=width, height=height)
-
-    def buffered(self, ybuff, xbuff) -> "GeoBox":
-        """
-        Produce a tile buffered by ybuff, xbuff (in CRS units)
-        """
-        by, bx = (
-            _round_to_res(buf, res) for buf, res in zip((ybuff, xbuff), self.resolution)
-        )
-        affine = self.affine * Affine.translation(-bx, -by)
-
-        return GeoBox(
-            width=self.width + 2 * bx,
-            height=self.height + 2 * by,
-            affine=affine,
-            crs=self.crs,
-        )
-
-    def __getitem__(self, roi) -> "GeoBox":
-        if isinstance(roi, int):
-            roi = (slice(roi, roi + 1), slice(None, None))
-
-        if isinstance(roi, slice):
-            roi = (roi, slice(None, None))
-
-        if len(roi) > 2:
-            raise ValueError("Expect 2d slice")
-
-        if not all(s.step is None or s.step == 1 for s in roi):
-            raise NotImplementedError("scaling not implemented, yet")
-
-        roi = roi_normalise(roi, self.shape)
-        ty, tx = (s.start for s in roi)
-        h, w = roi_shape(roi)
-
-        affine = self.affine * Affine.translation(tx, ty)
-
-        return GeoBox(width=w, height=h, affine=affine, crs=self.crs)
-
-    def __or__(self, other) -> "GeoBox":
-        """A geobox that encompasses both self and other."""
-        return geobox_union_conservative([self, other])
-
-    def __and__(self, other) -> "GeoBox":
-        """A geobox that is contained in both self and other."""
-        return geobox_intersection_conservative([self, other])
-
-    def is_empty(self) -> bool:
-        return self.width == 0 or self.height == 0
-
-    def __bool__(self) -> bool:
-        return not self.is_empty()
-
-    def __hash__(self):
-        return hash((*self.shape, self.crs, self.affine))
-
-    @property
-    def transform(self) -> Affine:
-        return self.affine
-
-    @property
-    def shape(self) -> Tuple[int, int]:
-        return self.height, self.width
-
-    @property
-    def crs(self) -> Optional[CRS]:
-        return self.extent.crs
-
-    @property
-    def dimensions(self) -> Tuple[str, str]:
-        """
-        List of dimension names of the GeoBox
-        """
-        crs = self.crs
-        if crs is None:
-            return ("y", "x")
-        return crs.dimensions
-
-    @property
-    def resolution(self) -> Tuple[float, float]:
-        """
-        Resolution in Y,X dimensions
-        """
-        return self.affine.e, self.affine.a
-
-    @property
-    def alignment(self) -> Tuple[float, float]:
-        """
-        Alignment of pixel boundaries in Y,X dimensions
-        """
-        return self.affine.yoff % abs(self.affine.e), self.affine.xoff % abs(
-            self.affine.a
-        )
-
-    @property
-    def coordinates(self) -> Dict[str, Coordinate]:
-        """
-        dict of coordinate labels
-        """
-        yres, xres = self.resolution
-        yoff, xoff = self.affine.yoff, self.affine.xoff
-
-        xs = numpy.arange(self.width) * xres + (xoff + xres / 2)
-        ys = numpy.arange(self.height) * yres + (yoff + yres / 2)
-
-        units = self.crs.units if self.crs is not None else ("1", "1")
-
-        return OrderedDict(
-            (dim, Coordinate(labels, units, res))
-            for dim, labels, units, res in zip(
-                self.dimensions, (ys, xs), units, (yres, xres)
-            )
-        )
-
-    def xr_coords(
-        self, with_crs: Union[bool, str] = False
-    ) -> Dict[Hashable, xr.DataArray]:
-        """Dictionary of Coordinates in xarray format
-
-        :param with_crs: If True include netcdf/cf style CRS Coordinate
-        with default name 'spatial_ref', if with_crs is a string then treat
-        the string as a name of the coordinate.
-
-        Returns
-        =======
-
-        OrderedDict name:str -> xr.DataArray
-
-        where names are either `y,x` for projected or `latitude, longitude` for geographic.
-
-        """
-        spatial_ref = "spatial_ref"
-        if isinstance(with_crs, str):
-            spatial_ref = with_crs
-            with_crs = True
-
-        attrs = {}
-        coords = self.coordinates
-        crs = self.crs
-        if crs is not None:
-            attrs["crs"] = str(crs)
-
-        coords: Dict[Hashable, xr.DataArray] = {
-            n: _coord_to_xr(n, c, **attrs) for n, c in coords.items()
-        }
-
-        if with_crs and crs is not None:
-            coords[spatial_ref] = _mk_crs_coord(crs, spatial_ref)
-
-        return coords
-
-    @property
-    def geographic_extent(self) -> Geometry:
-        """GeoBox extent in EPSG:4326"""
-        if self.crs is None or self.crs.geographic:
-            return self.extent
-        return self.extent.to_crs(CRS("EPSG:4326"))
-
-    coords = coordinates
-    dims = dimensions
-
-    def __str__(self):
-        return f"GeoBox({self.geographic_extent})"
-
-    def __repr__(self):
-        return f"GeoBox({self.width}, {self.height}, {self.affine!r}, {self.crs})"
-
-    def __eq__(self, other):
-        if not isinstance(other, GeoBox):
-            return False
-
-        return (
-            self.shape == other.shape
-            and self.transform == other.transform
-            and self.crs == other.crs
-        )
-
-
-def bounding_box_in_pixel_domain(geobox: GeoBox, reference: GeoBox) -> BoundingBox:
-    """
-    Returns the bounding box of `geobox` with respect to the pixel grid
-    defined by `reference` when their coordinate grids are compatible,
-    that is, have the same CRS, same pixel size and orientation, and
-    are related by whole pixel translation,
-    otherwise raises `ValueError`.
-    """
-    tol = 1.0e-8
-
-    if reference.crs != geobox.crs:
-        raise ValueError("Cannot combine geoboxes in different CRSs")
-
-    a, b, c, d, e, f, *_ = ~reference.affine * geobox.affine
-
-    if not (
-        numpy.isclose(a, 1)
-        and numpy.isclose(b, 0)
-        and is_almost_int(c, tol)
-        and numpy.isclose(d, 0)
-        and numpy.isclose(e, 1)
-        and is_almost_int(f, tol)
-    ):
-        raise ValueError("Incompatible grids")
-
-    tx, ty = round(c), round(f)
-    return BoundingBox(tx, ty, tx + geobox.width, ty + geobox.height)
-
-
-def geobox_union_conservative(geoboxes: List[GeoBox]) -> GeoBox:
-    """Union of geoboxes. Fails whenever incompatible grids are encountered."""
-    if len(geoboxes) == 0:
-        raise ValueError("No geoboxes supplied")
-
-    reference, *_ = geoboxes
-
-    bbox = bbox_union(
-        bounding_box_in_pixel_domain(geobox, reference=reference) for geobox in geoboxes
-    )
-
-    affine = reference.affine * Affine.translation(*bbox[:2])
-
-    return GeoBox(
-        width=bbox.width, height=bbox.height, affine=affine, crs=reference.crs
-    )
-
-
-def geobox_intersection_conservative(geoboxes: List[GeoBox]) -> GeoBox:
-    """
-    Intersection of geoboxes. Fails whenever incompatible grids are encountered.
-    """
-    if len(geoboxes) == 0:
-        raise ValueError("No geoboxes supplied")
-
-    reference, *_ = geoboxes
-
-    bbox = bbox_intersection(
-        bounding_box_in_pixel_domain(geobox, reference=reference) for geobox in geoboxes
-    )
-
-    # standardise empty geobox representation
-    if bbox.left > bbox.right:
-        bbox = BoundingBox(
-            left=bbox.left, bottom=bbox.bottom, right=bbox.left, top=bbox.top
-        )
-    if bbox.bottom > bbox.top:
-        bbox = BoundingBox(
-            left=bbox.left, bottom=bbox.bottom, right=bbox.right, top=bbox.bottom
-        )
-
-    affine = reference.affine * Affine.translation(*bbox[:2])
-
-    return GeoBox(
-        width=bbox.width, height=bbox.height, affine=affine, crs=reference.crs
-    )
-
-
-def scaled_down_geobox(src_geobox: GeoBox, scaler: int) -> GeoBox:
-    """Given a source geobox and integer scaler compute geobox of a scaled down image.
-
-    Output geobox will be padded when shape is not a multiple of scaler.
-    Example: 5x4, scaler=2 -> 3x2
-
-    NOTE: here we assume that pixel coordinates are 0,0 at the top-left
-          corner of a top-left pixel.
-
-    """
-    assert scaler > 1
-
-    H, W = (X // scaler + (1 if X % scaler else 0) for X in src_geobox.shape)
-
-    # Since 0,0 is at the corner of a pixel, not center, there is no
-    # translation between pixel plane coords due to scaling
-    A = src_geobox.transform * Affine.scale(scaler, scaler)
-
-    return GeoBox(W, H, A, src_geobox.crs)
-
-
-def _round_to_res(value: float, res: float) -> int:
-    res = abs(res)
-    return int(math.ceil((value - 0.1 * res) / res))
-
-
 def intersects(a: Geometry, b: Geometry) -> bool:
     """Returns True if geometries intersect, else False"""
     return a.intersects(b) and not a.touches(b)
@@ -1393,63 +786,6 @@ def bbox_intersection(bbs: Iterable[BoundingBox]) -> BoundingBox:
     return BoundingBox(L, B, R, T)
 
 
-def _mk_crs_coord(crs: CRS, name: str = "spatial_ref") -> xr.DataArray:
-    # pylint: disable=protected-access
-
-    if crs.projected:
-        grid_mapping_name = crs._crs.to_cf().get("grid_mapping_name")
-        if grid_mapping_name is None:
-            grid_mapping_name = "??"
-        grid_mapping_name = grid_mapping_name.lower()
-    else:
-        grid_mapping_name = "latitude_longitude"
-
-    epsg = 0 if crs.epsg is None else crs.epsg
-
-    return xr.DataArray(
-        numpy.asarray(epsg, "int32"),
-        name=name,
-        dims=(),
-        attrs={"spatial_ref": crs.wkt, "grid_mapping_name": grid_mapping_name},
-    )
-
-
-def _coord_to_xr(name: str, c: Coordinate, **attrs) -> xr.DataArray:
-    """Construct xr.DataArray from named Coordinate object, this can then be used
-    to define coordinates for xr.Dataset|xr.DataArray
-    """
-    attrs = dict(units=c.units, resolution=c.resolution, **attrs)
-    return xr.DataArray(c.values, coords={name: c.values}, dims=(name,), attrs=attrs)
-
-
-def crs_units_per_degree(
-    crs: SomeCRS,
-    lon: Union[float, Tuple[float, float]],
-    lat: float = 0,
-    step: float = 0.1,
-) -> float:
-    """Compute number of CRS units per degree for a projected CRS at a given location
-    in lon/lat.
-
-    Location can be supplied as a tuple or as two arguments.
-
-    Returns
-    -------
-    A floating number S such that `S*degrees -> meters`
-    """
-    if isinstance(lon, tuple):
-        lon, lat = lon
-
-    lon2 = lon + step
-    if lon2 > 180:
-        lon2 = lon - step
-
-    ll = line([(lon, lat), (lon2, lat)], "EPSG:4326")
-    xy = ll.to_crs(crs, resolution=math.inf)
-
-    return xy.length / step
-
-
 def lonlat_bounds(
     geom: Geometry, mode: str = "safe", resolution: Optional[float] = None
 ) -> BoundingBox:
@@ -1494,52 +830,6 @@ def lonlat_bounds(
     return BoundingBox.from_xy(xx_range, bbox.range_y)
 
 
-def assign_crs(
-    xx: Union[xr.DataArray, xr.Dataset],
-    crs: MaybeCRS = None,
-    crs_coord_name: str = "spatial_ref",
-) -> Union[xr.Dataset, xr.DataArray]:
-    """
-    Assign CRS for a non-georegistered array or dataset.
-
-    Returns a new object with CRS information populated.
-
-    Can also be called without ``crs`` argument on data that already has CRS
-    information but not in the format used by datacube, in this case CRS
-    metadata will be restructured into a shape used by datacube. This format
-    allows for better propagation of CRS information through various
-    computations.
-
-    .. code-block:: python
-
-        xx = datacube.utils.geometry.assign_crs(xr.open_rasterio("some-file.tif"))
-        print(xx.geobox)
-        print(xx.astype("float32").geobox)
-
-
-    :param xx:  Dataset or DataArray
-    :param crs: CRS to assign, if omitted try to guess from attributes
-    :param crs_coord_name: how to name crs corodinate (defaults to ``spatial_ref``)
-    """
-    if crs is None:
-        geobox = getattr(xx, "geobox", None)
-        if geobox is None:
-            raise ValueError("Failed to guess CRS for this object")
-        crs = geobox.crs
-
-    crs = _norm_crs_or_error(crs)
-    crs_coord = _mk_crs_coord(crs, name=crs_coord_name)
-    xx = xx.assign_coords({crs_coord.name: crs_coord})
-
-    xx.attrs.update(grid_mapping=crs_coord_name)
-
-    if isinstance(xx, xr.Dataset):
-        for band in xx.data_vars.values():
-            band.attrs.update(grid_mapping=crs_coord_name)
-
-    return xx
-
-
 def mid_longitude(geom: Geometry) -> float:
     """
     Compute longitude of the center point of a geometry
diff --git a/odc/geo/_gridspec.py b/odc/geo/_gridspec.py
index 33ed2e0..94942ec 100644
--- a/odc/geo/_gridspec.py
+++ b/odc/geo/_gridspec.py
@@ -2,22 +2,22 @@
 #
 # Copyright (c) 2015-2020 ODC Contributors
 # SPDX-License-Identifier: Apache-2.0
-"""
-Core classes used across modules.
-"""
+"""GridSpec class."""
 import math
 from typing import Iterator, Optional, Tuple
 
 from affine import Affine
 
-from ._base import CRS, BoundingBox, GeoBox, Geometry, intersects
+from ._crs import CRS
+from ._geobox import GeoBox
+from ._geom import BoundingBox, Geometry, intersects
 
 DEFAULT_SPATIAL_DIMS = ("y", "x")  # Used when product lacks grid_spec
 
 
 class GridSpec:
     """
-    Definition for a regular spatial grid
+    Definition for a regular spatial grid.
 
     >>> gs = GridSpec(crs=CRS('EPSG:4326'), tile_size=(1, 1), resolution=(-0.1, 0.1), origin=(-50.05, 139.95))
     >>> gs.tile_resolution
diff --git a/odc/geo/_warp.py b/odc/geo/_warp.py
index 90893f1..d7c3623 100644
--- a/odc/geo/_warp.py
+++ b/odc/geo/_warp.py
@@ -9,7 +9,7 @@ import rasterio.crs
 import rasterio.warp
 from affine import Affine
 
-from ._base import GeoBox
+from ._geobox import GeoBox
 
 Resampling = Union[str, int, rasterio.warp.Resampling]  # pylint: disable=invalid-name
 Nodata = Optional[Union[int, float]]  # pylint: disable=invalid-name
diff --git a/odc/geo/gbox.py b/odc/geo/gbox.py
index 012a7ac..60b500d 100644
--- a/odc/geo/gbox.py
+++ b/odc/geo/gbox.py
@@ -11,7 +11,8 @@ from typing import Dict, Iterable, Optional, Tuple
 
 from affine import Affine
 
-from ._base import BoundingBox, GeoBox, Geometry
+from ._geobox import GeoBox
+from ._geom import BoundingBox, Geometry
 from .math import clamp
 from .tools import align_up
 
@@ -140,9 +141,9 @@ class GeoboxTiles:
     """Partition GeoBox into sub geoboxes"""
 
     def __init__(self, box: GeoBox, tile_shape: Tuple[int, int]):
-        """Construct from a ``GeoBox``
+        """Construct from a ``GeoBox``.
 
-        :param box: source :class:`datacube.utils.geometry.GeoBox`
+        :param box: source :class:`~odc.geo.GeoBox`
         :param tile_shape: Shape of sub-tiles in pixels (rows, cols)
         """
         self._gbox = box
diff --git a/odc/geo/math.py b/odc/geo/math.py
index becd354..283472b 100644
--- a/odc/geo/math.py
+++ b/odc/geo/math.py
@@ -3,43 +3,12 @@
 # Copyright (c) 2015-2020 ODC Contributors
 # SPDX-License-Identifier: Apache-2.0
 from math import ceil, fmod
-from typing import Optional, Tuple, Union
+from typing import Tuple, Union
 
 import numpy
-import xarray as xr
 from affine import Affine
 
 
-def spatial_dims(
-    xx: Union[xr.DataArray, xr.Dataset], relaxed: bool = False
-) -> Optional[Tuple[str, str]]:
-    """Find spatial dimensions of `xx`.
-
-    Checks for presence of dimensions named:
-      y, x | latitude, longitude | lat, lon
-
-    Returns
-    =======
-    None -- if no dimensions with expected names are found
-    ('y', 'x') | ('latitude', 'longitude') | ('lat', 'lon')
-
-    If *relaxed* is True and none of the above dimension names are found,
-    assume that last two dimensions are spatial dimensions.
-    """
-    guesses = [("y", "x"), ("latitude", "longitude"), ("lat", "lon")]
-
-    _dims = [str(dim) for dim in xx.dims]
-    dims = set(_dims)
-    for guess in guesses:
-        if dims.issuperset(guess):
-            return guess
-
-    if relaxed and len(_dims) >= 2:
-        return _dims[-2], _dims[-1]
-
-    return None
-
-
 def maybe_zero(x: float, tol: float) -> float:
     """Turn almost zeros to actual zeros"""
     if abs(x) < tol:
diff --git a/odc/geo/xarray_geoextensions.py b/odc/geo/xarray_geoextensions.py
index 5a427a6..7447d23 100644
--- a/odc/geo/xarray_geoextensions.py
+++ b/odc/geo/xarray_geoextensions.py
@@ -7,14 +7,13 @@ Add geometric extensions to :class:`xarray.Dataset` and :class:`xarray.DataArray
 with Data Cube by Monkey Patching those classes.
 """
 import warnings
+from typing import Optional, Tuple, Union
 
 import xarray
 
-from ._base import CRS, CRSError, GeoBox
-from .math import affine_from_axis, spatial_dims
-
-# from datacube.utils import geometry, spatial_dims
-# from datacube.utils.math import affine_from_axis
+from ._crs import CRS, CRSError
+from ._geobox import GeoBox
+from .math import affine_from_axis
 
 
 def _norm_crs(crs):
@@ -130,6 +129,36 @@ def _get_crs_from_coord(obj, mode="strict"):
     raise ValueError(f"Mode needs to be: strict|any|all got {mode}")
 
 
+def spatial_dims(
+    xx: Union[xarray.DataArray, xarray.Dataset], relaxed: bool = False
+) -> Optional[Tuple[str, str]]:
+    """Find spatial dimensions of `xx`.
+
+    Checks for presence of dimensions named:
+      y, x | latitude, longitude | lat, lon
+
+    Returns
+    =======
+    None -- if no dimensions with expected names are found
+    ('y', 'x') | ('latitude', 'longitude') | ('lat', 'lon')
+
+    If *relaxed* is True and none of the above dimension names are found,
+    assume that last two dimensions are spatial dimensions.
+    """
+    guesses = [("y", "x"), ("latitude", "longitude"), ("lat", "lon")]
+
+    _dims = [str(dim) for dim in xx.dims]
+    dims = set(_dims)
+    for guess in guesses:
+        if dims.issuperset(guess):
+            return guess
+
+    if relaxed and len(_dims) >= 2:
+        return _dims[-2], _dims[-1]
+
+    return None
+
+
 def _xarray_affine_impl(obj):
     sdims = spatial_dims(obj, relaxed=True)
     if sdims is None:


Test patch:
diff --git a/tests/test_geometry.py b/tests/test_geometry.py
index ce69030..92d94fb 100644
--- a/tests/test_geometry.py
+++ b/tests/test_geometry.py
@@ -39,18 +39,16 @@ from odc.geo import (
     split_translation,
     w_,
 )
-from odc.geo._base import (
+from odc.geo._crs import _guess_crs_str, norm_crs, norm_crs_or_error
+from odc.geo._geobox import (
     _align_pix,
-    _guess_crs_str,
     _mk_crs_coord,
-    _norm_crs,
-    _norm_crs_or_error,
     _round_to_res,
     bounding_box_in_pixel_domain,
-    force_2d,
     geobox_intersection_conservative,
     geobox_union_conservative,
 )
+from odc.geo._geom import densify, force_2d
 from odc.geo.testutils.geom import (
     SAMPLE_WKT_WITHOUT_AUTHORITY,
     AlbersGS,
@@ -489,8 +487,6 @@ def test_boundingbox():
 
 
 def test_densify():
-    from odc.geo._base import densify
-
     s_x10 = [(0, 0), (10, 0)]
     assert densify(s_x10, 20) == s_x10
     assert densify(s_x10, 200) == s_x10
@@ -1599,10 +1595,10 @@ def test_base_internals():
     assert _round_to_res(0.0, 1.0) == 0
     assert _round_to_res(0.05, 1.0) == 0
 
-    assert _norm_crs(None) is None
+    assert norm_crs(None) is None
 
     with pytest.raises(ValueError):
-        _norm_crs_or_error(None)
+        norm_crs_or_error(None)
 
 
 def test_geom_clone():