geowatch.utils.util_raster module

Utilities for rasterio

SeeAlso

util_gdal.py

geowatch.utils.util_raster.mask(raster: DatasetReader | str, default_nodata=None, save=False, convex_hull=False, as_poly=True, tolerance=None, max_polys=None, use_overview=0)[source]

Compute a raster’s valid data mask in pixel coordinates.

Note that this is the rasterio mask, which for multi-band rasters is the binary OR of the individual band masks. This is different from the gdal mask, which is always per-band.

Parameters:

  • raster (str) – Path to a dataset (raster image file)

  • default_nodata (int) – if raster’s nodata value is None, default to this

  • save (bool) – if True and raster’s nodata value is None, write the default to it. If False, performance overhead is incurred from creating a tempfile

  • convex_hull (bool) – if True, return the convex hull of the mask image or poly

  • as_poly (bool) – if True, return the mask as a shapely Polygon or MultiPolygon instead of a raster image, in (w, h) order (opposite of Python convention).

  • tolerance (int) – if specified, simplifies the valid polygon.

  • use_overview (int) – if non-zero uses the closest overview if it is available. This increases computation time, but gives a better polygon when use_overview is closer to 0. Note, the polygon is rescaled to ensure it is returned in the input pixel space, not the overview space.

Returns:

If as_poly, a shapely Polygon or MultiPolygon bounding the valid data region(s) in pixel coordinates.

Else, a uint8 raster mask of the same shape as the input, where 255 == valid and 0 == invalid.

Example

>>> # xdoctest: +REQUIRES(--network)
>>> from geowatch.utils.util_raster import *
>>> # FIXME; this demo path no longer has any nodata values
>>> # Find a better demo with nodata
>>> from geowatch.demo.landsat_demodata import grab_landsat_product
>>> path = grab_landsat_product()['bands'][0]
>>> #
>>> mask_img = mask(path, as_poly=False)
>>> import kwimage as ki
>>> assert mask_img.shape == ki.load_image_shape(path)[:2]
>>> got = set(np.unique(mask_img))
>>> print(f'got={got}')
>>> # assert got == {0, 255}  # cant do this until nodata is fixed
>>> #
>>> mask_poly = mask(path, as_poly=True)
>>> import shapely
>>> assert isinstance(mask_poly, shapely.geometry.Polygon)
>>> # xdoctest: +REQUIRES(--show)
>>> import kwplot
>>> kwplot.autompl()
>>> kwplot.figure(fnum=1, doclf=True)
>>> imdata = kwimage.imread(path)
>>> imdata = kwimage.normalize_intensity(imdata)
>>> kw_poly = kwimage.Polygon.coerce(mask_poly.buffer(0).simplify(10))
>>> canvas = imdata.copy()
>>> mask_alpha = kwimage.ensure_alpha_channel(mask_img, alpha=(mask_img > 0))
>>> canvas = kwimage.overlay_alpha_layers([mask_alpha, canvas])
>>> canvas = kw_poly.scale(0.9, about='center').draw_on(canvas, color='green', alpha=0.6)
>>> kw_poly.scale(1.1, about='center').draw(alpha=0.5, color='red', setlim=True)
>>> kwplot.imshow(canvas)

Example

>>> # Test how the "save" functionality modifies the data
>>> import kwimage
>>> from geowatch.utils.util_raster import *
>>> import pathlib
>>> dpath = ub.Path.appdir('geowatch/tests/empty_raster').ensuredir()
>>> raster = dpath / 'empty.tif'
>>> ub.delete(raster)
>>> kwimage.imwrite(raster, np.zeros((3, 3, 5)))
>>> info1 = ub.cmd('gdalinfo {}'.format(raster))
>>> nodata = 0
>>> mask_img = mask(raster, as_poly=False)
>>> print('mask_img = {!r}'.format(mask_img))
>>> info2 = ub.cmd('gdalinfo {}'.format(raster))
>>> mask_poly = mask(raster, as_poly=True)
>>> info3 = ub.cmd('gdalinfo {}'.format(raster))
>>> print(info1['out'])
>>> print(info2['out'])
>>> print(info3['out'])
class geowatch.utils.util_raster.ResampledRaster(raster: str | DatasetReader, scale: float = 2, read: bool = True, resampling: Resampling = Resampling.bilinear)[source]

Bases: ExitStack

Context manager to rescale a raster on the fly using rasterio

This changes the number of pixels in the raster while maintaining its geographic bounds, that is, it changes the raster’s GSD.

Parameters:

  • raster – a DatasetReader (the object returned by rasterio.open) or path to a dataset

  • scale – factor to upscale the resolution, aka downscale the GSD, by

  • read – if True, read and return the resampled data (an expensive operation if scale>1) else, return the resampled dataset’s .profile attribute (metadata)

  • resampling – resampling algorithm, from rasterio.enums.Resampling [1]

Example

>>> # xdoctest: +REQUIRES(--slow)
>>> # xdoctest: +REQUIRES(--network)
>>> from geowatch.utils.util_raster import *
>>> from geowatch.demo.landsat_demodata import grab_landsat_product
>>> path = grab_landsat_product()['bands'][0]
>>> #
>>> current_gsd_meters = 60
>>> desired_gsd_meters = 10
>>> scale = current_gsd_meters / desired_gsd_meters
>>> #
>>> with rasterio.open(path) as f:
>>>     old_profile = f.profile
>>> #
>>> # can instantiate this class in a with-block
>>> with ResampledRaster(path, scale=scale, read=False) as f:
>>>     pass
>>> #
>>> # or have it stick around and change the resampling on the fly
>>> resampled = ResampledRaster(path, scale=scale, read=False)
>>> #
>>> # the computation only happens when you invoke 'with'
>>> with resampled as new_profile:
>>>     assert new_profile['width'] == int(old_profile['width'] * scale)
>>>     assert new_profile['crs'] == old_profile['crs']
>>> #
>>> resampled.scale = scale / 2
>>> resampled.read = True
>>> #
>>> with resampled as new:
>>>     assert new.profile['width'] == int(old_profile['width'] * scale / 2)
>>>     assert new.profile['crs'] == old_profile['crs']
>>>     # do other stuff with new

References

https://gis.stackexchange.com/a/329439 https://rasterio.readthedocs.io/en/latest/topics/reading.html https://rasterio.readthedocs.io/en/latest/topics/profiles.html [1] https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling

raster: str | DatasetReader
scale: float = 2
read: bool = True
resampling: Resampling = 1