in libs/solaris/tile/raster_tile.py [0:0]
def tile_generator(self, src, dest_dir=None, channel_idxs=None,
nodata=None, alpha=None, aoi_boundary=None,
restrict_to_aoi=False):
"""Create the tiled output imagery from input tiles.
Uses the arguments provided at initialization to generate output tiles.
First, tile locations are generated based on `Tiler.tile_size` and
`Tiler.size_in_meters` given the bounds of the input image.
Arguments
---------
src : `str` or :class:`Rasterio.DatasetReader`
The source data to tile from. If this is a "classic"
(non-cloud-optimized) GeoTIFF, the whole image will be loaded in;
if it's cloud-optimized, only the required portions will be loaded
during tiling unless ``force_load_cog=True`` was specified upon
initialization.
dest_dir : str, optional
The path to the destination directory to output images to. If the
path doesn't exist, it will be created. This argument is required
if it wasn't provided during initialization.
channel_idxs : list, optional
The list of channel indices to be included in the output array.
If not provided, all channels will be included. *Note:* per
``rasterio`` convention, indexing starts at ``1``, not ``0``.
nodata : int, optional
The value in `src` that specifies nodata. If this value is not
provided, solaris will attempt to infer the nodata value from the
`src` metadata.
alpha : int, optional
The band to specify as alpha. If not provided, solaris will attempt
to infer if an alpha band is present from the `src` metadata.
aoi_boundary : `list`-like or :class:`shapely.geometry.Polygon`, optional
AOI bounds can be provided either as a
``[left, bottom, right, top]`` :class:`list`-like or as a
:class:`shapely.geometry.Polygon`.
restrict_to_aoi : bool, optional
Should output tiles be restricted to the limits of the AOI? If
``True``, any tile that partially extends beyond the limits of the
AOI will not be returned. This is the inverse of the ``boundless``
argument for :class:`rasterio.io.DatasetReader` 's ``.read()``
method.
Yields
------
tile_data, mask, tile_bounds
tile_data : :class:`numpy.ndarray`
A list of lists of each tile's bounds in the order they were
created, to be used in tiling vector data. These data are also
stored as an attribute of the :class:`Tiler` instance named
`tile_bounds`.
"""
# parse arguments
if self.verbose:
print("Checking input data...")
# if isinstance(src, str):
# self.is_cog = cog_validate(src)
# else:
# self.is_cog = cog_validate(src.name)
# if self.verbose:
# print('COG: {}'.format(self.is_cog))
self.src = _check_rasterio_im_load(src)
if channel_idxs is None: # if not provided, include them all
channel_idxs = list(range(1, self.src.count + 1))
print(channel_idxs)
self.src_crs = _check_crs(self.src.crs, return_rasterio=True) # necessary to use rasterio crs for reproject
if self.verbose:
print('Source CRS: EPSG:{}'.format(self.src_crs.to_epsg()))
if self.dest_crs is None:
self.dest_crs = self.src_crs
if self.verbose:
print('Destination CRS: EPSG:{}'.format(self.dest_crs.to_epsg()))
self.src_path = self.src.name
self.proj_unit = raster_get_projection_unit(self.src) # for rounding
if self.verbose:
print("Inputs OK.")
if self.use_src_metric_size:
if self.verbose:
print("Checking if inputs are in metric units...")
if self.project_to_meters:
if self.verbose:
print("Input CRS is not metric. "
"Reprojecting the input to UTM.")
self.src = reproject(self.src,
resampling_method=self.resampling,
dest_path=os.path.join(self.dest_dir,
'tmp.tif'))
if self.verbose:
print('Done reprojecting.')
if nodata is None and self.nodata is None:
self.nodata = self.src.nodata
else:
self.nodata = nodata
# get index of alpha channel
if alpha is None and self.alpha is None:
mf_list = [rasterio.enums.MaskFlags.alpha in i for i in
self.src.mask_flag_enums] # list with True at idx of alpha c
try:
self.alpha = np.where(mf_list)[0] + 1
except IndexError: # if there isn't a True
self.alpha = None
else:
self.alpha = alpha
if getattr(self, 'tile_bounds', None) is None:
self.get_tile_bounds()
for tb in self.tile_bounds:
# removing the following line until COG functionality implemented
if True: # not self.is_cog or self.force_load_cog:
window = rasterio.windows.from_bounds(
*tb, transform=self.src.transform,
width=self.src_tile_size[1],
height=self.src_tile_size[0])
if self.src.count != 1:
src_data = self.src.read(
window=window,
indexes=channel_idxs, boundless=True)
else:
src_data = self.src.read(
window=window,
boundless=True)
dst_transform, width, height = calculate_default_transform(
self.src.crs, self.dest_crs,
self.src.width, self.src.height, *tb,
dst_height=self.dest_tile_size[0],
dst_width=self.dest_tile_size[1])
if self.dest_crs != self.src_crs and self.resampling_method is not None:
tile_data = np.zeros(shape=(src_data.shape[0], height, width),
dtype=src_data.dtype)
rasterio.warp.reproject(
source=src_data,
destination=tile_data,
src_transform=self.src.window_transform(window),
src_crs=self.src.crs,
dst_transform=dst_transform,
dst_crs=self.dest_crs,
resampling=getattr(Resampling, self.resampling))
elif self.dest_crs != self.src_crs and self.resampling_method is None:
print("Warning: You've set resampling to None but your "
"destination projection differs from the source "
"projection. Using bilinear resampling by default.")
tile_data = np.zeros(shape=(src_data.shape[0], height, width),
dtype=src_data.dtype)
rasterio.warp.reproject(
source=src_data,
destination=tile_data,
src_transform=self.src.window_transform(window),
src_crs=self.src.crs,
dst_transform=dst_transform,
dst_crs=self.dest_crs,
resampling=getattr(Resampling, "bilinear"))
else: # for the case where there is no resampling and no dest_crs specified, no need to reproject or resample
tile_data = src_data
if self.nodata:
mask = np.all(tile_data != nodata,
axis=0).astype(np.uint8) * 255
elif self.alpha:
mask = self.src.read(self.alpha, window=window)
else:
mask = None # placeholder
# else:
# tile_data, mask, window, aff_xform = read_cog_tile(
# src=self.src,
# bounds=tb,
# tile_size=self.dest_tile_size,
# indexes=channel_idxs,
# nodata=self.nodata,
# resampling_method=self.resampling
# )
profile = self.src.profile
profile.update(width=self.dest_tile_size[1],
height=self.dest_tile_size[0],
crs=self.dest_crs,
transform=dst_transform)
if len(tile_data.shape) == 2: # if there's no channel band
profile.update(count=1)
else:
profile.update(count=tile_data.shape[0])
yield tile_data, mask, profile, tb