-import os
-import rasterio
-from rasterio.warp import transform_bounds
-from rasterio.io import DatasetReader
-import math
-from rio_tiler.errors import TileOutsideBounds
-from . import utils
-import numpy as np
-
-
-[docs]def tile_utm_source(src, ll_x, ll_y, ur_x, ur_y, indexes=None, tilesize=256,
-
nodata=None, alpha=None, dst_crs='epsg:4326'):
-
"""
-
Create a UTM tile from a :py:class:`rasterio.Dataset` in memory.
-
-
Arguments
-
---------
-
src : :py:class:`rasterio.Dataset`
-
Source imagery dataset to tile.
-
ll_x : int or float
-
Lower left x position (i.e. Western bound).
-
ll_y : int or f
-
loat
-
Lower left y position (i.e. Southern bound).
-
ur_x : int or float
-
Upper right x position (i.e. Eastern bound).
-
ur_y : int or float
-
Upper right y position (i.e. Northern bound).
-
indexes : tuple of 3 ints, optional
-
Band indexes for the output. By default, extracts all of the
-
indexes from `src`.
-
tilesize : int, optional
-
Output image X and Y pixel extent. Defaults to ``256``.
-
nodata : int or float, optional
-
Value to use for `nodata` pixels during tiling. By default, uses
-
the existing `nodata` value in `src`.
-
alpha : int, optional
-
Alpha band index for tiling. By default, uses the same band as
-
specified by `src`.
-
dst_crs : str, optional
-
Coordinate reference system for output. Defaults to ``"epsg:4326"``.
-
-
Returns
-
-------
-
``(data, mask, window, window_transform)`` tuple.
-
data : :py:class:`numpy.ndarray`
-
int pixel values. Shape is ``(C, Y, X)`` if retrieving multiple
-
channels, ``(Y, X)`` otherwise.
-
mask : :py:class:`numpy.ndarray`
-
int mask indicating which pixels contain information and which are
-
`nodata`. Pixels containing data have value ``255``, `nodata`
-
pixels have value ``0``.
-
window : :py:class:`rasterio.windows.Window`
-
:py:class:`rasterio.windows.Window` object indicating the raster
-
location of the dataset subregion being returned in `data`.
-
window_transform : :py:class:`affine.Affine`
-
Affine transformation for the window.
-
-
"""
-
-
wgs_bounds = transform_bounds(
-
*[src.crs, dst_crs] + list(src.bounds), densify_pts=21)
-
-
indexes = indexes if indexes is not None else src.indexes
-
tile_bounds = (ll_x, ll_y, ur_x, ur_y)
-
if not utils.tile_exists_utm(wgs_bounds, tile_bounds):
-
raise TileOutsideBounds(
-
'Tile {}/{}/{}/{} is outside image bounds'.format(tile_bounds))
-
-
return utils.tile_read_utm(src, tile_bounds, tilesize, indexes=indexes,
-
nodata=nodata, alpha=alpha, dst_crs=dst_crs)
-
-
-[docs]def tile_utm(source, ll_x, ll_y, ur_x, ur_y, indexes=None, tilesize=256,
-
nodata=None, alpha=None, dst_crs='epsg:4326'):
-
"""
-
Create a UTM tile from a file or a :py:class:`rasterio.Dataset` in memory.
-
-
This function is a wrapper around :func:`tile_utm_source` to enable passing
-
of file paths instead of pre-loaded :py:class:`rasterio.Dataset` s.
-
-
Arguments
-
---------
-
source : :py:class:`rasterio.Dataset`
-
Source imagery dataset to tile.
-
ll_x : int or float
-
Lower left x position (i.e. Western bound).
-
ll_y : int or float
-
Lower left y position (i.e. Southern bound).
-
ur_x : int or float
-
Upper right x position (i.e. Eastern bound).
-
ur_y : int or float
-
Upper right y position (i.e. Northern bound).
-
indexes : tuple of 3 ints, optional
-
Band indexes for the output. By default, extracts all of the
-
indexes from `source` .
-
tilesize : :obj:`int`, optional
-
Output image X and Y pixel extent. Defaults to ``256``.
-
nodata : int or float, optional
-
Value to use for ``nodata`` pixels during tiling. By default, uses
-
the existing ``nodata`` value in `src`.
-
alpha : :obj:`int`, optional
-
Alpha band index for tiling. By default, uses the same band as
-
specified by `src`.
-
dst_crs : str, optional
-
Coordinate reference system for output. Defaults to ``"epsg:4326"``.
-
-
Returns
-
-------
-
``(data, mask, window, window_transform`` tuple.
-
data : :py:class:`numpy.ndarray`
-
int pixel values. Shape is ``(C, Y, X)`` if retrieving multiple
-
channels, ``(Y, X)`` otherwise.
-
mask : :class:`numpy.ndarray`
-
int mask indicating which pixels contain information and which are
-
`nodata`. Pixels containing data have value ``255``, `nodata` pixels
-
have value ``0``.
-
window : :py:class:`rasterio.windows.Window`
-
:py:class:`rasterio.windows.Window` object indicating the raster
-
location of the dataset subregion being returned in `data`.
-
window_transform : :py:class:`affine.Affine`
-
Affine transformation for the window.
-
-
"""
-
-
if isinstance(source, DatasetReader):
-
src = source
-
elif os.path.exists(source):
-
src = rasterio.open(source) # read in the file
-
else:
-
raise ValueError('Source is not a rasterio.Dataset or a valid path.')
-
-
return tile_utm_source(src, ll_x, ll_y, ur_x, ur_y, indexes=indexes,
-
tilesize=tilesize, nodata=nodata, alpha=alpha,
-
dst_crs=dst_crs)
-
-
-[docs]def get_chip(source, ll_x, ll_y, gsd,
-
utm_crs='',
-
indexes=None,
-
tilesize=256,
-
nodata=None,
-
alpha=None):
-
"""Get an image tile of specific pixel size.
-
-
This wrapper function permits passing of `ll_x`, `ll_y`, `gsd`, and
-
`tile_size_pixels` in place of boundary coordinates to extract an image
-
region of defined pixel extent.
-
-
Arguments
-
---------
-
source : :py:class:`rasterio.Dataset`
-
Source imagery dataset to tile.
-
ll_x : int or float
-
Lower left x position (i.e. Western bound).
-
ll_y : int or float
-
Lower left y position (i.e. Southern bound).
-
gsd : float
-
Ground sample distance of the source imagery in meter/pixel units.
-
utm_crs : :py:class:`rasterio.crs.CRS`, optional
-
UTM coordinate reference system string for the imagery. If not
-
provided, this is calculated using
-
:func:`cw_tiler.utils.get_wgs84_bounds` and
-
:func:`cw_tiler.utils.calculate_UTM_crs` .
-
indexes : tuple of 3 ints, optional
-
Band indexes for the output. By default, extracts all of the
-
indexes from `source`.
-
tilesize : int, optional
-
Output image X and Y pixel extent. Defaults to ``256`` .
-
nodata : int or float, optional
-
Value to use for `nodata` pixels during tiling. By default, uses
-
the existing `nodata` value in `source`.
-
alpha : int, optional
-
Alpha band index for tiling. By default, uses the same band as
-
specified by `source`.
-
-
Returns
-
-------
-
``(data, mask, window, window_transform`` tuple.
-
data : :class:`numpy.ndarray`
-
int pixel values. Shape is ``(C, Y, X)`` if retrieving multiple
-
channels, ``(Y, X)`` otherwise.
-
mask : :class:`numpy.ndarray`
-
int mask indicating which pixels contain information and which are
-
`nodata`. Pixels containing data have value ``255``, `nodata` pixels
-
have value ``0``.
-
window : :py:class:`rasterio.windows.Window`
-
:py:class:`rasterio.windows.Window` object indicating the raster
-
location of the dataset subregion being returned in `data`.
-
window_transform : :py:class:`affine.Affine`
-
Affine transformation for the window.
-
"""
-
-
ur_x = ll_x + gsd * tilesize
-
ur_y = ll_y + gsd * tilesize
-
-
if isinstance(source, DatasetReader):
-
src = source
-
else:
-
src = rasterio.open(source)
-
-
if not utm_crs:
-
wgs_bounds = utils.get_wgs84_bounds(src)
-
utm_crs = utils.calculate_UTM_crs(wgs_bounds)
-
-
return tile_utm(source, ll_x, ll_y, ur_x, ur_y, indexes=indexes,
-
tilesize=tilesize, nodata=nodata, alpha=alpha,
-
dst_crs=utm_crs)
-
-
-[docs]def calculate_anchor_points(utm_bounds, stride_size_meters=400, extend=False,
-
quad_space=False):
-
"""Get anchor point (lower left corner of bbox) for chips from a tile.
-
-
Arguments
-
---------
-
utm_bounds : tuple of 4 floats
-
A :obj:`tuple` of shape ``(min_x, min_y, max_x, max_y)`` that defines
-
the spatial extent of the tile to be split. Coordinates should be in
-
UTM.
-
stride_size_meters : int, optional
-
Stride size in both X and Y directions for generating chips. Defaults
-
to ``400``.
-
extend : bool, optional
-
Defines whether UTM boundaries should be rounded to the nearest integer
-
outward from `utm_bounds` (`extend` == ``True``) or
-
inward from `utm_bounds` (`extend` == ``False``). Defaults
-
to ``False`` (inward).
-
quad_space : bool, optional
-
If tiles will overlap by no more than half their X and/or Y extent in
-
each direction, `quad_space` can be used to split chip
-
anchors into four non-overlapping subsets. For example, if anchor
-
points are 400m apart and each chip will be 800m by 800m, `quad_space`
-
will generate four sets which do not internally overlap;
-
however, this would fail if tiles are 900m by 900m. Defaults to
-
``False``, in which case the returned ``anchor_point_list_dict`` will
-
comprise a single list of anchor points.
-
-
Returns
-
-------
-
anchor_point_list_dict : dict of list(s) of lists
-
-
If ``quad_space==True`` , `anchor_point_list_dict` is a
-
:obj:`dict` with four keys ``[0, 1, 2, 3]`` corresponding to the four
-
subsets of chips generated (see `quad_space` ). If
-
``quad_space==False`` , `anchor_point_list_dict` is a
-
:obj:`dict` with a single key, ``0`` , that corresponds to a list of all
-
of the generated anchor points. Each anchor point in the list(s) is an
-
``[x, y]`` pair of UTM coordinates denoting the SW corner of a chip.
-
-
"""
-
if extend:
-
min_x = math.floor(utm_bounds[0])
-
min_y = math.floor(utm_bounds[1])
-
max_x = math.ceil(utm_bounds[2])
-
max_y = math.ceil(utm_bounds[3])
-
else:
-
print("NoExtend")
-
print('UTM_Bounds: {}'.format(utm_bounds))
-
min_x = math.ceil(utm_bounds[0])
-
min_y = math.ceil(utm_bounds[1])
-
max_x = math.floor(utm_bounds[2])
-
max_y = math.floor(utm_bounds[3])
-
-
if quad_space:
-
print("quad_space")
-
row_cell = np.asarray([[0, 1], [2, 3]])
-
anchor_point_list_dict = {0: [], 1: [], 2: [], 3: []}
-
else:
-
anchor_point_list_dict = {0: []}
-
-
for rowidx, x in enumerate(np.arange(min_x, max_x, stride_size_meters)):
-
for colidx, y in enumerate(np.arange(min_y, max_y,
-
stride_size_meters)):
-
if quad_space:
-
anchor_point_list_dict[
-
row_cell[rowidx % 2, colidx % 2]].append([x, y])
-
else:
-
anchor_point_list_dict[0].append([x, y])
-
-
return anchor_point_list_dict
-
-
-[docs]def calculate_cells(anchor_point_list_dict, cell_size_meters, utm_bounds=[]):
-
""" Calculate boundaries for image cells (chips) from anchor points.
-
-
This function takes the output from :func:`calculate_anchor_points` as well
-
as a desired cell size (`cell_size_meters`) and outputs
-
``(W, S, E, N)`` tuples for generating cells.
-
-
Arguments
-
---------
-
anchor_point_list_dict : dict
-
Output of :func:`calculate_anchor_points`. See that function for
-
details.
-
cell_size_meters : int or float
-
Desired width and height of each cell in meters.
-
utm_bounds : list -like of float s, optional
-
A :obj:`list`-like of shape ``(W, S, E, N)`` that defines the limits
-
of an input image tile in UTM coordinates to ensure that no cells
-
extend beyond those limits. If not provided, all cells will be included
-
even if they extend beyond the UTM limits of the source imagery.
-
-
Returns
-
-------
-
cells_list_dict : dict of list(s) of lists
-
A dict whose keys are either ``0`` or ``[0, 1, 2, 3]`` (see
-
:func:`calculate_anchor_points` . ``quad_space`` ), and whose values are
-
:obj:`list` s of boundaries in the shape ``[W, S, E, N]`` . Boundaries
-
are in UTM coordinates.
-
-
"""
-
cells_list_dict = {}
-
for anchor_point_list_id, anchor_point_list in anchor_point_list_dict.items():
-
cells_list = []
-
for anchor_point in anchor_point_list:
-
if utm_bounds:
-
if (anchor_point[0] + cell_size_meters < utm_bounds[2]) and (
-
anchor_point[1] + cell_size_meters < utm_bounds[3]):
-
cells_list.append([anchor_point[0], anchor_point[1],
-
anchor_point[0] + cell_size_meters,
-
anchor_point[1] + cell_size_meters])
-
else:
-
pass
-
else: # include it regardless of extending beyond bounds
-
cells_list.append([anchor_point[0], anchor_point[1],
-
anchor_point[0] + cell_size_meters,
-
anchor_point[1] + cell_size_meters])
-
-
cells_list_dict[anchor_point_list_id] = cells_list
-
-
return cells_list_dict
-
-
-[docs]def calculate_analysis_grid(utm_bounds, stride_size_meters=300,
-
cell_size_meters=400, quad_space=False,
-
snapToGrid=False):
-
"""Wrapper for :func:`calculate_anchor_points` and :func:`calculate_cells`.
-
-
Based on UTM boundaries of an image tile, stride size, and cell size,
-
output a dictionary of boundary lists for analysis chips.
-
-
Arguments
-
---------
-
utm_bounds : list-like of shape ``(W, S, E, N)``
-
UTM coordinate limits of the input tile.
-
stride_size_meters : int, optional
-
Step size in both X and Y directions between cells in units of meters.
-
Defaults to ``300`` .
-
cell_size_meters : int, optional
-
Extent of each cell in both X and Y directions in units of meters.
-
Defaults to ``400`` .
-
quad_space : bool, optional
-
See :func:`calculate_anchor_points` . ``quad_space`` . Defaults to
-
``False`` .
-
snapToGrid : bool, optional
-
.. :deprecated: 0.2.0
-
This argument is deprecated and no longer has any effect.
-
-
Returns
-
-------
-
cells_list_dict : dict of list(s) of lists
-
A dict whose keys are either ``0`` or ``[0, 1, 2, 3]`` (see
-
:func:`calculate_anchor_points` . ``quad_space`` ), and whose values are
-
:obj:`list` s of boundaries in the shape ``[W, S, E, N]`` . Boundaries
-
are in UTM coordinates.
-
-
"""
-
anchor_point_list_dict = calculate_anchor_points(
-
utm_bounds, stride_size_meters=stride_size_meters,
-
quad_space=quad_space)
-
cells_list_dict = calculate_cells(anchor_point_list_dict, cell_size_meters,
-
utm_bounds=utm_bounds)
-
return cells_list_dict
-
-
-if __name__ == '__main__':
- utmX, utmY = 658029, 4006947
- cll_x = utmX
- cur_x = utmX + 500
- cll_y = utmY
- cur_y = utmY + 500
- stride_size_meters = 300
- cell_size_meters = 400
- ctile_size_pixels = 1600
- spacenetPath = "s3://spacenet-dataset/AOI_2_Vegas/srcData/rasterData/AOI_2_Vegas_MUL-PanSharpen_Cloud.tif"
- address = spacenetPath
-
- with rasterio.open(address) as src:
- cwgs_bounds = utils.get_wgs84_bounds(src)
- cutm_crs = utils.calculate_UTM_crs(cwgs_bounds)
- cutm_bounds = utils.get_utm_bounds(src, cutm_crs)
-
- #ccells_list = calculate_analysis_grid(cutm_bounds, stride_size_meters=stride_size_meters,
- # cell_size_meters=cell_size_meters)
-
- #random_cell = random.choice(ccells_list)
- #cll_x, cll_y, cur_x, cur_y = random_cell
- tile, mask, window, window_transform = tile_utm(src, cll_x, cll_y, cur_x, cur_y, indexes=None, tilesize=ctile_size_pixels, nodata=None, alpha=None,
- dst_crs=cutm_crs)
-