Source code for odc.geo.geobox

# 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
from __future__ import annotations

import importlib
import itertools
import math
from collections import OrderedDict, namedtuple
from typing import (
    Any,
    Callable,
    Dict,
    Iterator,
    List,
    Literal,
    Mapping,
    Optional,
    Sequence,
    Tuple,
    Union,
)

import numpy
from affine import Affine

from . import geom
from .crs import CRS, MaybeCRS, SomeCRS, norm_crs
from .geom import BoundingBox, Geometry, bbox_intersection, bbox_union
from .math import (
    clamp,
    extract_anchor,
    is_affine_st,
    is_almost_int,
    maybe_zero,
    resolution_from_affine,
    snap_affine,
    snap_grid,
    split_translation,
)
from .roi import (
    RoiTiles,
    align_up,
    clip_tiles,
    roi_boundary,
    roi_normalise,
    roi_shape,
    roi_tiles,
)
from .types import (
    ROI,
    XY,
    AnchorEnum,
    Chunks2d,
    MaybeInt,
    NormalizedROI,
    OutlineMode,
    Resolution,
    Shape2d,
    SomeIndex2d,
    SomeResolution,
    SomeShape,
    Unset,
    func2map,
    res_,
    shape_,
    xy_,
)

GeoboxAnchor = Union[
    AnchorEnum,
    XY[float],
    float,
    Literal["center"],
    Literal["centre"],
    Literal["edge"],
    Literal["floating"],
    Literal["default"],
]

# pylint: disable=invalid-name,too-many-public-methods,too-many-lines
Coordinate = namedtuple("Coordinate", ("values", "units", "resolution"))


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

    return {
        "center": AnchorEnum.CENTER,
        "centre": AnchorEnum.CENTER,
        "edge": AnchorEnum.EDGE,
        "floating": AnchorEnum.FLOATING,
        "default": AnchorEnum.EDGE,
    }[anchor]


class GeoBoxBase:
    """
    Defines the location and resolution of a rectangular grid of data,
    including it's :py:class:`~odc.geo.crs.CRS`.

    :param shape: Shape in pixels ``(ny, nx)``
    :param crs: Coordinate Reference System
    :param affine: Affine transformation defining the location of the geobox
    """

    __slots__ = ("_shape", "_affine", "_crs", "_extent", "_lazy_ui")

    def __init__(self, shape: SomeShape, affine: Affine, crs: MaybeCRS):
        shape = shape_(shape)

        self._shape = shape
        self._affine = affine
        self._crs = norm_crs(crs)
        self._extent: Optional[Geometry] = None
        self._lazy_ui = None

        if self.shape.x < 0 or self.shape.y < 0:
            raise ValueError(
                f"Got shape {self._shape.yx!r}: negative sizes are not allowed."
            )

    @property
    def width(self) -> int:
        """Width in pixels (nx)."""
        return self._shape.x

    @property
    def height(self) -> int:
        """Height in pixels (ny)."""
        return self._shape.y

    @property
    def shape(self) -> Shape2d:
        """Shape in pixels ``(height, width)``."""
        return self._shape

    @property
    def aspect(self) -> float:
        """Aspect ratio (X/Y in pixel space)."""
        return self._shape.aspect

    @property
    def crs(self) -> Optional[CRS]:
        """Coordinate Reference System of the GeoBox."""
        return self._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

    dims = dimensions

    def is_empty(self) -> bool:
        """Check if geobox is "empty"."""
        return 0 in self._shape

    def __bool__(self) -> bool:
        return not self.is_empty()

    @property
    def resolution(self) -> Resolution:
        """Resolution, pixel size in CRS units."""
        return resolution_from_affine(self._affine)

    def boundary(self, pts_per_side: int = 16) -> numpy.ndarray:
        """
        Boundary of a :py:class:`~odc.geo.geobox.GeoBox`.

        Construct a ring of points in pixel space along the edge of the geobox.

        :param pts_per_side: Number of points per side, default is 16.

        :return:
          Points in pixel space along the perimeter of a GeoBox as a ``Nx2`` array
          in pixel coordinates.
        """
        ny, nx = self._shape.yx
        return roi_boundary(numpy.s_[0:ny, 0:nx], pts_per_side)

    @property
    def alignment(self) -> XY[float]:
        """
        Alignment of pixel boundaries in CRS units.

        This is usally ``(0,0)``.

        .. seealso:: :py:meth:`~odc.geo.geobox.GeoBoxBase.anchor`
        """
        # TODO: deprecate alignment method
        rx, _, tx, _, ry, ty, *_ = self._affine
        return xy_(tx % abs(rx), ty % abs(ry))

    @property
    def anchor(self) -> XY[float] | AnchorEnum:
        """
        Anchor point in pixel coordinates.

        This is only valid for linear sources.
        """
        if not self.linear:
            return AnchorEnum.FLOATING
        return extract_anchor(self._affine)

    @property
    def linear(self) -> bool:
        return True

    def wld2pix(self, x, y):
        return (~self._affine) * (x, y)

    def pix2wld(self, x, y):
        return self._affine * (x, y)

    @property
    def extent(self) -> Geometry:
        """GeoBox footprint in native CRS."""
        if self._extent is not None:
            return self._extent

        if self.linear:
            _extent = geom.polygon_from_transform(self._shape, self._affine, self._crs)
        else:
            _extent = geom.polygon(self.boundary(16).tolist(), self._crs).transform(
                self.pix2wld
            )

        self._extent = _extent
        return _extent

    @property
    def boundingbox(self) -> BoundingBox:
        """GeoBox bounding box in the native CRS."""
        return BoundingBox.from_transform(self._shape, self._affine, crs=self._crs)

    def _reproject_resolution(self, npoints: int = 100):
        bbox = self.extent.boundingbox
        span = max(bbox.span_x, bbox.span_y)
        return span / npoints

    def footprint(
        self, crs: SomeCRS, buffer: float = 0, npoints: int = 100
    ) -> Geometry:
        """
        Compute footprint in foreign CRS.

        :param crs: CRS of the destination
        :param buffer: amount to buffer in source pixels before transforming
        :param npoints: number of points per-side to use, higher number
                        is slower but more accurate
        """
        assert self.crs is not None
        ext = self.extent
        if buffer != 0:
            buffer = buffer * max(*self.resolution.xy)
            ext = ext.buffer(buffer)

        return ext.to_crs(crs, resolution=self._reproject_resolution(npoints)).dropna()

    @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.footprint("epsg:4326")

    @property
    def _ui(self):
        # pylint: disable=import-outside-toplevel
        from .ui import PixelGridDisplay

        if self._lazy_ui is not None:
            return self._lazy_ui

        gsd = max(*self.resolution.map(abs).xy)
        self._lazy_ui = PixelGridDisplay(self, self.pix2wld, gsd)
        return self._lazy_ui

    def svg(
        self,
        scale_factor: float = 1.0,
        mode: OutlineMode = "auto",
        notch: float = 0.0,
        grid_stroke: str = "pink",
    ) -> str:
        """
        Produce SVG paths.

        :param mode: One of pixel, native, geo (default is geo)
        :return: SVG path
        """
        return self._ui.svg(
            scale_factor=scale_factor, mode=mode, notch=notch, grid_stroke=grid_stroke
        )

    def grid_lines(self, step: int = 0, mode: OutlineMode = "native") -> Geometry:
        """
        Construct pixel edge aligned grid lines.
        """
        return self._ui.grid_lines(step=step, mode=mode)

    def outline(self, mode: OutlineMode = "native", notch: float = 0.1) -> Geometry:
        return self._ui.outline(mode, notch=notch)

    def _repr_svg_(self):
        # pylint: disable=protected-access
        return self._ui._render_svg()

    def _repr_html_(self):
        # pylint: disable=protected-access
        return self._ui._repr_html_()

    def compute_crop(self, roi) -> Tuple[Shape2d, Affine]:
        if isinstance(roi, BoundingBox):
            roi = roi.polygon
        if isinstance(roi, GeoBoxBase):
            roi = roi.extent

        if isinstance(roi, Geometry):
            if roi.crs is not None:
                roi = self.project(roi)
            pix_bbox = roi.boundingbox.round() & BoundingBox(
                0, 0, self.width, self.height
            )
            nx, ny = (max(1, int(span)) for span in (pix_bbox.span_x, pix_bbox.span_y))
            tx, ty = map(int, pix_bbox.bbox[:2])
            roi = numpy.s_[ty : ty + ny, tx : tx + nx]

        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")

        roi = roi_normalise(roi, self._shape.shape)

        if not all(s.step is None or s.step == 1 for s in roi):
            raise NotImplementedError("scaling not implemented, yet")

        ty, tx = (s.start for s in roi)
        ny, nx = roi_shape(roi)

        affine = self._affine * Affine.translation(tx, ty)
        return shape_((ny, nx)), affine

    def compute_zoom_out(self, factor: float) -> Tuple[Shape2d, Affine]:
        ny, nx = (max(1, math.ceil(s / factor)) for s in self.shape)
        A = self._affine * Affine.scale(factor, factor)
        return (shape_((ny, nx)), A)

    def compute_zoom_to(
        self,
        shape: Union[SomeShape, int, float, None] = None,
        *,
        resolution: Optional[SomeResolution] = None,
    ) -> Tuple[Shape2d, Affine]:
        """
        Change GeoBox shape.

        When supplied a single integer scale longest dimension to match that.

        :returns:
          GeoBox covering the same region but with different number of pixels and therefore resolution.
        """
        if shape is None:
            if resolution is None:
                raise ValueError("Have to supply shape or resolution")
            new_geobox = GeoBox.from_bbox(
                self.boundingbox, resolution=resolution, tight=True
            )
            return new_geobox.shape, new_geobox.affine

        if isinstance(shape, (int, float)):
            nmax = max(*self._shape)
            return self.compute_zoom_out(nmax / shape)

        shape = shape_(shape)
        sy, sx = (N / float(n) for N, n in zip(self._shape, shape.shape))
        A = self._affine * Affine.scale(sx, sy)
        return (shape, A)

    def project(self, g: Geometry) -> Geometry:
        """
        Map Geometry between world and pixel coords.

        When input geometry has no CRS, map from pixels to the world.

        When input geometry has CRS (can be different from GeoBox), project
        geometry into pixel coordinates, note that result is not clipped to the
        image bounds.
        """
        if g.crs is None:  # assume pixel plane
            g = g.transform(self.pix2wld)
            return Geometry(g.geom, self._crs)

        assert self._crs is not None
        if g.crs != self._crs:
            g = g.to_crs(self._crs)
        g = g.transform(self.wld2pix)
        return Geometry(g.geom, crs=None)

    def qr2sample(
        self,
        n: int,
        padding: Optional[float] = None,
        with_edges: bool = False,
        offset: int = 0,
    ) -> Geometry:
        """
        Generate quasi-random sample of image locations.

        :param n:
           Number of points

        :param padding:
           In pixels, minimal distance from the edge

        :param offset:
           Offset into quasi-random sequence from where to start

        :param edges:
           Also include samples along the edge and corners

        References: http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
        """
        nx, ny = self._shape.xy
        return BoundingBox(0, 0, nx, ny, None).qr2sample(
            n, padding=padding, with_edges=with_edges, offset=offset
        )

    def explore(
        self,
        map: Optional[Any] = None,
        grid_lines: bool = True,
        tiles: Any = "OpenStreetMap",
        attr: Optional[str] = None,
        map_kwds: Optional[Dict[str, Any]] = None,
        **kwargs: Any,
    ) -> Any:
        """
        Plot a visual representation of the GeoBox's extent on an
        interactive :py:mod:`folium` leaflet map.

        :param map:
            An optional existing :py:mod:`folium` map object to plot into.
            By default, a new map object will be created.
        :param: grid_lines:
            Whether to plot pixel edge aligned grid lines over the GeoBox
            extent.
        :param tiles:
            Map tileset to use for the map basemap. Supports any option
            supported by :py:mod:`folium`, including "OpenStreetMap",
            "CartoDB positron", "CartoDB dark_matter" or a custom XYZ URL.
        :param attr:
            Map tile attribution; only required if passing custom tile URL.
        :param map_kwds:
            Additional keyword arguments to pass to :py:class:`folium.Map`.
        :param kwargs:
            Additional keyword arguments to pass to :py:class:`folium.GeoJson`.

        :return: A :py:mod:`folium` map containing the plotted GeoBox.
        """
        # pylint: disable=redefined-builtin

        # Add outline to map
        map = self.outline(mode="geo", notch=0.0).explore(
            map=map, tiles=tiles, attr=attr, map_kwds=map_kwds, **kwargs
        )

        # Optionally overlay grid lines
        if grid_lines:
            self.grid_lines(mode="geo").explore(
                map=map,
                style_function=lambda feature: {
                    "weight": 1,
                    "fillOpacity": 0.3,
                },
            )
        return map

    def __getitem__(self, roi) -> "GeoBoxBase":
        raise NotImplementedError()


[docs] class GeoBox(GeoBoxBase): """ Defines the location and resolution of a rectangular grid of data, including it's :py:class:`~odc.geo.crs.CRS`. :param shape: Shape in pixels ``(ny, nx)`` :param crs: Coordinate Reference System :param affine: Affine transformation defining the location of the geobox """ __slots__ = ()
[docs] def __init__(self, shape: SomeShape, affine: Affine, crs: MaybeCRS): GeoBoxBase.__init__(self, shape, affine, crs)
[docs] @staticmethod def from_bbox( bbox: Union[BoundingBox, Tuple[float, float, float, float]], crs: MaybeCRS = None, *, tight: bool = False, shape: Union[SomeShape, int, None] = None, resolution: Optional[SomeResolution] = None, anchor: GeoboxAnchor = "default", tol: float = 0.01, ) -> "GeoBox": """ Construct :py:class:`~odc.geo.geobox.GeoBox` from a bounding box. :param bbox: Bounding box in CRS units, lonlat is assumed when ``crs`` is not supplied :param crs: CRS of the bounding box (defaults to EPSG:4326) :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. :param resolution: Use specified resolution :param tight: Supplying ``tight=True`` turns off pixel snapping. :param anchor: By default snaps grid such that pixel edges fall on X/Y axis. Ignored when tight mode is used. :param tol: Fraction of a pixel that can be ignored, defaults to 1/100. Bounding box of the output geobox is allowed to be smaller than supplied bounding box by that amount. :return: :py:class:`~odc.geo.geobox.GeoBox` that covers supplied bounding box. """ # pylint: disable=too-many-locals, too-many-branches anchor = _norm_anchor(anchor) _snap: Optional[XY[float]] = None if tight: anchor = AnchorEnum.FLOATING if isinstance(anchor, XY): _snap = anchor if anchor == AnchorEnum.EDGE: _snap = xy_(0, 0) elif anchor == AnchorEnum.CENTER: _snap = xy_(0.5, 0.5) def _norm_bbox( bbox: Tuple[float, float, float, float], crs: MaybeCRS ) -> BoundingBox: if isinstance(crs, str): if crs.lower().startswith("utm"): return BoundingBox(*bbox, crs="epsg:4326").to_crs(crs) return BoundingBox(*bbox, crs=(crs or "epsg:4326")) if not isinstance(bbox, BoundingBox): bbox = _norm_bbox(bbox, crs) elif bbox.crs is None: bbox = _norm_bbox(bbox.bbox, crs) if isinstance(shape, (int, float)): if bbox.aspect > 1: resolution = bbox.span_x / shape else: resolution = bbox.span_y / shape shape = None if resolution is not None: rx, ry = res_(resolution).xy if _snap is None: offx, nx = snap_grid(bbox.left, bbox.right, rx, None, tol=tol) offy, ny = snap_grid(bbox.bottom, bbox.top, ry, None, tol=tol) else: offx, nx = snap_grid(bbox.left, bbox.right, rx, _snap.x, tol=tol) offy, ny = snap_grid(bbox.bottom, bbox.top, ry, _snap.y, tol=tol) affine = Affine.translation(offx, offy) * Affine.scale(rx, ry) return GeoBox((ny, nx), crs=bbox.crs, affine=affine) if shape is None: raise ValueError("Must supply shape or resolution") shape = shape_(shape) nx, ny = shape.wh if nx == 0 or ny == 0: raise ValueError( f"Got zero size in ({nx}, {ny}): cannot calculate resolution." ) if nx < 0 or ny < 0: raise ValueError( f"Got negative size in ({nx}, {ny}): negative dimensions are not allowed." ) rx = bbox.span_x / nx ry = -bbox.span_y / ny if _snap is None: offx, offy = bbox.left, bbox.top else: offx, _ = snap_grid(bbox.left, bbox.right, rx, _snap.x, tol=tol) offy, _ = snap_grid(bbox.bottom, bbox.top, ry, _snap.y, tol=tol) affine = Affine.translation(offx, offy) * Affine.scale(rx, ry) return GeoBox((ny, nx), crs=bbox.crs, affine=affine)
[docs] @staticmethod def from_geopolygon( geopolygon: Geometry, resolution: Optional[SomeResolution] = None, crs: MaybeCRS = None, align: Optional[XY[float]] = None, *, shape: Union[SomeShape, int, None] = None, tight: bool = False, anchor: GeoboxAnchor = "default", tol: float = 0.01, ) -> "GeoBox": """ Construct :py:class:`~odc.geo.geobox.GeoBox` from a polygon. :param resolution: Either a single number or a :py:class:`~odc.geo.types.Resolution` object. :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. :param crs: CRS to use, if different from the geopolygon :param align: Deprecated: please switch to ``anchor=`` :param anchor: By default snaps grid such that pixel edges fall on X/Y axis. :param tol: Fraction of a pixel that can be ignored, defaults to 1/100. Bounding box of the output geobox is allowed to be smaller than supplied bounding box by that amount. :param tight: Supplying ``tight=True`` turns off pixel snapping. """ if align is not None: # support old-style "align", which is basically anchor but in CRS units ax, ay = align.xy if ax == 0 and ay == 0: anchor = AnchorEnum.EDGE else: assert resolution is not None resolution = res_(resolution) anchor = xy_(ax / abs(resolution.x), ay / abs(resolution.y)) if crs is None or isinstance(crs, Unset): crs = geopolygon.crs else: geopolygon = geopolygon.to_crs(crs) return GeoBox.from_bbox( geopolygon.boundingbox, crs, shape=shape, resolution=resolution, anchor=anchor, tol=tol, tight=tight, )
@staticmethod def from_rio(rdr) -> "GeoBox": """ Construct GeoBox from rasterio. :param rdr: Openned :py:class:`rasterio.DatasetReader` :returns: :py:class:`~odc.geo.geobox.GeoBox` """ return GeoBox(rdr.shape, rdr.transform, rdr.crs)
[docs] def buffered(self, xbuff: float, ybuff: Optional[float] = None) -> "GeoBox": """ Produce a tile buffered by ``xbuff, ybuff`` (in CRS units). """ if ybuff is None: ybuff = xbuff by, bx = ( _round_to_res(buf, res) for buf, res in zip((ybuff, xbuff), self.resolution.yx) ) affine = self._affine * Affine.translation(-bx, -by) ny, nx = (sz + 2 * b for sz, b in zip(self._shape, (by, bx))) return GeoBox( (ny, nx), affine=affine, crs=self._crs, )
[docs] def enclosing(self, region: Union[Geometry, BoundingBox]) -> "GeoBox": """ Construct compatible geobox covering given ``region``. Output GeoBox shares exactly the same pixel grid as source, but has different shape and location in the world. :param region: Region to be covered by the new GeoBox. """ if isinstance(region, BoundingBox): region = region.polygon if region.crs is None: raise ValueError("Must supply geo-resgistered region") pix_bbox = self.project(region).boundingbox.round() nx, ny = (max(1, int(span)) for span in (pix_bbox.span_x, pix_bbox.span_y)) tx, ty, *_ = pix_bbox.bbox A = self.translate_pix(tx, ty).affine return GeoBox(shape_((ny, nx)), A, self._crs)
def __getitem__(self, roi) -> "GeoBox": _shape, _affine = self.compute_crop(roi) return GeoBox(shape=_shape, 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 __hash__(self): return hash((*self._shape, self._crs, self._affine))
[docs] def overlap_roi(self, other: "GeoBox", tol: float = 1e-8) -> NormalizedROI: """ Compute overlap as ROI. Figure out slice into this geobox that shares pixels with the ``other`` geobox with consistent pixel grid. :raises: :py:class:`ValueError` when two geoboxes are not pixel-aligned. """ nx, ny = self._shape.xy x0, y0, x1, y1 = map(int, bounding_box_in_pixel_domain(other, self, tol)) x0, y0 = max(0, x0), max(0, y0) x1, y1 = min(x1, nx), min(y1, ny) return numpy.s_[y0:y1, x0:x1]
@property def transform(self) -> Affine: """Linear mapping from pixel space to CRS.""" return self._affine @property def affine(self) -> Affine: """ Linear mapping from pixel space to CRS. alias for :py:attr:`~odc.geo.geobox.GeoBox.transform` """ return self._affine @property def coordinates(self) -> Dict[str, Coordinate]: """ Query coordinates. This method only works with axis-aligned boxes. It will raise :py:class:`ValueError` if called on non-axis aligned :py:class:`~odc.geo.geobox.GeoBox`. :raises: :py:class:`ValueError` if not axis aligned. :return: Mapping from coordinate name to :py:class:`~odc.geo.geobox.Coordinate`. """ self._confirm_axis_aligned("Only axis aligned GeoBox can do this.") rx, _, tx, _, ry, ty, *_ = self._affine ny, nx = self._shape xs = numpy.arange(nx) * rx + (tx + rx / 2) ys = numpy.arange(ny) * ry + (ty + ry / 2) crs_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), crs_units, (ry, rx) ) ) coords = coordinates
[docs] def map_bounds(self) -> Tuple[Tuple[float, float], Tuple[float, float]]: """ Query bounds in folium/ipyleaflet style. Returns SW, and NE corners in lat/lon order. ``((lat_w, lon_s), (lat_e, lon_n))``. """ if self._crs is not None: (x0, y0), _, (x1, y1) = self.extent.exterior.to_crs("epsg:4326").points[:3] else: (x0, y0), _, (x1, y1) = self.extent.exterior.points[:3] return (y0, x0), (y1, x1)
[docs] def to_crs( self, 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": """ Compute GeoBox covering the same region in a different projection. :param crs: Desired CRS of the output :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" :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`` :return: Similar resolution, axis aligned geobox that fully encloses this one but in a different projection. """ # pylint: disable=import-outside-toplevel,too-many-arguments # can't be up-top due to circular imports issues from .overlap import compute_output_geobox return compute_output_geobox( self, crs, resolution=resolution, shape=shape, tight=tight, anchor=anchor, round_resolution=round_resolution, tol=tol, )
def __str__(self): return self.__repr__() def __repr__(self): return f"GeoBox({self._shape.yx!r}, {self._affine!r}, {self._crs!r})" def __eq__(self, other): if not isinstance(other, GeoBox): return False return ( self._shape == other._shape and self._affine == other._affine and self._crs == other._crs )
[docs] def __rmul__(self, transform: Affine) -> "GeoBox": """ Apply affine transform on CRS side. This has effect of transforming footprint of the source via ``transform``. :param transform: Affine matrix that shifts footprint of the source geobox. :return: :py:class:`~odc.geo.gebox.GeoBox` of the same pixel shape but covering different region. """ return GeoBox(self._shape, transform * self._affine, self._crs)
[docs] def __mul__(self, transform: Affine) -> "GeoBox": """ Apply affine transform on pixel side. ``X_old_pix = transform * X_new_pix`` :param transform: Affine matrix mapping from new pixel coordinate space to pixel coordinate space of input geobox. :returns: :py:class:`~odc.geo.gebox.GeoBox` of the same pixel shape but covering different region. Pixel coordinates in the output relate to input coordinates via ``transform``. """ return GeoBox(self._shape, self._affine * transform, self._crs)
[docs] def snap_to(self, other: "GeoBox") -> "GeoBox": """ Snap pixel grid to ``other``. Find smallest sub-pixel translation to apply to this geobox such that pixel edges align with ``other``. :param other: GeoBox to snap to, must be related by translation only, no change in scale or rotation. :raises: ``ValueError`` when ``other`` is in a different projection or has different resolution or orientation. """ _, subpix = split_translation(pixel_translation(other, self)) tx, ty = subpix.map(lambda x: maybe_zero(x, 1e-8)).xy return self.translate_pix(tx, ty)
[docs] def pad(self, padx: int, pady: MaybeInt = None) -> "GeoBox": """ Pad geobox. Expand GeoBox by fixed number of pixels on each side """ # false positive for -pady, it's never None by the time it runs # pylint: disable=invalid-unary-operand-type pady = padx if pady is None else pady ny, nx = self._shape.yx A = self._affine * Affine.translation(-padx, -pady) shape = (ny + pady * 2, nx + padx * 2) return GeoBox(shape, A, self._crs)
[docs] def pad_wh(self, alignx: int = 16, aligny: MaybeInt = None) -> "GeoBox": """ Possibly expand :py:class:`~odc.geo.geobox.GeoBox` by a few pixels. Find nearest ``width``/``height`` that are multiples of the desired factor. And return a new geobox that is slighly taller and/or wider covering roughly the same region. The new geobox will have the same CRS and transform but possibly larger shape. """ aligny = alignx if aligny is None else aligny ny, nx = (align_up(sz, n) for sz, n in zip(self._shape.yx, (aligny, alignx))) return GeoBox((ny, nx), self._affine, self._crs)
def crop(self, shape: SomeShape) -> "GeoBox": """ Crop or expand to a given shape. :returns: New :py:class:`~odc.geo.geobox.GeoBox` with a new shape, top left pixel remains at the same location and scale. """ return GeoBox(shape, self._affine, self._crs) expand = crop
[docs] def zoom_out(self, factor: float) -> "GeoBox": """ Compute :py:class:`~odc.geo.geobox.GeoBox` with changed resolution. - ``factor > 1`` implies smaller width/height, fewer but bigger pixels - ``factor < 1`` implies bigger width/height, more but smaller pixels :returns: GeoBox covering the same region but with different pixels (i.e. lower or higher resolution) """ _shape, _affine = self.compute_zoom_out(factor) return GeoBox(_shape, _affine, self._crs)
[docs] def zoom_to( self, shape: Union[SomeShape, int, float, None] = None, *, resolution: Optional[SomeResolution] = None, ) -> "GeoBox": """ Change GeoBox shape. When supplied a single integer scale longest dimension to match that. :returns: GeoBox covering the same region but with different number of pixels and therefore resolution. """ _shape, _affine = self.compute_zoom_to(shape, resolution=resolution) return GeoBox(_shape, _affine, self._crs)
[docs] def flipy(self) -> "GeoBox": """ Flip along Y axis. :returns: GeoBox covering the same region but with Y-axis flipped """ ny, _ = self._shape A = Affine.translation(0, ny) * Affine.scale(1, -1) return self * A
[docs] def flipx(self) -> "GeoBox": """ Flip along X axis. :returns: GeoBox covering the same region but with X-axis flipped """ _, nx = self._shape A = Affine.translation(nx, 0) * Affine.scale(-1, 1) return self * A
[docs] def translate_pix(self, tx: float, ty: float) -> "GeoBox": """ Shift GeoBox in pixel plane. ``(0,0)`` of the new GeoBox will be at the same location as pixel ``(tx, ty)`` in the original GeoBox. """ return self * Affine.translation(tx, ty)
@property def left(self) -> "GeoBox": """Same size geobox to the left of this one.""" return self.translate_pix(-self.shape.x, 0) @property def right(self) -> "GeoBox": """Same size geobox to the right of this one.""" return self.translate_pix(self.shape.x, 0) @property def top(self) -> "GeoBox": """Same size geobox directly above this one.""" return self.translate_pix(0, -self.shape.y) @property def bottom(self) -> "GeoBox": """Same size geobox directly below this one.""" return self.translate_pix(0, self.shape.y)
[docs] def rotate(self, deg: float) -> "GeoBox": """ Rotate GeoBox around the center. It's as if you stick a needle through the center of the GeoBox footprint and rotate it counter clock wise by supplied number of degrees. Note that from the pixel point of view image rotates the other way. If you have source image with an arrow pointing right, and you rotate GeoBox 90 degrees, in that view arrow should point down (this is assuming usual case of inverted y-axis) """ ny, nx = self._shape c0 = self._affine * (nx * 0.5, ny * 0.5) return Affine.rotation(deg, c0) * self
def _confirm_axis_aligned(self, raise_error: Optional[str] = None) -> bool: if is_affine_st(self._affine): return True if raise_error is not None: raise ValueError(raise_error) return False @property def axis_aligned(self): """ Check if Geobox is axis-aligned (not rotated). """ return self._confirm_axis_aligned() @property def center_pixel(self) -> "GeoBox": """ GeoBox of a center pixel. """ return self[self.shape.map(lambda x: x // 2).yx] @property def compat(self): """ Convert to :py:class:`datacube.utils.geometry.GeoBox` """ try: dc_geom = importlib.import_module("datacube.utils.geometry") except ModuleNotFoundError: return None w, h = self.shape.wh return dc_geom.GeoBox(w, h, self._affine, str(self._crs)) def __dask_tokenize__(self): return ( "odc.geo.geobox.GeoBox", str(self.crs), *self._shape.yx, *self._affine[:6], )
[docs] def gbox_boundary(gbox: GeoBoxBase, pts_per_side: int = 16) -> numpy.ndarray: """Alias for :py:meth:`odc.geo.geobox.GeoBox.boundary`.""" return gbox.boundary(pts_per_side)
def pixel_translation(a: GeoBox, b: GeoBox) -> XY[float]: """ Compute pixel translation ``a -> b``. Both geoboxes should have the same CRS and resolution. """ if a.crs != b.crs: raise ValueError("Geobox CRSs must match") # compute pixel-to-pixel transform # expect it to be a pure, pixel aligned translation # 1 0 tx # 0 1 ty # 0 0 1 # Such that tx,ty are almost integer. sx, z1, tx, z2, sy, ty, *_ = ~b.affine * a.affine if not ( numpy.isclose(sx, 1) and numpy.isclose(z1, 0) and numpy.isclose(z2, 0) and numpy.isclose(sy, 1) ): raise ValueError("Incompatible grids") return xy_(tx, ty)
[docs] def bounding_box_in_pixel_domain( geobox: GeoBox, reference: GeoBox, tol: float = 1e-8 ) -> BoundingBox: """ Bounding box of ``geobox`` in pixel space of ``reference``. :return: The bounding box of ``geobox`` with respect to the pixel grid defined by ``reference`` when their coordinate grids are compatible. Two geoboxes are compatible when they have the same CRS, same pixel size and orientation, and are related by whole pixel translation. :raises: :py:class:`ValueError` when two geoboxes are not pixel-aligned. """ # offset of ``geobox`` in ``reference`` pixels tx, ty = pixel_translation(geobox, reference).xy if not (is_almost_int(tx, tol) and is_almost_int(ty, tol)): raise ValueError("Incompatible grids") tx, ty = round(tx), round(ty) ny, nx = geobox.shape return BoundingBox(tx, ty, tx + nx, ty + ny, None)
[docs] def geobox_union_conservative(geoboxes: List[GeoBox]) -> GeoBox: """ Union of geoboxes as a geobox. 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(shape=bbox.shape, affine=affine, crs=reference.crs)
[docs] 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, crs=bbox.crs, ) if bbox.bottom > bbox.top: bbox = BoundingBox( left=bbox.left, bottom=bbox.bottom, right=bbox.right, top=bbox.bottom, crs=bbox.crs, ) affine = reference.affine * Affine.translation(*bbox[:2]) return GeoBox(shape=bbox.shape, affine=affine, crs=reference.crs)
[docs] def scaled_down_geobox(src_geobox: GeoBox, scaler: int) -> GeoBox: """ Compute :py:class:`~odc.geo.geobox.GeoBox` of a zoomed image. Given a source geobox and an 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:: We assume that pixel coordinates are ``0,0`` at the top-left corner of a top-left pixel. """ assert scaler > 1 ny, nx = (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((ny, nx), 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))
[docs] def flipy(gbox: GeoBox) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.flipy`.""" return gbox.flipy()
[docs] def flipx(gbox: GeoBox) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.flipx`.""" return gbox.flipx()
[docs] def translate_pix(gbox: GeoBox, tx: float, ty: float) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.GeoBox.translate_pix`.""" return gbox.translate_pix(tx, ty)
[docs] def pad(gbox: GeoBox, padx: int, pady: MaybeInt = None) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.GeoBox.pad`.""" return gbox.pad(padx, pady)
[docs] def pad_wh(gbox: GeoBox, alignx: int = 16, aligny: MaybeInt = None) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.GeoBox.pad_wh`.""" return gbox.pad_wh(alignx, aligny)
[docs] def zoom_out(gbox: GeoBox, factor: float) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.GeoBox.zoom_out`.""" return gbox.zoom_out(factor)
[docs] def zoom_to( gbox: GeoBox, shape: Union[SomeShape, int, float, None] = None, *, resolution: Optional[SomeResolution] = None, ) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.GeoBox.zoom_to`.""" return gbox.zoom_to(shape, resolution=resolution)
[docs] def rotate(gbox: GeoBox, deg: float) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.GeoBox.`.""" return gbox.rotate(deg)
[docs] def affine_transform_pix(gbox: GeoBox, transform: Affine) -> GeoBox: """Alias for :py:meth:`odc.geo.geobox.GeoBox.__mul__`.""" return gbox * transform
[docs] class GeoboxTiles: """Partition GeoBox into sub geoboxes.""" __slots__ = ("_gbox", "_tiles")
[docs] def __init__( self, box: GeoBoxBase, tile_shape: Union[SomeShape, Chunks2d, None], *, _tiles: Optional[RoiTiles] = None, ): """ Construct from a :py:class:`~odc.geo.GeoBox`. :param box: source :py:class:`~odc.geo.GeoBox` :param tile_shape: Shape of sub-tiles in pixels ``(rows, cols)`` """ self._gbox = box if _tiles is not None: self._tiles = _tiles else: assert tile_shape is not None self._tiles = roi_tiles(box.shape, tile_shape)
@property def base(self) -> GeoBoxBase: """Access base Geobox""" return self._gbox @property def shape(self) -> Shape2d: """Number of tiles along each dimension.""" return self._tiles.shape @property def roi(self) -> RoiTiles: """ Access ROI covered by tile. .. code-block:: python gbt = GeoboxTiles(..) roi = gbt.roi[0, 3] """ return self._tiles
[docs] def chunk_shape(self, idx: SomeIndex2d) -> Shape2d: """ Query chunk shape for a given chunk. :param idx: ``(row, col)`` chunk index :returns: ``(nrows, ncols)`` shape of a tile (edge tiles might be smaller) :raises: :py:class:`IndexError` when index is outside of ``[(0,0) -> .shape)``. """ return self._tiles.tile_shape(idx)
@property def chunks(self) -> Chunks2d: return self._tiles.chunks def _crop(self, roi: ROI) -> "GeoboxTiles": gbox_new = self.base[self._tiles[roi]] return GeoboxTiles(gbox_new, (0, 0), _tiles=self._tiles.crop(roi)) def clip( self, selection: Sequence[Tuple[int, int]] ) -> Tuple["GeoboxTiles", List[Tuple[int, int]]]: """ Crop to a set of tiles. Returns cropped version of :py:class:`GeoboxTiles` and a list of tile coordinates in the new, cropped space. """ tiles, roi, new_idx = clip_tiles(self._tiles, selection) return GeoboxTiles(self[roi], None, _tiles=tiles), new_idx @property def crop(self) -> Mapping[ROI, "GeoboxTiles"]: return func2map(self._crop) def __getitem__(self, idx: Union[SomeIndex2d, ROI]) -> GeoBoxBase: """ Lookup tile by index, index is in matrix access order: ``(row, col)``. :param idx: ``(row, col)`` index :returns: GeoBox of a tile :raises: IndexError when index is outside of ``[(0,0) -> .shape)`` """ return self._gbox[self._tiles[idx]] def pix_bbox(self, idx: Union[SomeIndex2d, ROI]) -> BoundingBox: """ BoundingBox in pixel space of a given tile. """ ry, rx = self._tiles[idx] return BoundingBox(rx.start, ry.start, rx.stop, ry.stop)
[docs] def range_from_bbox(self, bbox: BoundingBox) -> Tuple[range, range]: """ Intersect with a bounding box. Compute rows and columns overlapping with a given :py:class:`~odc.geo.geom.BoundingBox`. """ if bbox.crs is not None: bbox = self._gbox.project(bbox.polygon).boundingbox def _clamp(span: Tuple[float, float], N: int): a1, a2 = span a1 = int(clamp(math.floor(a1), 0, N - 1)) a2 = int(clamp(math.ceil(a2), 1, N)) - 1 return a1, a2 NY, NX = self._gbox.shape.yx x1, x2 = _clamp(bbox.range_x, NX) y1, y2 = _clamp(bbox.range_y, NY) y1, x1 = self._tiles.locate((y1, x1)) y2, x2 = self._tiles.locate((y2, x2)) return range(y1, y2 + 1), range(x1, x2 + 1)
def _tiles_from_pix_bbox(self, bbox: BoundingBox) -> Iterator[Tuple[int, int]]: yy, xx = self.range_from_bbox(bbox) yield from itertools.product(yy, xx)
[docs] def tiles(self, query: Union[Geometry, BoundingBox]) -> Iterator[Tuple[int, int]]: """Return tile indexes overlapping with a given geometry.""" target_crs = self._gbox.crs if isinstance(query, BoundingBox): if query.crs is None: # special case for bounding box in pixel domain yield from self._tiles_from_pix_bbox(query) return poly = query.polygon else: poly = query if target_crs is not None and poly.crs != target_crs: poly = poly.to_crs(target_crs, check_and_fix=True) yy, xx = self.range_from_bbox(poly.boundingbox) for idx in itertools.product(yy, xx): gbox = self[idx] if not poly.disjoint(gbox.extent): yield idx
def _check_linear(self, src: "GeoboxTiles") -> Optional[Affine]: if src.base.crs != self.base.crs: return None if not isinstance(self.base, GeoBox): return None if not isinstance(src.base, GeoBox): return None # src_pix = A*dst_pix A = snap_affine((~src.base.transform) * self.base.transform) if is_affine_st(A): return A return None def _all_tiles(self) -> Iterator[Tuple[int, int]]: for iy, ix in numpy.ndindex(self.shape.shape): yield (iy, ix) def _grid_intersect_linear( self, src: "GeoboxTiles", A: Affine, ) -> Dict[Tuple[int, int], List[Tuple[int, int]]]: deps: Dict[Tuple[int, int], List[Tuple[int, int]]] = {} for idx in self._all_tiles(): bbox = self.pix_bbox(idx).transform(A).round() src_idx = list(src.tiles(bbox)) deps[idx] = src_idx return deps
[docs] def grid_intersect( self, src: "GeoboxTiles" ) -> Dict[Tuple[int, int], List[Tuple[int, int]]]: """ Figure out tile to tile overlap graph between two grids. For every tile in this :py:class:`GeoboxTiles` find every tile in ``other`` that intersects with this ``tile``. """ A = self._check_linear(src) if A is not None: return self._grid_intersect_linear(src, A) if src.base.crs == self.base.crs: src_footprint = src.base.extent else: # compute "robust" source footprint in CRS of self via espg:4326 src_footprint = ( src.base.footprint(4326, 2) & self.base.footprint(4326, 2) ).to_crs(self.base.crs) xy_chunks_with_data = list(self.tiles(src_footprint)) deps: Dict[Tuple[int, int], List[Tuple[int, int]]] = {} for idx in xy_chunks_with_data: geobox = self[idx] deps[idx] = list(src.tiles(geobox.extent)) return deps
def __dask_tokenize__(self): return ( "odc.geo.geobox.GeoboxTiles", *self._gbox.__dask_tokenize__()[1:], *self._tiles.__dask_tokenize__()[1:], ) def __str__(self): return str(self.roi) def __eq__(self, __value: object) -> bool: if not isinstance(__value, GeoboxTiles): return False if self is __value: return True return self._tiles == __value._tiles and self._gbox == __value._gbox __repr__ = __str__