DEM Processing¶

In [1]:

from __future__ import annotations

import shutil
from pathlib import Path
from tempfile import TemporaryDirectory

import geopandas as gpd
import pywbt

import seamless_3dep as s3dep

from __future__ import annotations import shutil from pathlib import Path from tempfile import TemporaryDirectory import geopandas as gpd import pywbt import seamless_3dep as s3dep

Let's start by getting a HUC8 geometry from GeoConnex web service for St. Vrain region in Colorado. Note that we need to make sure that geometry is in 4326 projection.

In [2]:

url = "https://reference.geoconnex.us/collections/hu08/items/10190005"
vrain = gpd.read_file(url)
geom_org = vrain.to_crs(4326).union_all()
geom = vrain.to_crs(3857).buffer(5e3).to_crs(4326).union_all()

url = "https://reference.geoconnex.us/collections/hu08/items/10190005" vrain = gpd.read_file(url) geom_org = vrain.to_crs(4326).union_all() geom = vrain.to_crs(3857).buffer(5e3).to_crs(4326).union_all()

We can use get_dem to get the DEM but this function only returns DEM in 4326 projection. Since we want to analyze the DEM using WhiteboxTools and it's better to have the input DEM be in a projected CRS, instead, we use get_map to get the DEM in 3857 projection. Note that since we're getting the DEM in 3857 projection, it's a good idea to add some padding to the bounding box to make sure the obtained maps will not end up having NaN values within our original bounding box. In this case, we add a 5-km buffer to the bounding box.

In [3]:

data_dir = Path("data")
tiff_files = s3dep.get_map("DEM", geom.bounds, data_dir, 10)

data_dir = Path("data") tiff_files = s3dep.get_map("DEM", geom.bounds, data_dir, 10)

We then use build_vrt to build a VRT file from the obtained GeoTIFF files so we can use rioxarray to read the data into an xarray.DataArray.

In [4]:

dem = s3dep.tiffs_to_da(tiff_files, geom_org.bounds, crs=4326)
dem.size

dem = s3dep.tiffs_to_da(tiff_files, geom_org.bounds, crs=4326) dem.size

Out[4]:

35503248

In [5]:

ax = dem.plot.imshow(robust=True)
ax.figure.savefig("images/dem.png")

ax = dem.plot.imshow(robust=True) ax.figure.savefig("images/dem.png")

Now, we use PyWBT package to compute Topographic Wetness Index (TWI).

In [6]:

twi_files = [data_dir / fname.name.replace("dem", "twi") for fname in tiff_files]
slope_files = [data_dir / fname.name.replace("dem", "slope") for fname in tiff_files]
for dname, tname, sname in zip(tiff_files, twi_files, slope_files, strict=False):
    with TemporaryDirectory(dir=data_dir) as temp:
        shutil.copy(dname, temp)
        wbt_args = {
            "BreachDepressions": [f"-i={dname.name}", "--fill_pits", "-o=dem_corr.tiff"],
            "D8Pointer": ["-i=dem_corr.tiff", "-o=fdir.tiff"],
            "D8FlowAccumulation": [
                "-i=fdir.tiff",
                "--pntr",
                "--out_type='specific contributing area'",
                "-o=sca.tiff",
            ],
            "Slope": ["-i=dem_corr.tiff", "--units=degrees", f"-o={sname.name}"],
            "WetnessIndex": ["--sca=sca.tiff", f"--slope={sname.name}", f"-o={tname.name}"],
        }
        pywbt.whitebox_tools(temp, wbt_args, [sname.name, tname.name], temp)
        shutil.copy(Path(temp) / tname.name, data_dir)
        shutil.copy(Path(temp) / sname.name, data_dir)

twi_files = [data_dir / fname.name.replace("dem", "twi") for fname in tiff_files] slope_files = [data_dir / fname.name.replace("dem", "slope") for fname in tiff_files] for dname, tname, sname in zip(tiff_files, twi_files, slope_files, strict=False): with TemporaryDirectory(dir=data_dir) as temp: shutil.copy(dname, temp) wbt_args = { "BreachDepressions": [f"-i={dname.name}", "--fill_pits", "-o=dem_corr.tiff"], "D8Pointer": ["-i=dem_corr.tiff", "-o=fdir.tiff"], "D8FlowAccumulation": [ "-i=fdir.tiff", "--pntr", "--out_type='specific contributing area'", "-o=sca.tiff", ], "Slope": ["-i=dem_corr.tiff", "--units=degrees", f"-o={sname.name}"], "WetnessIndex": ["--sca=sca.tiff", f"--slope={sname.name}", f"-o={tname.name}"], } pywbt.whitebox_tools(temp, wbt_args, [sname.name, tname.name], temp) shutil.copy(Path(temp) / tname.name, data_dir) shutil.copy(Path(temp) / sname.name, data_dir)

In [7]:

twi = s3dep.tiffs_to_da(twi_files, geom_org.bounds, crs=4326)
ax = twi.plot.imshow(robust=True)
ax.figure.savefig("images/twi.png")

twi = s3dep.tiffs_to_da(twi_files, geom_org.bounds, crs=4326) ax = twi.plot.imshow(robust=True) ax.figure.savefig("images/twi.png")

In [8]:

slope = s3dep.tiffs_to_da(slope_files, geom_org.bounds, crs=4326)
ax = slope.plot.imshow(robust=True)
ax.figure.savefig("images/slope.png")

slope = s3dep.tiffs_to_da(slope_files, geom_org.bounds, crs=4326) ax = slope.plot.imshow(robust=True) ax.figure.savefig("images/slope.png")

We can directly get slope using get_map function and passing map_type="Slope Degrees"

In [9]:

tiff_files = s3dep.get_map("Slope Degrees", geom.bounds, data_dir, 10)
slope = s3dep.tiffs_to_da(tiff_files, geom_org.bounds, crs=4326)
ax = slope.plot.imshow(robust=True)
ax.figure.savefig("images/slope_dynamic.png")

tiff_files = s3dep.get_map("Slope Degrees", geom.bounds, data_dir, 10) slope = s3dep.tiffs_to_da(tiff_files, geom_org.bounds, crs=4326) ax = slope.plot.imshow(robust=True) ax.figure.savefig("images/slope_dynamic.png")