Repository: opendatacube/odc-geo (python)

Issue:
Output resolution arguments are dropped in `[Dataset].odc.reproject(...)`
I have a 1 metre LiDAR dataset that I am loading from a raster file via `rioxarray`. I would like to reproject this dataset into a specific CRS (Australian Albers) and resolution (10 metres), similar to using `output_crs` and `resolution` in `datacube.load()`. 

I can easily reproject to the CRS like this:
```
ds.odc.reproject(how="EPSG:3577")
```
However, specifying `resolution` has no effect:
```
ds.odc.reproject(how="EPSG:3577", resolution=10)
```

Although one workaround would be to manually create a new GeoBox with my custom CRS and resolution, I would love to be able to simply set both CRS and resolution directly within `.odc.reproject()` for ease of use - especially since in this example I don't even need a particularly precise output grid (e.g. snapping/alignment), just a quick and dirty transformation to a specific CRS and lower resolution so I can more easily process a very large LIDAR dataset). 

Full code example:
```
import rioxarray
import odc.geo.xr

ds = rioxarray.open_rasterio(
    filename="https://dotazprdauegisextpubst01.blob.core.windows.net/transport-wa-public/bathymetry/rasters/HD20180814_Laser.bag",
    band_as_variable=True,
    chunks={},
)

ds_albers = ds.odc.reproject(how="EPSG:3577")
ds_albers_10m = ds.odc.reproject(how="EPSG:3577", resolution=10)

print(ds_albers.odc.geobox.resolution)
print(ds_albers_10m.odc.geobox.resolution)
```


Pull request description:
Unify reproject behaviour on Dataset and DataArray
Expose all the possible `GeoBox.from_bbox` parameters on `.reproject`, `.output_geobox`,  `.to_crs`

- Closes #122

Patch:
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 702543c..2720365 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -306,7 +306,7 @@ jobs:
       - name: Deploy to Netlify
         id: netlify
         if: github.event_name == 'pull_request'
-        uses: nwtgck/actions-netlify@v1.2
+        uses: nwtgck/actions-netlify@v2
         with:
           production-branch: "main"
           publish-dir: "docs/_build/html"
diff --git a/docs/api-geobox.rst b/docs/api-geobox.rst
index fbf675c..8073e65 100644
--- a/docs/api-geobox.rst
+++ b/docs/api-geobox.rst
@@ -6,7 +6,6 @@ GeoBox
    :toctree: _api/
 
    GeoBox
-   AnchorEnum
    GeoBox.from_geopolygon
    GeoBox.from_bbox
    GeoBox.affine
diff --git a/docs/api.rst b/docs/api.rst
index cc7a55c..60d5425 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -175,6 +175,7 @@ Basic types used for unambigous specification of X/Y values.
    Resolution
    Shape2d
    Index2d
+   AnchorEnum
 
    xy_
    yx_
diff --git a/odc/geo/__init__.py b/odc/geo/__init__.py
index fc65e75..4aa865c 100644
--- a/odc/geo/__init__.py
+++ b/odc/geo/__init__.py
@@ -11,6 +11,7 @@
 from ._version import __version__
 
 from .types import (
+    AnchorEnum,
     XY,
     Index2d,
     Shape2d,
@@ -46,6 +47,7 @@ from .geom import (
 
 __all__ = [
     "__version__",
+    "AnchorEnum",
     "XY",
     "Index2d",
     "Shape2d",
diff --git a/odc/geo/_xr_interop.py b/odc/geo/_xr_interop.py
index e714058..fa40da6 100644
--- a/odc/geo/_xr_interop.py
+++ b/odc/geo/_xr_interop.py
@@ -15,6 +15,7 @@ from typing import (
     Dict,
     Hashable,
     List,
+    Literal,
     Optional,
     Set,
     Tuple,
@@ -30,12 +31,12 @@ from ._interop import have, is_dask_collection
 from ._rgba import colorize, to_rgba
 from .crs import CRS, CRSError, SomeCRS, norm_crs_or_error
 from .gcp import GCPGeoBox, GCPMapping
-from .geobox import Coordinate, GeoBox
+from .geobox import Coordinate, GeoBox, GeoboxAnchor
 from .geom import Geometry
 from .math import affine_from_axis, maybe_int, resolution_from_affine
 from .overlap import compute_output_geobox
 from .roi import roi_is_empty
-from .types import Resolution, xy_
+from .types import Resolution, SomeResolution, SomeShape, xy_
 
 # pylint: disable=import-outside-toplevel
 # pylint: disable=too-many-lines
@@ -275,7 +276,7 @@ def mask(
         A :py:class:`~xarray.Dataset` or :py:class:`~xarray.DataArray`
         masked by ``poly``.
 
-    ..seealso:: :py:meth:`odc.geo.xr.rasterize`
+    .. seealso:: :py:meth:`odc.geo.xr.rasterize`
     """
     # Rasterise `poly` into geobox of `xx`
     rasterized = rasterize(
@@ -325,7 +326,7 @@ def crop(
         A :py:class:`~xarray.Dataset` or :py:class:`~xarray.DataArray`
         cropped and optionally masked to the spatial extent of ``poly``.
 
-    ..seealso:: :py:meth:`odc.geo.xr.mask`
+    .. seealso:: :py:meth:`odc.geo.xr.mask`
     """
     # Create new geobox with pixel grid of `xx` but enclosing `poly`.
     poly_geobox = xx.odc.geobox.enclosing(poly)
@@ -582,13 +583,78 @@ def xr_reproject(
     *,
     resampling: Union[str, int] = "nearest",
     dst_nodata: Optional[float] = None,
+    resolution: Union[SomeResolution, Literal["auto", "fit", "same"]] = "auto",
+    shape: Union[SomeShape, int, None] = None,
+    tight: bool = False,
+    anchor: GeoboxAnchor = "default",
+    tol: float = 0.01,
+    round_resolution: Union[None, bool, Callable[[float, str], float]] = None,
     **kw,
 ) -> XrT:
     """
     Reproject raster to different projection/resolution.
 
+    :param src:
+      :py:class:`~xarray.Dataset` or :py:class:`~xarray.DataArray` to reproject.
+
+    :param how:
+      How to reproject the raster. Can be a GeoBox or a CRS (e.g. CRS object or
+      an "ESPG:XXXX" string/integer). If a CRS is provided, the output pixel
+      grid can be customised further via ``resolution``, ``shape``, ``tight``,
+      ``anchor``, ``tol``, ``round_resolution``.
+
+    :param resampling:
+      Resampling method to use when reprojecting the raster. Defaults to
+      "nearest", also supports "average", "bilinear", "cubic", "cubic_spline",
+      "lanczos", "mode", "gauss", "max", "min", "med", "q1", "q3".
+
+    :param dst_nodata:
+      Set a custom nodata value for the output resampled raster.
+
+    :param resolution:
+
+       * "same" use exactly the same resolution as src
+       * "fit" use center pixel to determine scale change between the two
+       * | "auto" is to use the same resolution on the output if CRS units are
+         | the same between the source and destination and otherwise use "fit"
+       * Ignored if ``shape=`` is supplied
+       * Else resolution in the units of the output crs
+
+    :param shape:
+      Span that many pixels, if it's a single number then span that many pixels
+      along the longest dimension, other dimension will be computed to maintain
+      roughly square pixels. Takes precedence over ``resolution=`` parameter.
+
+    :param tight:
+      By default output pixel grid is adjusted to align pixel edges to X/Y axis,
+      suppling ``tight=True`` produces unaligned geobox on the output.
+
+    :param anchor:
+      Control pixel snapping, default is to snap pixel edge to ``X=0,Y=0``.
+      Ignored when ``tight=True`` is supplied.
+
+    :param tol:
+       Fraction of the output pixel that can be ignored, defaults to 1/100.
+       Bounding box of the output geobox is allowed to be smaller by that amount
+       than transformed footprint of the original.
+
+    :param round_resolution:
+      ``round_resolution(res: float, units: str) -> float``
+
     This method uses :py:mod:`rasterio`.
+
+    .. seealso:: :py:meth:`odc.geo.overlap.compute_output_geobox`
+
     """
+    kw = {
+        "shape": shape,
+        "resolution": resolution,
+        "tight": tight,
+        "anchor": anchor,
+        "tol": tol,
+        "round_resolution": round_resolution,
+        **kw,
+    }
     if isinstance(src, xarray.DataArray):
         return _xr_reproject_da(
             src, how, resampling=resampling, dst_nodata=dst_nodata, **kw
@@ -598,6 +664,15 @@ def xr_reproject(
     )
 
 
+def _extract_output_geobox_params(kw):
+    # NOTE: modifies input, removes keys
+    out = {}
+    for k in ("tight", "anchor", "resolution", "shape", "tol", "round_resolution"):
+        if k in kw:
+            out[k] = kw.pop(k)
+    return out
+
+
 def _xr_reproject_ds(
     src: Any,
     how: Union[SomeCRS, GeoBox],
@@ -615,10 +690,12 @@ def _xr_reproject_ds(
     if src.odc.geobox is None:
         raise ValueError("Can not reproject non-georegistered array.")
 
+    kw_gbox = _extract_output_geobox_params(kw)
+
     if isinstance(how, GeoBox):
         dst_geobox = how
     else:
-        dst_geobox = src.odc.output_geobox(how)
+        dst_geobox = src.odc.output_geobox(how, **kw_gbox)
 
     def _maybe_reproject(dv: xarray.DataArray):
         if dv.odc.geobox is None:
@@ -642,11 +719,6 @@ def _xr_reproject_da(
     dst_nodata: Optional[float] = None,
     **kw,
 ) -> xarray.DataArray:
-    """
-    Reproject raster to different projection/resolution.
-
-    This method uses :py:mod:`rasterio`.
-    """
     # pylint: disable=too-many-locals
     assert isinstance(src, xarray.DataArray)
 
@@ -659,13 +731,12 @@ def _xr_reproject_da(
     if src_gbox is None or src_gbox.crs is None:
         raise ValueError("Can not reproject non-georegistered array.")
 
-    resolution = kw.pop("resolution", "auto")
-    tight = kw.pop("tight", False)
+    kw_gbox = _extract_output_geobox_params(kw)
 
     if isinstance(how, GeoBox):
         dst_geobox = how
     else:
-        dst_geobox = src.odc.output_geobox(how, resolution=resolution, tight=tight)
+        dst_geobox = src.odc.output_geobox(how, **kw_gbox)
 
     # compute destination shape by replacing spatial dimensions shape
     ydim = src.odc.ydim
@@ -776,7 +847,7 @@ class ODCExtension:
         """
         Compute geobox of this data in other projection.
 
-        ..seealso:: :py:meth:`odc.geo.overlap.compute_output_geobox`
+        .. seealso:: :py:meth:`odc.geo.overlap.compute_output_geobox`
         """
         gbox = self.geobox
         if gbox is None:
@@ -797,6 +868,7 @@ class ODCExtension:
 
     if have.rasterio:
         explore = _wrap_op(explore)
+        reproject = _wrap_op(xr_reproject)
 
 
 @xarray.register_dataarray_accessor("odc")
@@ -850,7 +922,6 @@ class ODCExtensionDa(ODCExtension):
         write_cog = _wrap_op(write_cog)
         to_cog = _wrap_op(to_cog)
         compress = _wrap_op(compress)
-        reproject = _wrap_op(_xr_reproject_da)
         add_to = _wrap_op(add_to)
 
 
@@ -882,9 +953,6 @@ class ODCExtensionDs(ODCExtension):
     ) -> xarray.DataArray:
         return to_rgba(self._xx, bands=bands, vmin=vmin, vmax=vmax)
 
-    if have.rasterio:
-        reproject = _wrap_op(_xr_reproject_ds)
-
 
 ODCExtensionDs.to_rgba.__doc__ = to_rgba.__doc__
 
diff --git a/odc/geo/geobox.py b/odc/geo/geobox.py
index b51684c..6255235 100644
--- a/odc/geo/geobox.py
+++ b/odc/geo/geobox.py
@@ -6,7 +6,6 @@ import importlib
 import itertools
 import math
 from collections import OrderedDict, namedtuple
-from enum import Enum
 from typing import (
     Callable,
     Dict,
@@ -48,6 +47,7 @@ from .roi import (
 from .types import (
     ROI,
     XY,
+    AnchorEnum,
     Chunks2d,
     MaybeInt,
     NormalizedROI,
@@ -64,26 +64,40 @@ from .types import (
     xy_,
 )
 
+GeoboxAnchor = Union[
+    AnchorEnum,
+    XY[float],
+    float,
+    Literal["center"],
+    Literal["centre"],
+    Literal["edge"],
+    Literal["floating"],
+    Literal["default"],
+]
 
-class AnchorEnum(Enum):
-    """
-    Defines which way to snap geobox pixel grid.
-    """
-
-    EDGE = 0
-    """Snap pixel edges to multiples of pixel size."""
-
-    CENTER = 1
-    """Snap pixel centers to multiples of pixel size."""
-
-    FLOATING = 2
-    """Turn off pixel snapping."""
+# pylint: disable=invalid-name,too-many-public-methods,too-many-lines
+Coordinate = namedtuple("Coordinate", ("values", "units", "resolution"))
 
 
-GeoboxAnchor = Union[AnchorEnum, XY[float]]
+def _norm_anchor(anchor: GeoboxAnchor) -> Union[AnchorEnum, XY[float]]:
+    if isinstance(anchor, AnchorEnum):
+        return anchor
+    if isinstance(anchor, (float, int)):
+        if anchor == 0:
+            return AnchorEnum.EDGE
+        if anchor == 0.5:
+            return AnchorEnum.CENTER
+        return xy_(float(anchor), float(anchor))
+    if isinstance(anchor, XY):
+        return anchor
 
-# pylint: disable=invalid-name,too-many-public-methods,too-many-lines
-Coordinate = namedtuple("Coordinate", ("values", "units", "resolution"))
+    return {
+        "center": AnchorEnum.CENTER,
+        "centre": AnchorEnum.CENTER,
+        "edge": AnchorEnum.EDGE,
+        "floating": AnchorEnum.FLOATING,
+        "default": AnchorEnum.EDGE,
+    }[anchor]
 
 
 class GeoBoxBase:
@@ -435,7 +449,7 @@ class GeoBox(GeoBoxBase):
         tight: bool = False,
         shape: Union[SomeShape, int, None] = None,
         resolution: Optional[SomeResolution] = None,
-        anchor: GeoboxAnchor = AnchorEnum.EDGE,
+        anchor: GeoboxAnchor = "default",
         tol: float = 0.01,
     ) -> "GeoBox":
         """
@@ -460,6 +474,7 @@ class GeoBox(GeoBoxBase):
         """
         # pylint: disable=too-many-locals, too-many-branches
 
+        anchor = _norm_anchor(anchor)
         _snap: Optional[XY[float]] = None
 
         if tight:
@@ -531,7 +546,7 @@ class GeoBox(GeoBoxBase):
         *,
         shape: Union[SomeShape, int, None] = None,
         tight: bool = False,
-        anchor: GeoboxAnchor = AnchorEnum.EDGE,
+        anchor: GeoboxAnchor = "default",
         tol: float = 0.01,
     ) -> "GeoBox":
         """
@@ -733,7 +748,10 @@ class GeoBox(GeoBoxBase):
         crs: SomeCRS,
         *,
         resolution: Literal["auto", "fit", "same"] = "auto",
+        shape: Union[SomeShape, int, None] = None,
         tight: bool = False,
+        anchor: GeoboxAnchor = "default",
+        tol: float = 0.01,
         round_resolution: Union[None, bool, Callable[[float, str], float]] = None,
     ) -> "GeoBox":
         """
@@ -749,15 +767,32 @@ class GeoBox(GeoBoxBase):
            * | "auto" is to use the same resolution on the output if CRS units are the same
              |  between the source and destination and otherwise use "fit"
 
+        :param shape:
+          Span that many pixels, if it's a single number then span that many pixels
+          along the longest dimension, other dimension will be computed to maintain
+          roughly square pixels. Takes precedence over ``resolution=`` parameter.
+
         :param tight:
-          By default output pixel grid is adjusted to align pixel edges to X/Y axis, suppling
-          ``tight=True`` produces unaligned geobox on the output.
+          By default output pixel grid is adjusted to align pixel edges to X/Y
+          axis, suppling ``tight=True`` produces unaligned geobox on the output.
+
+        :param anchor:
+            Control pixel snapping, default is to snap pixel edge to
+            ``X=0,Y=0``. Ignored when ``tight=True`` is supplied.
+
+        :param tol:
+            Fraction of the output pixel that can be ignored, defaults to 1/100.
+            Bounding box of the output geobox is allowed to be smaller by that
+            amount than transformed footprint of the original.
+
+        :param round_resolution:
+          ``round_resolution(res: float, units: str) -> float``
 
         :return:
-           Similar resolution, axis aligned geobox that fully encloses this one but in a different
-           projection.
+           Similar resolution, axis aligned geobox that fully encloses this one
+           but in a different projection.
         """
-        # pylint: disable=import-outside-toplevel
+        # pylint: disable=import-outside-toplevel,too-many-arguments
         # can't be up-top due to circular imports issues
         from .overlap import compute_output_geobox
 
@@ -765,8 +800,11 @@ class GeoBox(GeoBoxBase):
             self,
             crs,
             resolution=resolution,
+            shape=shape,
             tight=tight,
+            anchor=anchor,
             round_resolution=round_resolution,
+            tol=tol,
         )
 
     def __str__(self):
diff --git a/odc/geo/overlap.py b/odc/geo/overlap.py
index 7230141..bf7d233 100644
--- a/odc/geo/overlap.py
+++ b/odc/geo/overlap.py
@@ -11,7 +11,7 @@ from affine import Affine
 
 from .crs import SomeCRS
 from .gcp import GCPGeoBox
-from .geobox import GeoBox, GeoBoxBase, gbox_boundary
+from .geobox import GeoBox, GeoboxAnchor, GeoBoxBase, gbox_boundary
 from .math import (
     affine_from_pts,
     decompose_rws,
@@ -563,14 +563,17 @@ def compute_output_geobox(
     crs: SomeCRS,
     *,
     resolution: Union[SomeResolution, Literal["auto", "fit", "same"]] = "auto",
+    shape: Union[SomeShape, int, None] = None,
     tight: bool = False,
+    anchor: GeoboxAnchor = "default",
+    tol: float = 0.01,
     round_resolution: Union[None, bool, Callable[[float, str], float]] = None,
 ) -> GeoBox:
     """
     Compute output ``GeoBox``.
 
-    Find best fitting, axis aligned GeoBox in a different coordinate reference given source
-    ``GeoBox`` on input.
+    Find best fitting, axis aligned GeoBox in a different coordinate reference
+    given source ``GeoBox`` on input.
 
     :param gbox:
        Source geobox.
@@ -582,24 +585,41 @@ def compute_output_geobox(
 
        * "same" use exactly the same resolution as src
        * "fit" use center pixel to determine scale change between the two
-       * | "auto" is to use the same resolution on the output if CRS units are the same
-         |  between the source and destination and otherwise use "fit"
+       * | "auto" is to use the same resolution on the output if CRS units are
+         | the same between the source and destination and otherwise use "fit"
+       * Ignored if ``shape=`` is supplied
        * Else resolution in the units of the output crs
 
+    :param shape:
+      Span that many pixels, if it's a single number then span that many pixels
+      along the longest dimension, other dimension will be computed to maintain
+      roughly square pixels. Takes precedence over ``resolution=`` parameter.
+
     :param tight:
-      By default output pixel grid is adjusted to align pixel edges to X/Y axis, suppling
-      ``tight=True`` produces unaligned geobox on the output.
+      By default output pixel grid is adjusted to align pixel edges to X/Y axis,
+      suppling ``tight=True`` produces unaligned geobox on the output.
+
+    :param anchor:
+      Control pixel snapping, default is to snap pixel edge to ``X=0,Y=0``.
+      Ignored when ``tight=True`` is supplied.
+
+    :param tol:
+       Fraction of the output pixel that can be ignored, defaults to 1/100.
+       Bounding box of the output geobox is allowed to be smaller by that amount
+       than transformed footprint of the original.
 
-    :param round_resolution: ``round_resolution(res: float, units: str) -> float``
+    :param round_resolution:
+      ``round_resolution(res: float, units: str) -> float``
 
     :return:
-       Similar resolution, axis aligned geobox that fully encloses source one but in a different
-       projection.
+       Similar resolution, axis aligned geobox that fully encloses source one
+       but in a different projection.
     """
-    # pylint: disable=too-many-locals
+    # pylint: disable=too-many-locals,too-many-arguments
     src_crs = gbox.crs
     assert src_crs is not None
 
+    # figure out "true" dts_crs (handles "utm" -> actual CRS)
     bbox = gbox.footprint(crs, buffer=0.9, npoints=100).boundingbox
     dst_crs = bbox.crs
     assert dst_crs is not None
@@ -607,13 +627,17 @@ def compute_output_geobox(
     if (
         dst_crs == src_crs
         and resolution in ("auto", "same")
+        and shape is None
+        and anchor == "default"
         and isinstance(gbox, GeoBox)
     ):
         return gbox
 
     same_units = src_crs.units == dst_crs.units
 
-    if resolution == "same":
+    if shape is not None:
+        res: Optional[SomeResolution] = None
+    elif resolution == "same":
         res = gbox.resolution
     elif resolution == "auto" and same_units:
         res = gbox.resolution
@@ -648,4 +672,12 @@ def compute_output_geobox(
 
         res = res_(resolution)
 
-    return GeoBox.from_bbox(bbox, dst_crs, resolution=res, tight=tight)
+    return GeoBox.from_bbox(
+        bbox,
+        dst_crs,
+        shape=shape,
+        resolution=res,
+        tight=tight,
+        anchor=anchor,
+        tol=tol,
+    )
diff --git a/odc/geo/types.py b/odc/geo/types.py
index 455767c..92b73ca 100644
--- a/odc/geo/types.py
+++ b/odc/geo/types.py
@@ -1,4 +1,5 @@
 """Basic types."""
+from enum import Enum
 from typing import (
     Callable,
     Generic,
@@ -458,6 +459,21 @@ OutlineMode = Union[
 ]
 
 
+class AnchorEnum(Enum):
+    """
+    Defines which way to snap geobox pixel grid.
+    """
+
+    EDGE = 0
+    """Snap pixel edges to multiples of pixel size."""
+
+    CENTER = 1
+    """Snap pixel centers to multiples of pixel size."""
+
+    FLOATING = 2
+    """Turn off pixel snapping."""
+
+
 class _Func2Map(Mapping[TK, TV]):
     """
     Turn ``f(K)`` into ``ff[K]``.


Test patch:
diff --git a/tests/test_geobox.py b/tests/test_geobox.py
index 052bc44..ce4b70c 100644
--- a/tests/test_geobox.py
+++ b/tests/test_geobox.py
@@ -412,14 +412,46 @@ def test_from_bbox():
         bbox, shape=shape, anchor=AnchorEnum.FLOATING
     ) == GeoBox.from_bbox(bbox, shape=shape, tight=True)
 
+    assert GeoBox.from_bbox(
+        bbox, shape=shape, anchor=AnchorEnum.FLOATING
+    ) == GeoBox.from_bbox(bbox, shape=shape, anchor="floating")
+
     assert GeoBox.from_bbox(
         bbox, shape=shape, anchor=AnchorEnum.EDGE
     ) == GeoBox.from_bbox(bbox, shape=shape, anchor=xy_(0, 0))
 
+    assert GeoBox.from_bbox(
+        bbox, shape=shape, anchor=AnchorEnum.EDGE
+    ) != GeoBox.from_bbox(bbox, shape=shape, anchor=xy_(0.3, 0))
+
+    assert GeoBox.from_bbox(bbox, shape=shape, anchor=0.3) == GeoBox.from_bbox(
+        bbox, shape=shape, anchor=xy_(0.3, 0.3)
+    )
+
+    assert GeoBox.from_bbox(
+        bbox, shape=shape, anchor=AnchorEnum.EDGE
+    ) == GeoBox.from_bbox(bbox, shape=shape, anchor=0)
+
+    assert GeoBox.from_bbox(
+        bbox, shape=shape, anchor=AnchorEnum.EDGE
+    ) == GeoBox.from_bbox(bbox, shape=shape)
+
     assert GeoBox.from_bbox(
         bbox, shape=shape, anchor=AnchorEnum.CENTER
     ) == GeoBox.from_bbox(bbox, shape=shape, anchor=xy_(0.5, 0.5))
 
+    assert GeoBox.from_bbox(
+        bbox, shape=shape, anchor=AnchorEnum.CENTER
+    ) == GeoBox.from_bbox(bbox, shape=shape, anchor=0.5)
+
+    assert GeoBox.from_bbox(
+        bbox, shape=shape, anchor=AnchorEnum.CENTER
+    ) == GeoBox.from_bbox(bbox, shape=shape, anchor="center")
+
+    assert GeoBox.from_bbox(
+        bbox, shape=shape, anchor=AnchorEnum.CENTER
+    ) == GeoBox.from_bbox(bbox, shape=shape, anchor="centre")
+
     assert GeoBox.from_bbox((1, 40, 2, 43), "utm", shape=(12, 12)).crs == "epsg:32631"
 
     # one of resolution= or shape= must be supplied
diff --git a/tests/test_overlap.py b/tests/test_overlap.py
index 1a53467..ff89548 100644
--- a/tests/test_overlap.py
+++ b/tests/test_overlap.py
@@ -9,7 +9,7 @@ import numpy as np
 import pytest
 from affine import Affine
 
-from odc.geo import CRS, geom, res_, resyx_, wh_, xy_
+from odc.geo import CRS, AnchorEnum, geom, res_, resyx_, wh_, xy_
 from odc.geo.geobox import GeoBox, scaled_down_geobox
 from odc.geo.gridspec import GridSpec
 from odc.geo.math import affine_from_pts, decompose_rws, is_affine_st, stack_xy
@@ -487,6 +487,15 @@ def test_compute_output_geobox():
     npix_change = (src.shape[0] * src.shape[1]) / (dst.shape[0] * dst.shape[1])
     assert 0.8 < npix_change < 1.1
 
+    # check anchor= having effect
+    assert compute_output_geobox(
+        src, 4326, anchor=AnchorEnum.CENTER
+    ) != compute_output_geobox(src, 4326, anchor=AnchorEnum.EDGE)
+
+    # check shape= having effect
+    assert compute_output_geobox(src, 4326, shape=(10, 20)).shape == (10, 20)
+    assert max(compute_output_geobox(src, 4326, shape=100, tight=True).shape) == 100
+
     # go back from 4326
     _src = dst
     dst = compute_output_geobox(_src, src.crs)