diff --git a/digger/make/reservoir.py b/digger/make/reservoir.py index 1698bfa217a1a31d41aceeaf504fa0f761bbf2f6..e98d0eaeefa3ff6b8a59bb5ad345386f2f4b5111 100644 --- a/digger/make/reservoir.py +++ b/digger/make/reservoir.py @@ -30,6 +30,7 @@ import numpy as np import rasterio from pysheds.grid import Grid from pysheds.sview import Raster, ViewFinder +from rasterio.warp import Resampling, reproject from digger import _plot from digger._utils import write @@ -46,6 +47,7 @@ def reservoir( q4_prefix: str = "q4", eta_prefix: str = "eta", dam_dict: dict = {}, + coarsen_factor: int = 1, write_tt3: bool = False, write_tif: bool = False, fig_path: str = "dam.png", @@ -76,7 +78,12 @@ def reservoir( - **y**: The y-coordinate of the downstream end of the dam - **volume**: The total volume of material behind the dam - **m0**: The solid volume fraction of the material behind the dam - + coarsen_factor : int + A factor by which to coarsen the DEM before further processing. This + may be useful in cases for which fine topography is available but not + needed for reservoir construction (the reservoir building algorithm + is slow for small grid cell sizes). Default is 1, which results in no + coarsening. write_tt3 : bool Whether to write out topotype 3 files of q1, eta, q4, and b. If written, the `*_prefix` for each output has '.tt3' added to the end. @@ -113,11 +120,51 @@ def reservoir( z = src.read(1) meta = src.meta transform = src.transform - metadata = { - "shape": z.shape, - "affine": transform, - "nodata": src.meta["nodata"] or np.nan, - } + + # Optionally coarsen the DEM. This is useful when you have a 1-m dem but don't + # need the DEM used for this step to be so fine. + # This algorithm really slows down between 10 and 3 m in typical usage. + if coarsen_factor < 1: + raise ValueError("coarsen_factor must be greater than 1") + + if coarsen_factor > 1: + + dst_row = int(np.floor(z.shape[0] / coarsen_factor)) + dst_col = int(np.floor(z.shape[1] / coarsen_factor)) + + # print(f"Coarsening from ({z.shape}) to ({dst_row}, {dst_col})") + dst_transform = transform * transform.scale(coarsen_factor) + + if np.abs(dst_transform[0]) > 20: + raise ValueError( + "coarsen_factor must not result in topography with grid cell sizes greater than 20 m" + ) + + dst_crs = meta["crs"] + z_warp = np.zeros((dst_row, dst_col)) + + reproject( + z, + z_warp, + src_transform=transform, + src_crs=meta["crs"], + dst_transform=dst_transform, + dst_crs=dst_crs, + resampling=Resampling.bilinear, + ) + + z = z_warp + transform = dst_transform + meta["width"] = dst_col + meta["height"] = dst_row + meta["transform"] = dst_transform + + # set metadata for rasterio + metadata = { + "shape": z.shape, + "affine": transform, + "nodata": src.meta["nodata"] or np.nan, + } view = ViewFinder(**metadata) dem = Raster(z, view) diff --git a/digger/make/scoops.py b/digger/make/scoops.py index a2423a185ee5f23b1aca73c858b8ad52d62fa10e..fedbd26c9df29eb7a195dc1f334e127f4920af74 100644 --- a/digger/make/scoops.py +++ b/digger/make/scoops.py @@ -95,7 +95,7 @@ Center, 2009, Prince William Sound, Alaska 8 arc-second MHHW coastal digital elevation model: National Oceanic and Atmospheric Administration [NOAA], National Centers for Environmental Information web page, `<https://data.noaa.gov/metaview/page?xml=NOAA/NESDIS/NGDC/MGG/DEM//iso/xml/ -638.xml&view=getDataView&header=none>`_. """ +638.xml&view=getDataView&header=none>`_.""" import warnings from typing import Dict diff --git a/digger/make/slabs.py b/digger/make/slabs.py index 93ff498e2d7f46870ea4ae850f4ca0f9072cc378..5c382acffd00e1c6dd2424f20c26705dfd849bfb 100644 --- a/digger/make/slabs.py +++ b/digger/make/slabs.py @@ -17,7 +17,7 @@ shapefile, the following information is expected in the attribute table The resulting output files support a number of additional options, including adding material of thickness h on top of the provided topography or subtracting -it from the provided topography. """ +it from the provided topography.""" import fiona import fiona.crs diff --git a/examples/pre-run/schultz_fire/make_reservoir.py b/examples/pre-run/schultz_fire/make_reservoir.py index c2b220aa25267737fbb4404d0689d2580142da0e..403fd3446c2951b1677887d62044e02c681f2d27 100644 --- a/examples/pre-run/schultz_fire/make_reservoir.py +++ b/examples/pre-run/schultz_fire/make_reservoir.py @@ -17,6 +17,7 @@ regions = make.reservoir( q1_prefix="q1_sch_reservoir", q4_prefix="q4_sch_reservoir", eta_prefix="eta_sch_reservoir", + coarsen_factor=1, dam_dict={"1": {"x": 444065, "y": 3909870, "volume": 100000, "m0": 0.61}}, write_tt3=True, write_tif=True, diff --git a/pyproject.toml b/pyproject.toml index 910a86bcf8af13dc5bb6e7fa63b1a26cf4a8d740..1b793ebe7e9d763d3fa7765b9764351ec2394a27 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.core.masonry.api" [tool.poetry] name = "digger" -version = "1.1.1" +version = "1.1.2" description = "digger: Utility tools for landslide runout modeling" authors = [ "Katherine R. Barnhart <krbarnhart@usgs.gov>", @@ -57,7 +57,7 @@ coverage = "*" pytest-cov = "*" pytest-datadir = "*" isort = "*" -black = "*" +black = "25.1.0" flake8 = "*" safety = "^2.1.1" sphinx = ">=7.2.6"