add raster clipping options

[jhbravo]

Hello @mdbartos pysheds is a great package, thanks! it is possible to add clipping ...? by extension [lonmin, lonmax, latmin, latmax] or by mask (a shapefile)

[mdbartos]

Thank you, great suggestion. Pysheds does offer the ability to do clipping and other spatial referencing operations through the pysheds.view module---though it could use some more helper functions for convenience. These documentation pages might be helpful:

https://mattbartos.com/pysheds/raster.html https://mattbartos.com/pysheds/views.html

Essentially, you can do any clipping or masking operation you want by using a ViewFinder object. I will post some examples below in a followup comment.

[mdbartos]

Here's an example showing how to use 'views' to perform clipping and alignment of different datasets. This example uses the dem.tif file in the ./data directory, and the impervious_area.tiff file linked in the README. Note that these two datasets have differing extents.

# Import modules
import numpy as np
from pysheds.grid import Grid
import matplotlib.pyplot as plt

# Load data
grid = Grid.from_raster('dem.tif')
dem = grid.read_raster('dem.tif')
imperv = grid.read_raster('impervious_area.tiff')

# Spatial reference of dem dataset
dem.viewfinder

Output...

``` 'affine' : Affine(0.0008333333333333, 0.0, -97.4849999999961, 0.0, -0.0008333333333333, 32.82166666666536) 'shape' : (359, 367) 'nodata' : -32768 'crs' : Proj('+proj=longlat +datum=WGS84 +no_defs', preserve_units=True) 'mask' : array([[ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], ..., [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True]]) ```

# Spatial reference of impervious area dataset
imperv.viewfinder

Output...

``` 'affine' : Affine(30.0, 0.0, -167475.0, 0.0, -30.0, 1155045.0) 'shape' : (4471, 5469) 'nodata' : 0 'crs' : Proj('+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs', preserve_units=True) 'mask' : array([[ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], ..., [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True]]) ```

# Plot dem data
plt.imshow(dem, cmap='terrain', extent=dem.extent)
plt.title('DEM')

# Plot impervious area data
plt.imshow(imperv, cmap='bone', extent=imperv.extent)
plt.title('Impervious area')

# Convert CRS of impervious area data to DEM CRS
imperv_p = imperv.to_crs(dem.crs)
# Plot projected impervious area data
plt.imshow(imperv_p, cmap='bone', extent=imperv_p.extent)
plt.title('Impervious area (projected)')

# Set viewfinder of grid to DEM's viewfinder
grid.viewfinder = dem.viewfinder
# Create view of impervious area at DEM's spatial extent
imperv_view = grid.view(imperv_p)
# Plot view of impervious area
plt.imshow(imperv_view, cmap='bone', extent=imperv_view.extent)
plt.title('Impervious area (at DEM extent)')

[mdbartos]

The following example extends the example above by showing how a custom viewfinder can be created to provide a view corresponding to a given bounding box:

# Import modules
from affine import Affine
from pysheds.view import ViewFinder

# Create a custom viewfinder based on a given bbox
xmin, xmax, ymin, ymax = -97., -96.4, 32.6, 33.2
affine = Affine(grid.affine.a, grid.affine.b, xmin,
                grid.affine.d, grid.affine.e, ymax)
cols = int((xmax - xmin) // abs(affine.a))
rows = int((ymax - ymin) // abs(affine.e))
viewfinder = ViewFinder(affine=affine, shape=(rows, cols), crs=grid.crs)

viewfinder

Output...

``` 'affine' : Affine(0.0008333333333333, 0.0, -97.0, 0.0, -0.0008333333333333, 33.2) 'shape' : (720, 720) 'nodata' : 0 'crs' : Proj('+proj=longlat +datum=WGS84 +no_defs', preserve_units=True) 'mask' : array([[ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], ..., [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True]]) ```

# Set viewfinder of grid to custom viewfinder
grid.viewfinder = viewfinder
# View impervious area with new viewfinder
new_imperv_view = grid.view(imperv_p)
# Plot view of impervious area
plt.imshow(new_imperv_view, cmap='bone', extent=new_imperv_view.extent)
plt.title('Impervious area (new view)')

More info about the affine module and its format can be found here: https://www.perrygeo.com/python-affine-transforms.html

[mdbartos]

Feel free to suggest any particular helper functions that you think would make this process easier for you.

MDB