Source code for odc.geo.gridspec

# 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
"""GridSpec class."""
import math
from typing import Any, Dict, Iterator, Optional, Tuple

from affine import Affine

from . import geom
from .crs import SomeCRS, norm_crs_or_error
from .geobox import GeoBox
from .geom import BoundingBox, Geometry, intersects
from .math import Bin1D
from .types import (
    XY,
    Index2d,
    Shape2d,
    SomeIndex2d,
    SomeResolution,
    SomeShape,
    ixy_,
    res_,
    resyx_,
    shape_,
    xy_,
    yx_,
)


[docs]class GridSpec: """ Definition for a regular spatial grid. :param crs: Coordinate System used to define the grid :param tile_shape: Size of each tile in pixels :param resolution: Size of each data point in the grid, in CRS units. ``Y`` will usually be negative. :param origin: Coordinate of a bottom-left corner of the ``(0,0)`` tile in CRS units. Default is ``xy_(0.0, 0.0)`` :param flipx: when ``True`` grid index for X axis increments left to right :param flipy: when ``True`` grid index for Y axis increments top to bottom """
[docs] def __init__( self, crs: SomeCRS, tile_shape: SomeShape, resolution: SomeResolution, origin: Optional[XY[float]] = None, flipx: bool = False, flipy: bool = False, ): tile_shape = shape_(tile_shape) resolution = res_(resolution) if origin is None: origin = xy_(0.0, 0.0) else: assert isinstance(origin, XY) self.crs = norm_crs_or_error(crs) self._shape = tile_shape self.resolution = resolution self.tile_size = xy_( tile_shape.x * abs(resolution.x), tile_shape.y * abs(resolution.y), ) self.origin = origin ox, oy = origin.xy self._ybin = Bin1D(self.tile_size.y, oy, -1 if flipy else 1) self._xbin = Bin1D(self.tile_size.x, ox, -1 if flipx else 1)
def __eq__(self, other): if not isinstance(other, GridSpec): return False return ( self._shape == other._shape and self._ybin == other._ybin and self._xbin == other._xbin and self.crs == other.crs ) @property def dimensions(self) -> Tuple[str, str]: """List of dimension names of the grid spec.""" return self.crs.dimensions @property def alignment(self) -> XY[float]: """Pixel boundary alignment.""" y, x = ( orig % abs(res) for orig, res in zip(self.origin.yx, self.resolution.yx) ) return yx_(y, x) @property def tile_shape(self) -> Shape2d: """Tile shape in pixels (Y,X order, like numpy).""" return self._shape
[docs] def pt2idx(self, x: float, y: float) -> Index2d: """ Compute tile index from a point. :param x: X coordinate of a point in CRS units :param y: Y coordinate of a point in CRS units :return: ``(ix, iy)``, index of a tile containing given point """ return ixy_(self._xbin.bin(x), self._ybin.bin(y))
def _tile_txy(self, tile_index: Index2d) -> XY[float]: """Location of 0,0 pixel in CRS units.""" ix, iy = tile_index.xy rx, ry = self.resolution.xy x0, x1 = self._xbin[ix] tx = x0 if rx > 0 else x1 y0, y1 = self._ybin[iy] ty = y0 if ry > 0 else y1 return xy_(tx, ty)
[docs] def tile_geobox(self, tile_index: SomeIndex2d) -> GeoBox: """ Tile geobox. :param tile_index: ``(ix, iy)`` """ tile_index = ixy_(tile_index) tx, ty = self._tile_txy(tile_index).xy rx, ry = self.resolution.xy return GeoBox(self._shape, crs=self.crs, affine=Affine(rx, 0, tx, 0, ry, ty))
[docs] def __getitem__(self, idx: SomeIndex2d) -> GeoBox: """Lookup :py:class:`~odc.geo.geobox.GeoBox` of a given tile.""" return self.tile_geobox(idx)
def idx_bounds(self, bounds: BoundingBox) -> BoundingBox: """ Convert bounds from CRS to index space. :param bounds: Query bounding box :return: Bounding box in tile index space """ x1, y1, x2, y2 = bounds ix1, iy1 = self.pt2idx(x1, y1).xy ix2, iy2 = self.pt2idx(x2, y2).xy ix1, ix2 = sorted([ix1, ix2]) iy1, iy2 = sorted([iy1, iy2]) return BoundingBox(ix1, iy1, ix2, iy2)
[docs] def tiles( self, bounds: BoundingBox, geobox_cache: Optional[dict] = None ) -> Iterator[Tuple[Tuple[int, int], GeoBox]]: """ Query tiles overlapping with bounding box. Output is a sequence of ``tile_index``, :py:class:`odc.geo.geobox.GeoBox` tuples. .. note:: Grid cells are referenced by coordinates ``(x, y)``, which is the opposite to the usual CRS dimension order. :param bounds: Boundary coordinates of the required grid :param geobox_cache: Optional cache to re-use geoboxes instead of creating new one each time :return: Iterator of tuples of grid index and corresponding :py:class:`odc.geo.geobox.GeoBox` """ def geobox(tile_index): if geobox_cache is None: return self.tile_geobox(tile_index) gbox = geobox_cache.get(tile_index) if gbox is None: gbox = self.tile_geobox(tile_index) geobox_cache[tile_index] = gbox return gbox ix1, iy1, ix2, iy2 = self.idx_bounds(bounds) for iy in range(iy1, iy2 + 1): for ix in range(ix1, ix2 + 1): tile_index = (ix, iy) yield tile_index, geobox(tile_index)
[docs] def tiles_from_geopolygon( self, geopolygon: Geometry, geobox_cache: Optional[dict] = None, ) -> Iterator[Tuple[Tuple[int, int], GeoBox]]: """ Query tiles overlapping with a given polygon. Output is a sequence of ``tile_index``, :py:class:`odc.geo.geobox.GeoBox` tuples. :param geopolygon: Polygon to tile :param geobox_cache: Optional cache to re-use geoboxes instead of creating new one each turn: iterator of grid cells with :py:class:`odc.geo.geobox.GeoBox` tiles """ geopolygon = geopolygon.to_crs(self.crs) bbox = geopolygon.boundingbox for tile_index, tile_geobox in self.tiles(bbox, geobox_cache): if intersects(tile_geobox.extent, geopolygon): yield (tile_index, tile_geobox)
def __str__(self) -> str: return f"GridSpec(crs={self.crs}, tile_shape={self._shape}, resolution={self.resolution})" def __repr__(self) -> str: return self.__str__()
[docs] def geojson( self, *, bbox: Optional[BoundingBox] = None, geopolygon: Optional[Geometry] = None, ) -> Dict[str, Any]: """ Render to GeoJSON. :param bbox: Limit output to tiles overlapping with the given bounding box (native CRS of the grid). :param geopolygon: Limit output to tiles overlapping with the given geometry (any CRS) :return: GeoJSON representation of the grid spec """ if geopolygon is not None: _tiles = self.tiles_from_geopolygon(geopolygon) elif bbox is not None: _tiles = self.tiles(bbox) else: valid_region = self.crs.valid_region if valid_region is None: valid_region = geom.box(-180, -90, 180, 90, "epsg:4326") _tiles = self.tiles( valid_region.buffer(-0.05).to_crs(self.crs, resolution=0.5).boundingbox ) props = { "native_crs": str(self.crs), "tile_shape": self._shape, "resolution": self.resolution, } features = [] for (ix, iy), geobox in _tiles: features.append(geobox.extent.geojson(idx=f"{ix},{iy}")) return {"type": "FeatureCollection", "features": features, "properties": props}
[docs] @staticmethod def from_sample_tile( box: Geometry, *, shape: SomeShape = (-1, -1), idx: SomeIndex2d = (0, 0), flipx: bool = False, flipy: bool = False, ) -> "GridSpec": """ Construct :py:class:`odc.geo.gridspec.GridSpec` from a sample tile. Bounding box of one tile, it's index, and a shape of the tile in pixels fully define a grid. This method could be more convenient than canonical representation. :param box: Geometry of the tile in some CRS :param idx: ``ix, iy`` index of the sample tile, default ``(0, 0)`` :param shape: ``height, width`` of the tile, must be supplied :param flipx: when ``True`` grid index for X axis increments left to right :param flipy: when ``True`` grid index for Y axis increments top to bottom """ if shape == (-1, -1): raise ValueError("Must specify shape of the tile in pixels") shape = shape_(shape) idx = ixy_(idx) crs = norm_crs_or_error(box.crs) ix, iy = idx.xy nx, ny = shape.xy bbox = box.boundingbox xbin = Bin1D.from_sample_bin(ix, bbox.range_x, -1 if flipx else 1) ybin = Bin1D.from_sample_bin(iy, bbox.range_y, -1 if flipy else 1) origin = xy_(xbin.origin, ybin.origin) resolution = resyx_(-ybin.sz / ny, xbin.sz / nx) return GridSpec( crs, shape, resolution=resolution, origin=origin, flipx=flipx, flipy=flipy )
[docs] @staticmethod def web_tiles(zoom: int, npix: int = 256) -> "GridSpec": """ Construct :py:class:`~odc.geo.gridspec.GridSpec` that matches slippy tiles. Tile with index ``(0, 0)`` is at the top left corner of the map, and tile with index ``(2^zoom - 1, 2^zoom - 1)`` is at the bottom right. :param zoom: Zoom level, ``0`` is one single tile, ``1`` is ``2x2``, ``3`` is ``8x8``... :param npix: Usually tiles are ``256x256`` pixels wide. But you can override that. :return: Grid spec that encodes slippy tiles scheme in ``EPSG:3857``. """ R = 6_378_137 pi = math.pi tsz = pi * R * (2 ** (1 - zoom)) # in meters x, y = -pi * R, pi * R # top-left corner of tile 0,0 tile0 = geom.box(x, y - tsz, x + tsz, y, "epsg:3857") shape = (npix, npix) return GridSpec.from_sample_tile(tile0, shape=shape, idx=(0, 0), flipy=True)