Coverage for stems/gis/grids.py : 88%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
""" Create and work with geospatial data tiling schemes
A tile scheme defines the grid which many products are based on. Tiling schemes define a projection, usually a projection applicable over a wide area, the size (in pixels) of each tile, and georeferencing information that defines the upper left coordinate and pixel sizes (thus defining the posting of each pixel). Tiles are defined by the tile grid coordinate (horizontal & vertical indices) that determines the number of tiles offset from the upper left coordinate of the tiling scheme, such that one can retrieve the real-world coordinate of a pixel if the tile grid index and pixel row/column within the tile is known. """
""" A tile grid specification for gridding data
Attributes ---------- ul : tuple Upper left X/Y coordinates crs : rasterio.crs.CRS Coordinate system information res : tuple Pixel X/Y resolution size : tuple Number of pixels in X/Y dimensions for each tile limits : tuple[tuple, tuple] Maximum and minimum rows (vertical) and columns (horizontal) given as ((row_start, row_stop), (col_start, col_stop)). Used to limit access to Tiles beyond domain. name : str, optional Name of this tiling scheme """
""" Return this TileGrid as a dictionary (e.g., to serialize)
Returns ------- dict TileGrid attributes needed to re-initialize class """ 'ul': tuple(self.ul), 'crs': self.crs_wkt, 'res': tuple(self.res), 'size': tuple(self.size), 'limits': tuple(self.limits), 'name': self.name }
def from_dict(cls, d): """ Return a ``TileGrid`` from a dictionary
Parameters ---------- d : dict Dictionary of class attributes (see __init__)
Returns ------- TileGrid A new TileGrid according to parameters in ``d`` """
f'<{self.__class__.__name__} at {hex(id(self))}>', f' * name: {self.name}', f' * ul={self.ul}', f' * crs={self.crs_wkt}', f' * res={self.res}', f' * size={self.size}', f' * limits={self.limits}' ])
warnings.warn( f"{self.__class__.__name__} '{self.name}' does not specify a " f"``limits``, so has an unlimited number of rows/columns. " f"Defaulting to {_DEFAULT_TILEGRID_UNSIZED_LIMITS}" ) return ((0, _DEFAULT_TILEGRID_UNSIZED_LIMITS, ), ) * 2 else:
def rows(self): """ list[int]: the vertical/row numbers for all tiles """
def cols(self): """ list[int]: the horizontal/column numbers for all tiles """
def nrow(self): """ int: the number of rows in the TileGrid """
def ncol(self): """ int: the number of columns in the TileGrid """
def transform(self): """ affine.Affine : the affine transform for this TileGrid """ 0, -self.res[1], self.ul[1])
rows=None, cols=None, rfc7946=False, skip_invalid=False): """ Returns this grid of tiles as GeoJSON
Parameters ---------- crs : rasterio.crs.CRS Coordinate reference system of output. Defaults to EPSG:4326 per GeoJSON standard (RFC7946). If ``None``, will return geometries in TileGrid's CRS rows : Sequence[int], optional If this TileGrid was not given ``limits`` or if you want a subset of the tiles, specify the rows to map cols : Sequence[int], optional If this TileGrid was not given ``limits`` or if you want a subset of the tiles, specify the rows to map rfc7946 : bool, optional Return GeoJSON compliant with RFC7946. Helps fix GeoJSON that crosses the anti-meridian/datelines by splitting polygons into multiple as needed skip_invalid : bool, optional If ``True``, checks for tile bounds for invalid geometries and will only include valid tile geometries.
Returns ------- dict GeoJSON """ rows_ = rows or self.rows cols_ = cols or self.cols tile_rc = itertools.product(rows_, cols_)
features = [] for r, c in tile_rc: tile = self[r, c] feat = tile.geojson(crs=crs) if skip_invalid and geom.is_null(feat['geometry']): continue features.append(feat)
geojson = { "type": "FeatureCollection", "features": features } if rfc7946: return geom.fix_geojson_rfc7946(geojson) else: return geojson
# `Mapping` ABC requirement """ Return a Tile for the grid row/column specified by index """
# TILE ACCESS METHODS """ Return a Tile containing a given point (x, y)
Parameters ---------- point : tuple X/Y coordinates. Coordinates must be in the same CRS as the TileGrid.
Returns ------- tile : stems.gis.grids.Tile The intersecting :class`Tile` """
""" Yield Tile objects for this grid within a given bounds
Parameters ---------- bounds : BoundingBox Input bounds. Bounds must be in the same CRS as the TileGrid.
Yields ------ tile : Tile Tiles that are within within provided bounds """
""" Yield tiles within a Region of Interest (ROI) ``shapely`` geometry
Parameters ---------- roi : Polygon, MultiPolygon, etc. A shapely geometry. Must be in the same CRS as the TileGrid.
Yields ------ iterable[Tile] Yields ``Tile`` objects within provided Region of Interest """
# HELPERS """ Return Tile footprint bounds for given index
Parameters ---------- index : tuple Tile row/column index
Returns ------- bbox : BoundingBox Bounding box of tile """ left=self.ul[0] + index[1] * self.size[0] * self.res[0], right=self.ul[0] + (index[1] + 1) * self.size[0] * self.res[0], top=self.ul[1] - index[0] * self.size[1] * self.res[1], bottom=self.ul[1] - (index[0] + 1) * self.size[1] * self.res[1] )
""" Return a valid tile index, or raise an error """ # TODO: this could be simpler if we store (h, v) in a np.array # because we could borrow fancy indexing to support returning # a scalar (row/col), or a Sequence (e.g., grid[:, 0], # grid[[0, 1], [0]], etc)
# Only support row/col for now raise IndexError('TileSpec only has two dimensions (row/col)') # Need to be ints, and we'll cast as built-in int else:
# Check if index is outside of limits col < min(col_lim) or col > max(col_lim)): f'"{self.name}" limits ({self.limits}).')
""" Return the Tile for given index
Parameters ---------- index : tuple Tile row/column index
Returns ------- tile : Tile Tile at index
Raises ------ IndexError Raise if requested Tile is out of domain (outside ``limits``) """ self.res, self.size)
# yield `Tile`s that intersect # TODO - this should really a check that we can turn off tile.bbox.touches(bbox_poly)):
# return y/x tile index that intersect `bounds` range(min_grid_x, max_grid_x + 1))
""" A Tile
Attributes ---------- index : tuple[int, int] The (row, column) index of this tile in the grid crs : rasterio.crs.CRS The tile coordinate reference system bounds : BoundingBox The bounding box of the tile res : tuple[float, float] Pixel X/Y resolution size : tuple[int, int] Number of columns and rows in tile """
return False
return hash((self.index, self.crs.wkt, self.bounds, self.res, self.size))
""" Return a ``dict`` representing this Tile
Returns ------- dict This Tile as a dict (CRS will be a WKT for portability) """ 'index': self.index, 'crs': self.crs.wkt, 'bounds': self.bounds, 'res': self.res, 'size': self.size }
def from_dict(cls, d): """ Create a Tile from a dict
Parameters ---------- d : dict Tile info, including "index", "crs" (a WKT), "bounds", "res", and "size" """ convert.to_crs(d['crs']), d['bounds'], tuple(d['res']), tuple(d['size']))
def vertical(self): """ int: The horizontal index of this tile in its tile specification """
def horizontal(self): """ int: The horizontal index of this tile in its tile specification """
def transform(self): """ affine.Affine: The ``Affine`` transform for the tile """ return Affine(self.res[0], 0, self.bounds.left, 0, -self.res[1], self.bounds.top)
def bbox(self): """ shapely.geometry.Polygon: This tile's bounding box geometry """
def width(self): """ int : The number of columns in this Tile """
def height(self): """ int : The number of columns in this Tile """
""" Return y/x pixel coordinates
Parameters ---------- center : bool, optional Return coordinates for pixel centers (default)
Returns ------- np.ndarray Y coordinates np.ndarray X coordinates """ return transform_to_coords(self.transform, width=self.width, height=self.height, center=center)
""" Return this Tile's geometry as GeoJSON
Parameters ---------- crs : rasterio.crs.CRS Coordinate reference system of output. Defaults to EPSG:4326 per GeoJSON standard (RFC7946). If ``None``, will return geometries in Tile's CRS
Returns ------- dict This tile's geometry and crs represented as GeoJSON
References ---------- .. [1] https://tools.ietf.org/html/rfc7946#page-12 """
else: logger.debug('Not reprojecting GeoJSON since output CRS ' 'is the same as the Tile CRS')
'type': 'Feature', 'properties': { 'horizontal': self.horizontal, 'vertical': self.vertical }, 'geometry': gj }
""" Retrieve tile grid specifications from a YAML file
Parameters ---------- filename : str Filename of YAML data containing specifications. If ``None``, will load grids packaged with module
Returns ------- dict Mapping of (name, TileGrid) pairs for tile grid specifications. By default, returns tile specifications included with this package """
# ints of many types, if they work... else: |