From bf18986715472f7f072338af5f54616351e9636b Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Wed, 23 Oct 2024 16:54:11 -0400 Subject: [PATCH] Update whitebox notebook (#126) * Update whitebox notebook * Update the whitebox notebook * Add py3dep * Add pyvista * Remove pyvista * Set quiet download --- .github/workflows/build.yml | 4 +- .github/workflows/deploy.yml | 4 +- .gitignore | 8 + _config.yml | 1 - book/geospatial/whitebox.ipynb | 888 +++++++++++++++++++++++++-------- book/geospatial/whitebox.md | 434 ++++++++++++---- requirements.txt | 1 + 7 files changed, 1032 insertions(+), 308 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 7926983..6fe4d54 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -28,7 +28,9 @@ jobs: run: uv python install ${{ matrix.python-version }} - name: Install dependencies - run: uv sync --python ${{ matrix.python-version }} + run: | + uv sync --python ${{ matrix.python-version }} + uv pip install py3dep - name: Install optional dependencies run: | diff --git a/.github/workflows/deploy.yml b/.github/workflows/deploy.yml index 849322e..bc264d0 100644 --- a/.github/workflows/deploy.yml +++ b/.github/workflows/deploy.yml @@ -28,7 +28,9 @@ jobs: run: uv python install ${{ matrix.python-version }} - name: Install dependencies - run: uv sync --python ${{ matrix.python-version }} + run: | + uv sync --python ${{ matrix.python-version }} + uv pip install py3dep - name: Install optional dependencies run: | diff --git a/.gitignore b/.gitignore index 6a5d877..a2efcd7 100644 --- a/.gitignore +++ b/.gitignore @@ -37,3 +37,11 @@ book/geospatial/nyc_boroughs.* # pixi environments .pixi *.egg-info +*.zip +*.las +longest_flowpath.* +pour_point* +streams* +watershed* +histogram* +book/geospatial/cache/* \ No newline at end of file diff --git a/_config.yml b/_config.yml index 22c5e92..401399f 100644 --- a/_config.yml +++ b/_config.yml @@ -23,7 +23,6 @@ execute: "book/geospatial/vector_viz.ipynb", "book/geospatial/raster_viz.ipynb", "book/geospatial/maplibre.ipynb", - "book/geospatial/whitebox.ipynb", "book/geospatial/geemap.ipynb", "book/geospatial/samgeo.ipynb", "book/geospatial/hypercoast.ipynb", diff --git a/book/geospatial/whitebox.ipynb b/book/geospatial/whitebox.ipynb index 6f393fb..d289c1e 100644 --- a/book/geospatial/whitebox.ipynb +++ b/book/geospatial/whitebox.ipynb @@ -5,7 +5,7 @@ "id": "0", "metadata": {}, "source": [ - "# WhiteboxTools\n", + "# Whitebox\n", "\n", "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/giswqs/geog-312/blob/main/book/geospatial/whitebox.ipynb)" ] @@ -15,9 +15,62 @@ "id": "1", "metadata": {}, "source": [ + "## Overview\n", + "\n", + "In this lecture, we will explore the use of [**WhiteboxTools**](https://github.com/jblindsay/whitebox-tools), a powerful open-source library for performing geospatial analysis. Specifically, we will focus on two key applications: **watershed analysis** and **LiDAR data analysis**. You will learn how to manipulate geospatial data using Python, conduct hydrological analysis, and derive digital elevation models (DEMs) and canopy height models (CHMs) from LiDAR data.\n", + "\n", + "This lecture is structured into two main sections:\n", + "1. **Watershed Analysis**: Using DEMs and hydrological tools to delineate watersheds, calculate flow accumulation, and extract stream networks.\n", + "2. **LiDAR Data Analysis**: Processing LiDAR point cloud data to derive DEMs, DSMs, and CHMs while removing outliers and improving data quality.\n", + "\n", + "By the end of this session, you will have gained hands-on experience with **WhiteboxTools** and **leafmap**, allowing you to perform a wide range of geospatial and hydrological analyses.\n", + "\n", + "## Learning Objectives\n", + "\n", + "By the end of this lecture, you will be able to:\n", + "- Install and configure **WhiteboxTools** and **leafmap** for geospatial analysis.\n", + "- Create interactive maps to visualize basemaps and geospatial datasets.\n", + "- Perform watershed analysis by delineating watersheds, flow directions, and stream networks.\n", + "- Manipulate and analyze Digital Elevation Models (DEMs) to conduct hydrological modeling.\n", + "- Process and analyze LiDAR data to generate **Digital Surface Models (DSMs)**, **Digital Elevation Models (DEMs)**, and **Canopy Height Models (CHMs)**.\n", + "- Integrate **WhiteboxTools** with Python workflows to automate geospatial analysis.\n", + "\n", + "\n", + "## Introduction\n", + "\n", + "Below is a brief introduction to Whitebox, a powerful open-source library for geospatial analysis. For more information, please refer to the whiteboxgeo website: https://www.whiteboxgeo.com.\n", + "\n", + "### What is Whitebox?\n", + "\n", + "Whitebox is geospatial data analysis software originally developed at the University of Guelph‘s Geomorphometry and Hydrogeomatics Research Group ([GHRG](https://jblindsay.github.io/ghrg/index.html)) directed by Dr. John Lindsay. Whitebox contains over 550 tools for processing many types of geospatial data. It has many great features such as its extensive use of parallel computing, it doesn’t need other libraries to be installed (e.g., GDAL), it can be used from scripting environments, and it easily plugs into other geographical information system (GIS) software. The Whitebox Toolset Extension provides even more power.\n", + "\n", + "### What can Whitebox do?\n", + "\n", + "Whitebox can be used to perform common GIS and remote sensing analysis tasks. Whitebox also contains advanced tooling for spatial hydrological analysis and LiDAR data processing. Whitebox is not a cartographic or spatial data visualization package; instead it’s meant to serve as an analytical back-end for other data visualization software, like QGIS and ArcGIS. \n", + "\n", + "### How is Whitebox different?\n", + "\n", + "Whitebox doesn’t compete with QGIS, ArcGIS/Pro, and ArcPy but rather it extends them. You can plug WhiteboxTools into QGIS and ArcGIS and it’ll provide hundreds of additional tools for analyzing all kinds of geospatial data. You can also call Whitebox functions from Python scripts using [Whitebox Workflows](https://www.whiteboxgeo.com/whitebox-workflows-for-python) (WbW). Combine WbW with ArcPy to more effectively automate your data analysis workflows and streamline your geoprocessing solutions.\n", + "\n", + "There are many tools in Whitebox that you won’t find elsewhere. You can think of Whitebox as a portable, cross-platform GIS analysis powerhouse, allowing you to extend your preferred GIS or to embed Whitebox capabilities into your automated scripted workflows. Oh, and it’s fast, really fast!\n", + "\n", + "\n", + "## Useful Resources for Whitebox\n", + "\n", + "- GitHub Repository: https://github.com/jblindsay/whitebox-tools\n", + "- WhiteboxGeo: https://www.whiteboxgeo.com\n", + "- User Manual: https://whiteboxgeo.com/manual/wbt_book/preface.html\n", + "- Whitebox Workflows for Python: https://www.whiteboxgeo.com/whitebox-workflows-for-python\n", + "- Whitebox Workflows for Python Pro: https://www.whiteboxgeo.com/whitebox-workflows-professional\n", + "- Python Frontend: https://github.com/opengeos/whitebox-python\n", + "- Jupyter Frontend: https://github.com/opengeos/whiteboxgui\n", + "- R Frontend: https://github.com/opengeos/whiteboxR\n", + "- ArcGIS Toolbox: https://github.com/opengeos/WhiteboxTools-ArcGIS\n", + "- QGIS Plugin: https://plugins.qgis.org/plugins/whitebox_workflows_for_qgis\n", + "\n", "## Installation\n", "\n", - "Uncomment and run the following cell to install necessary packages for this notebook, including leafmap, geopandas, localtileserver, rio-cogeo, pynhd, py3dep." + "To get started, we need to install the required packages, such as `leafmap`, and `whitebox`. Uncomment the following lines to install the packages." ] }, { @@ -27,13 +80,23 @@ "metadata": {}, "outputs": [], "source": [ - "# %pip install \"leafmap[raster]\" geopandas" + "# %pip install \"leafmap[raster]\" \"leafmap[lidar]\" mapclassify" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "id": "3", "metadata": {}, + "outputs": [], + "source": [ + "# %pip install numpy==1.26.4" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, "source": [ "## Import libraries" ] @@ -41,49 +104,49 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "5", "metadata": {}, "outputs": [], "source": [ "import os\n", - "import leafmap" + "import leafmap\n", + "import numpy as np" ] }, { "cell_type": "markdown", - "id": "5", + "id": "6", "metadata": {}, "source": [ - "## Create interactive maps\n", - "\n", - "Specify the map center, zoom level, and height." + "Some of the raster datasets generated by whitebox will be int32 type with a nodata value of -32768. To make it easier to visualize these datasets, we set the nodata value as an environment variable, which will be used by leafmap to set the nodata value when visualizing the raster datasets." ] }, { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "7", "metadata": {}, "outputs": [], "source": [ - "m = leafmap.Map(center=[40, -100], zoom=4, height=\"600px\")\n", - "m" + "os.environ[\"NODATA\"] = \"-32768\"" ] }, { "cell_type": "markdown", - "id": "7", + "id": "8", "metadata": {}, "source": [ - "## Add basemaps\n", + "## Part 1: Watershed Analysis\n", + "\n", + "### Create Interactive Maps\n", "\n", - "Add OpenTopoMap, USGS 3DEP Elevation, and USGS Hydrography basemaps." + "To perform watershed analysis, we first create an interactive map using **leafmap**. This step allows us to visualize different basemaps." ] }, { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -96,63 +159,45 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "10", "metadata": {}, "source": [ - "Add NLCD land cover map and legend." + "### Download Watershed Data\n", + "\n", + "Next, we download watershed data for the **Calapooia River basin** in Oregon. We'll use the latitude and longitude of a point in the basin to extract watershed boundary data." ] }, { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "11", "metadata": {}, "outputs": [], "source": [ - "m = leafmap.Map(center=[40, -100], zoom=4)\n", - "m.add_basemap(\"HYBRID\")\n", - "m.add_basemap(\"NLCD 2019 CONUS Land Cover\")\n", - "m.add_legend(builtin_legend=\"NLCD\", title=\"NLCD Land Cover Type\")\n", + "lat = 44.361169\n", + "lon = -122.821802\n", + "\n", + "m = leafmap.Map(center=[lat, lon], zoom=10)\n", + "m.add_marker([lat, lon])\n", "m" ] }, { "cell_type": "markdown", - "id": "11", + "id": "12", "metadata": {}, "source": [ - "Add WMS layers." + "Download the watershed data and visualize it:" ] }, { "cell_type": "code", "execution_count": null, - "id": "12", - "metadata": {}, - "outputs": [], - "source": [ - "m = leafmap.Map(center=[40, -100], zoom=4)\n", - "m.add_basemap(\"Esri.WorldImagery\")\n", - "url = \"https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?\"\n", - "m.add_wms_layer(\n", - " url,\n", - " layers=\"NLCD_2019_Land_Cover_L48\",\n", - " name=\"NLCD 2019 CONUS Land Cover\",\n", - " format=\"image/png\",\n", - " transparent=True,\n", - ")\n", - "m.add_legend(builtin_legend=\"NLCD\", title=\"NLCD Land Cover Type\")\n", - "m" - ] - }, - { - "cell_type": "markdown", "id": "13", "metadata": {}, + "outputs": [], "source": [ - "## Get watershed data\n", - "\n", - "Let's download watershed data for the Calapooia River basin in Oregon." + "geometry = {\"x\": lon, \"y\": lat}" ] }, { @@ -162,7 +207,8 @@ "metadata": {}, "outputs": [], "source": [ - "gdf = leafmap.get_nhd_basins(feature_ids=23763529, fsource=\"comid\", simplified=False)" + "gdf = leafmap.get_wbd(geometry, digit=10, return_geometry=True)\n", + "gdf.explore()" ] }, { @@ -170,7 +216,7 @@ "id": "15", "metadata": {}, "source": [ - "Plot the watershed boundary on the map." + "Save the watershed boundary to a **GeoJSON** file:" ] }, { @@ -180,9 +226,7 @@ "metadata": {}, "outputs": [], "source": [ - "m = leafmap.Map()\n", - "m.add_gdf(gdf, layer_name=\"Catchment\", info_mode=None)\n", - "m" + "gdf.to_file(\"basin.geojson\")" ] }, { @@ -190,82 +234,62 @@ "id": "17", "metadata": {}, "source": [ - "Save the watershed boundary to a GeoJSON or shapefile." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "18", - "metadata": {}, - "outputs": [], - "source": [ - "gdf.to_file(\"basin.geojson\", driver=\"GeoJSON\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "19", - "metadata": {}, - "outputs": [], - "source": [ - "gdf.to_file(\"basin.shp\")" - ] - }, - { - "cell_type": "markdown", - "id": "20", - "metadata": {}, - "source": [ - "## Download DEM\n", + "### Download and Display DEM\n", "\n", - "Download a digital elevation model (DEM) for the watershed from the USGS 3DEP Elevation service. Convert the DEM to a Cloud Optimized GeoTIFF (COG)." + "We download a **Digital Elevation Model (DEM)** from the USGS 3DEP Elevation service to analyze the terrain of the watershed. The DEM will be used to delineate watersheds, calculate flow accumulation, and extract stream networks. The `leafmap.get_3dep_dem()` function returns the DEM as an `xarray.DataArray` object. Optionally, you can save the DEM to a GeoTIFF file by setting the `output` parameter." ] }, { "cell_type": "code", "execution_count": null, - "id": "21", + "id": "18", "metadata": {}, "outputs": [], "source": [ - "leafmap.get_3dep_dem(\n", - " gdf, resolution=30, output=\"dem.tif\", dst_crs=\"EPSG:3857\", to_cog=True\n", - ")" + "array = leafmap.get_3dep_dem(\n", + " gdf,\n", + " resolution=30,\n", + " output=\"dem.tif\",\n", + " dst_crs=\"EPSG:3857\",\n", + " to_cog=True,\n", + " overwrite=True,\n", + ")\n", + "array" ] }, { "cell_type": "markdown", - "id": "22", + "id": "19", "metadata": {}, "source": [ - "Display the DEM on the map." + "Visualize the DEM on the map:" ] }, { "cell_type": "code", "execution_count": null, - "id": "23", + "id": "20", "metadata": {}, "outputs": [], "source": [ - "m.add_raster(\"dem.tif\", palette=\"terrain\", layer_name=\"DEM\")\n", + "m.add_raster(\"dem.tif\", palette=\"terrain\", nodata=np.nan, layer_name=\"DEM\")\n", "m" ] }, { "cell_type": "markdown", - "id": "24", + "id": "21", "metadata": {}, "source": [ - "## Get DEM metadata" + "### Get DEM metadata\n", + "\n", + "We can get the metadata of the DEM, such as the spatial resolution, bounding box, and coordinate reference system (CRS)." ] }, { "cell_type": "code", "execution_count": null, - "id": "25", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -275,47 +299,48 @@ }, { "cell_type": "markdown", - "id": "26", + "id": "23", "metadata": {}, "source": [ "Get a summary statistics of the DEM." ] }, { - "cell_type": "code", - "execution_count": null, - "id": "27", + "cell_type": "markdown", + "id": "24", "metadata": {}, - "outputs": [], "source": [ - "metadata[\"bands\"]" + "### Add colorbar\n", + "\n", + "Add a colorbar to the map to show the elevation values. Use the `image_min_max()` function to get the minimum and maximum values of the DEM." ] }, { - "cell_type": "markdown", - "id": "28", + "cell_type": "code", + "execution_count": null, + "id": "25", "metadata": {}, + "outputs": [], "source": [ - "## Add colorbar" + "leafmap.image_min_max(\"dem.tif\")" ] }, { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "26", "metadata": {}, "outputs": [], "source": [ - "m.add_colormap(cmap=\"terrain\", vmin=\"60\", vmax=1500, label=\"Elevation (m)\")\n", - "m" + "m.add_colormap(cmap=\"terrain\", vmin=\"60\", vmax=1500, label=\"Elevation (m)\")" ] }, { "cell_type": "markdown", - "id": "30", + "id": "27", "metadata": {}, "source": [ - "## Initialize WhiteboxTools\n", + "### Initialize WhiteboxTools\n", "\n", "Initialize the WhiteboxTools class." ] @@ -323,7 +348,7 @@ { "cell_type": "code", "execution_count": null, - "id": "31", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -332,7 +357,7 @@ }, { "cell_type": "markdown", - "id": "32", + "id": "29", "metadata": {}, "source": [ "Check the WhiteboxTools version." @@ -341,7 +366,7 @@ { "cell_type": "code", "execution_count": null, - "id": "33", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -350,16 +375,16 @@ }, { "cell_type": "markdown", - "id": "34", + "id": "31", "metadata": {}, "source": [ - "Display the WhiteboxTools interface." + "Display the WhiteboxTools interface, which lists all available tools. You can use this interface to search for specific tools. You can also run any tool using the interface. However, we will use the Python API to run the tools." ] }, { "cell_type": "code", "execution_count": null, - "id": "35", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -368,16 +393,18 @@ }, { "cell_type": "markdown", - "id": "36", + "id": "33", "metadata": {}, "source": [ - "## Set working directory" + "### Set working directory\n", + "\n", + "Set the working directory to save the intermediate and output files. Set `wbt.version=False` to suppress the Whitebox processing log." ] }, { "cell_type": "code", "execution_count": null, - "id": "37", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -387,18 +414,20 @@ }, { "cell_type": "markdown", - "id": "38", + "id": "35", "metadata": {}, "source": [ - "## Smooth DEM\n", + "### Smooth DEM\n", "\n", - "All WhiteboxTools functions will return 0 if they are successful, and 1 if they are not." + "All WhiteboxTools functions will return 0 if they are successful, and 1 if they are not.\n", + "\n", + "Smoothing the DEM enhances the quality of subsequent hydrological analysis." ] }, { "cell_type": "code", "execution_count": null, - "id": "39", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -407,37 +436,49 @@ }, { "cell_type": "markdown", - "id": "40", + "id": "37", + "metadata": {}, + "source": [ + "Visualize the smoothed DEM and watershed boundary on the map." + ] + }, + { + "cell_type": "markdown", + "id": "38", "metadata": {}, "source": [ - "Display the smoothed DEM and watershed boundary on the map." + "Visualize the smoothed DEM:" ] }, { "cell_type": "code", "execution_count": null, - "id": "41", + "id": "39", "metadata": {}, "outputs": [], "source": [ "m = leafmap.Map()\n", - "m.add_raster(\"smoothed.tif\", palette=\"terrain\", layer_name=\"Smoothed DEM\")\n", + "m.add_basemap(\"Satellite\")\n", + "m.add_raster(\"smoothed.tif\", colormap=\"terrain\", layer_name=\"Smoothed DEM\")\n", "m.add_geojson(\"basin.geojson\", layer_name=\"Watershed\", info_mode=None)\n", + "m.add_basemap(\"USGS Hydrography\", show=False)\n", "m" ] }, { "cell_type": "markdown", - "id": "42", + "id": "40", "metadata": {}, "source": [ - "## Create hillshade" + "### Create hillshade\n", + "\n", + "Create a hillshade from the smoothed DEM." ] }, { "cell_type": "code", "execution_count": null, - "id": "43", + "id": "41", "metadata": {}, "outputs": [], "source": [ @@ -446,7 +487,7 @@ }, { "cell_type": "markdown", - "id": "44", + "id": "42", "metadata": {}, "source": [ "Overlay the hillshade on the smoothed DEM with transparency." @@ -455,7 +496,7 @@ { "cell_type": "code", "execution_count": null, - "id": "45", + "id": "43", "metadata": {}, "outputs": [], "source": [ @@ -465,18 +506,18 @@ }, { "cell_type": "markdown", - "id": "46", + "id": "44", "metadata": {}, "source": [ - "## Find no-flow cells\n", + "### Find no-flow cells\n", "\n", - "Find cells with undefined flow, i.e. no valid flow direction, based on the D8 flow direction algorithm" + "Find cells with undefined flow, i.e. no valid flow direction, based on the D8 flow direction algorithm." ] }, { "cell_type": "code", "execution_count": null, - "id": "47", + "id": "45", "metadata": {}, "outputs": [], "source": [ @@ -485,35 +526,36 @@ }, { "cell_type": "markdown", - "id": "48", + "id": "46", "metadata": {}, "source": [ - "Display the no-flow cells on the map." + "Visualize the no-flow cells on the map." ] }, { "cell_type": "code", "execution_count": null, - "id": "49", + "id": "47", "metadata": {}, "outputs": [], "source": [ - "m.add_raster(\"noflow.tif\", layer_name=\"No Flow Cells\")\n", - "m" + "m.add_raster(\"noflow.tif\", layer_name=\"No Flow Cells\")" ] }, { "cell_type": "markdown", - "id": "50", + "id": "48", "metadata": {}, "source": [ - "## Fill depressions" + "### Fill depressions\n", + "\n", + "First, we fill any depressions in the DEM to ensure proper flow simulation." ] }, { "cell_type": "code", "execution_count": null, - "id": "51", + "id": "49", "metadata": {}, "outputs": [], "source": [ @@ -522,7 +564,7 @@ }, { "cell_type": "markdown", - "id": "52", + "id": "50", "metadata": {}, "source": [ "Alternatively, you can use depression breaching to fill the depressions." @@ -531,7 +573,7 @@ { "cell_type": "code", "execution_count": null, - "id": "53", + "id": "51", "metadata": {}, "outputs": [], "source": [ @@ -541,7 +583,7 @@ { "cell_type": "code", "execution_count": null, - "id": "54", + "id": "52", "metadata": {}, "outputs": [], "source": [ @@ -551,26 +593,29 @@ { "cell_type": "code", "execution_count": null, - "id": "55", + "id": "53", "metadata": {}, "outputs": [], "source": [ + "m.layers[-1].visible = False\n", "m.add_raster(\"noflow2.tif\", layer_name=\"No Flow Cells after Breaching\")\n", "m" ] }, { "cell_type": "markdown", - "id": "56", + "id": "54", "metadata": {}, "source": [ - "## Delineate flow direction" + "### Delineate flow direction\n", + "\n", + "Next, we delineate the flow direction based on the D8 algorithm." ] }, { "cell_type": "code", "execution_count": null, - "id": "57", + "id": "55", "metadata": {}, "outputs": [], "source": [ @@ -579,16 +624,18 @@ }, { "cell_type": "markdown", - "id": "58", + "id": "56", "metadata": {}, "source": [ - "## Calculate flow accumulation" + "### Calculate flow accumulation\n", + "\n", + "Now, calculate **flow accumulation** to understand how water collects across the landscape." ] }, { "cell_type": "code", "execution_count": null, - "id": "59", + "id": "57", "metadata": {}, "outputs": [], "source": [ @@ -598,26 +645,27 @@ { "cell_type": "code", "execution_count": null, - "id": "60", + "id": "58", "metadata": {}, "outputs": [], "source": [ - "m.add_raster(\"flow_accum.tif\", layer_name=\"Flow Accumulation\")\n", - "m" + "m.add_raster(\"flow_accum.tif\", layer_name=\"Flow Accumulation\")" ] }, { "cell_type": "markdown", - "id": "61", + "id": "59", "metadata": {}, "source": [ - "## Extract streams" + "### Extract streams\n", + "\n", + "Extract the stream network using the flow accumulation data." ] }, { "cell_type": "code", "execution_count": null, - "id": "62", + "id": "60", "metadata": {}, "outputs": [], "source": [ @@ -627,25 +675,26 @@ { "cell_type": "code", "execution_count": null, - "id": "63", + "id": "61", "metadata": {}, "outputs": [], "source": [ + "m.layers[-1].visible = False\n", "m.add_raster(\"streams.tif\", layer_name=\"Streams\")" ] }, { "cell_type": "markdown", - "id": "64", + "id": "62", "metadata": {}, "source": [ - "## Calculate distance to outlet" + "### Calculate distance to outlet" ] }, { "cell_type": "code", "execution_count": null, - "id": "65", + "id": "63", "metadata": {}, "outputs": [], "source": [ @@ -657,7 +706,7 @@ { "cell_type": "code", "execution_count": null, - "id": "66", + "id": "64", "metadata": {}, "outputs": [], "source": [ @@ -666,16 +715,16 @@ }, { "cell_type": "markdown", - "id": "67", + "id": "65", "metadata": {}, "source": [ - "## Vectorize streams" + "### Vectorize streams" ] }, { "cell_type": "code", "execution_count": null, - "id": "68", + "id": "66", "metadata": {}, "outputs": [], "source": [ @@ -686,16 +735,16 @@ }, { "cell_type": "markdown", - "id": "69", + "id": "67", "metadata": {}, "source": [ - "The raster_streams_to_vector tool has a bug. The output vector file is missing the coordinate system. Use leafmap.vector_set_crs() to set the coordinate system." + "The raster_streams_to_vector tool has a bug. The output vector file is missing the coordinate system. Use `leafmap.vector_set_crs()` to set the coordinate system." ] }, { "cell_type": "code", "execution_count": null, - "id": "70", + "id": "68", "metadata": {}, "outputs": [], "source": [ @@ -705,56 +754,51 @@ { "cell_type": "code", "execution_count": null, - "id": "71", + "id": "69", "metadata": {}, "outputs": [], "source": [ "m.add_shp(\n", - " \"streams.shp\", layer_name=\"Streams Vector\", style={\"color\": \"#ff0000\", \"weight\": 3}\n", + " \"streams.shp\",\n", + " layer_name=\"Streams Vector\",\n", + " style={\"color\": \"#ff0000\", \"weight\": 3},\n", + " info_mode=None,\n", ")\n", "m" ] }, { "cell_type": "markdown", - "id": "72", + "id": "70", "metadata": {}, "source": [ - "## Delineate basins" + "You can turn on the USGS Hydrography basemap to visualize the stream network and compare it with the extracted stream network." ] }, { - "cell_type": "code", - "execution_count": null, - "id": "73", + "cell_type": "markdown", + "id": "71", "metadata": {}, - "outputs": [], "source": [ - "wbt.basins(\"flow_direction.tif\", \"basins.tif\")" + "### Delineate the longest flow path\n", + "\n", + "You can delineate the longest flow path in the watershed." ] }, { "cell_type": "code", "execution_count": null, - "id": "74", + "id": "72", "metadata": {}, "outputs": [], "source": [ - "m.add_raster(\"basins.tif\", layer_name=\"Basins\")" - ] - }, - { - "cell_type": "markdown", - "id": "75", - "metadata": {}, - "source": [ - "## Delineate the longest flow path" + "wbt.basins(\"flow_direction.tif\", \"basins.tif\")" ] }, { "cell_type": "code", "execution_count": null, - "id": "76", + "id": "73", "metadata": {}, "outputs": [], "source": [ @@ -765,7 +809,7 @@ }, { "cell_type": "markdown", - "id": "77", + "id": "74", "metadata": {}, "source": [ "Select only the longest flow path." @@ -774,7 +818,7 @@ { "cell_type": "code", "execution_count": null, - "id": "78", + "id": "75", "metadata": {}, "outputs": [], "source": [ @@ -786,7 +830,7 @@ { "cell_type": "code", "execution_count": null, - "id": "79", + "id": "76", "metadata": {}, "outputs": [], "source": [ @@ -800,38 +844,44 @@ }, { "cell_type": "markdown", - "id": "80", + "id": "77", "metadata": {}, "source": [ - "## Generate a pour point" + "### Generate a pour point\n", + "\n", + "To delineate a watershed, you need to specify a pour point. You can use the outlet of the longest flow path as the pour point or specify a different point. Use the drawing tool to place a marker on the map to specify the pour point. If no marker is placed, the default pour point below will be used." ] }, { "cell_type": "code", "execution_count": null, - "id": "81", + "id": "78", "metadata": {}, "outputs": [], "source": [ "if m.user_roi is not None:\n", " m.save_draw_features(\"pour_point.shp\", crs=\"EPSG:3857\")\n", "else:\n", - " coords = [-122.613559, 44.284383]\n", - " leafmap.coords_to_vector(coords, output=\"pour_point.shp\", crs=\"EPSG:3857\")" + " lat = 44.284642\n", + " lon = -122.611217\n", + " leafmap.coords_to_vector([lon, lat], output=\"pour_point.shp\", crs=\"EPSG:3857\")\n", + " m.add_marker([lat, lon])" ] }, { "cell_type": "markdown", - "id": "82", + "id": "79", "metadata": {}, "source": [ - "## Snap pour point to stream" + "### Snap pour point to stream\n", + "\n", + "Snap the pour point to the nearest stream." ] }, { "cell_type": "code", "execution_count": null, - "id": "83", + "id": "80", "metadata": {}, "outputs": [], "source": [ @@ -840,38 +890,56 @@ ")" ] }, + { + "cell_type": "markdown", + "id": "81", + "metadata": {}, + "source": [ + "Visualize the snapped pour point on the map." + ] + }, { "cell_type": "code", "execution_count": null, - "id": "84", + "id": "82", "metadata": {}, "outputs": [], "source": [ - "m.add_shp(\"pour_point_snapped.shp\", layer_name=\"Pour Point\")" + "m.add_shp(\"pour_point_snapped.shp\", layer_name=\"Pour Point\", info_mode=False)" ] }, { "cell_type": "markdown", - "id": "85", + "id": "83", "metadata": {}, "source": [ - "## Delineate watershed" + "### Delineate watershed\n", + "\n", + "Delineate the watershed using a **pour point** and the flow direction data." ] }, { "cell_type": "code", "execution_count": null, - "id": "86", + "id": "84", "metadata": {}, "outputs": [], "source": [ "wbt.watershed(\"flow_direction.tif\", \"pour_point_snapped.shp\", \"watershed.tif\")" ] }, + { + "cell_type": "markdown", + "id": "85", + "metadata": {}, + "source": [ + "Visualize the watershed boundary on the map." + ] + }, { "cell_type": "code", "execution_count": null, - "id": "87", + "id": "86", "metadata": {}, "outputs": [], "source": [ @@ -881,22 +949,32 @@ }, { "cell_type": "markdown", - "id": "88", + "id": "87", "metadata": {}, "source": [ - "## Convert watershed raster to vector" + "### Convert watershed raster to vector\n", + "\n", + "You can convert the watershed raster to a vector file." ] }, { "cell_type": "code", "execution_count": null, - "id": "89", + "id": "88", "metadata": {}, "outputs": [], "source": [ "wbt.raster_to_vector_polygons(\"watershed.tif\", \"watershed.shp\")" ] }, + { + "cell_type": "markdown", + "id": "89", + "metadata": {}, + "source": [ + "Visualize the watershed boundary on the map." + ] + }, { "cell_type": "code", "execution_count": null, @@ -904,17 +982,403 @@ "metadata": {}, "outputs": [], "source": [ + "m.layers[-1].visible = False\n", "m.add_shp(\n", " \"watershed.shp\",\n", " layer_name=\"Watershed Vector\",\n", " style={\"color\": \"#ffff00\", \"weight\": 3},\n", + " info_mode=False,\n", ")" ] + }, + { + "cell_type": "markdown", + "id": "91", + "metadata": {}, + "source": [ + "## Part 2: LiDAR Data Analysis\n", + "\n", + "In this section, we will process LiDAR data to derive **Digital Surface Models (DSMs)**, **Digital Elevation Models (DEMs)**, and **Canopy Height Models (CHMs)**. We will also remove outliers and improve the quality of the LiDAR data.\n", + "\n", + "### Set up whitebox\n", + "\n", + "First, we set up the WhiteboxTools and leafmap libraries." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92", + "metadata": {}, + "outputs": [], + "source": [ + "import leafmap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93", + "metadata": {}, + "outputs": [], + "source": [ + "wbt = leafmap.WhiteboxTools()\n", + "wbt.set_working_dir(os.getcwd())\n", + "wbt.verbose = False" + ] + }, + { + "cell_type": "markdown", + "id": "94", + "metadata": {}, + "source": [ + "### Download a sample dataset\n", + "\n", + "We download a sample LiDAR dataset for **Madison**." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "95", + "metadata": {}, + "outputs": [], + "source": [ + "url = \"https://github.com/opengeos/datasets/releases/download/lidar/madison.zip\"\n", + "filename = \"madison.las\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96", + "metadata": {}, + "outputs": [], + "source": [ + "leafmap.download_file(url, \"madison.zip\", quiet=True)" + ] + }, + { + "cell_type": "markdown", + "id": "97", + "metadata": {}, + "source": [ + "### Read LAS/LAZ data\n", + "\n", + "Load and inspect the LiDAR data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "98", + "metadata": {}, + "outputs": [], + "source": [ + "laz = leafmap.read_lidar(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "99", + "metadata": {}, + "outputs": [], + "source": [ + "laz" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "100", + "metadata": {}, + "outputs": [], + "source": [ + "str(laz.header.version)" + ] + }, + { + "cell_type": "markdown", + "id": "101", + "metadata": {}, + "source": [ + "### Upgrade file version\n", + "\n", + "Upgrade the LAS file version to 1.4." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "102", + "metadata": {}, + "outputs": [], + "source": [ + "las = leafmap.convert_lidar(laz, file_version=\"1.4\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "103", + "metadata": {}, + "outputs": [], + "source": [ + "str(las.header.version)" + ] + }, + { + "cell_type": "markdown", + "id": "104", + "metadata": {}, + "source": [ + "### Write LAS data\n", + "\n", + "Save the LAS data to a new file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "105", + "metadata": {}, + "outputs": [], + "source": [ + "leafmap.write_lidar(las, \"madison.las\")" + ] + }, + { + "cell_type": "markdown", + "id": "106", + "metadata": {}, + "source": [ + "### Histogram analysis\n", + "\n", + "Plot the histogram of the LiDAR data. The histogram shows the distribution of the LiDAR points based on their elevation values. The tool generates a histogram of the LiDAR data and saves it to an HTM file. You can open the HTM file in a web browser to view the histogram." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "107", + "metadata": {}, + "outputs": [], + "source": [ + "wbt.lidar_histogram(\"madison.las\", \"histogram.html\")" + ] + }, + { + "cell_type": "markdown", + "id": "108", + "metadata": {}, + "source": [ + "### Visualize LiDAR data\n", + "\n", + "Run the `view_lidar()` function to visualize the LiDAR data in 3D. Note that the `view_lidar()` function may not work in some environments, such as Google Colab." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "109", + "metadata": {}, + "outputs": [], + "source": [ + "leafmap.view_lidar(\"madison.las\")" + ] + }, + { + "cell_type": "markdown", + "id": "110", + "metadata": {}, + "source": [ + "### Remove outliers\n", + "\n", + "Remove outliers from the LiDAR dataset:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "111", + "metadata": {}, + "outputs": [], + "source": [ + "wbt.lidar_elevation_slice(\"madison.las\", \"madison_rm.las\", minz=0, maxz=450)" + ] + }, + { + "cell_type": "markdown", + "id": "112", + "metadata": {}, + "source": [ + "### Visualize LiDAR data after removing outliers\n", + "\n", + "We can visualize the LiDAR data after removing the outliers." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "113", + "metadata": {}, + "outputs": [], + "source": [ + "leafmap.view_lidar(\"madison_rm.las\", cmap=\"terrain\")" + ] + }, + { + "cell_type": "markdown", + "id": "114", + "metadata": {}, + "source": [ + "### Create DSM\n", + "\n", + "Using the LiDAR data to create a **Digital Surface Model (DSM)**." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "115", + "metadata": {}, + "outputs": [], + "source": [ + "wbt.lidar_digital_surface_model(\n", + " \"madison_rm.las\", \"dsm.tif\", resolution=1.0, minz=0, maxz=450\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "116", + "metadata": {}, + "source": [ + "The DSM generated by whitebox is missing the coordinate system. Use `leafmap.raster_set_crs()` to set the coordinate system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "117", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "outputs": [], + "source": [ + "leafmap.add_crs(\"dsm.tif\", epsg=2255)" + ] + }, + { + "cell_type": "markdown", + "id": "118", + "metadata": {}, + "source": [ + "### Visualize DSM\n", + "\n", + "Visualize the DSM on the map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "119", + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_basemap(\"Satellite\")\n", + "m.add_raster(\"dsm.tif\", colormap=\"terrain\", layer_name=\"DSM\")\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "120", + "metadata": {}, + "source": [ + "### Create DEM\n", + "\n", + "We can create a bare-earth DEM from the DSM. The tool is typically applied to LiDAR DEMs which frequently contain numerous off-terrain objects (OTOs) such as buildings, trees and other vegetation, cars, fences and other anthropogenic objects. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "121", + "metadata": {}, + "outputs": [], + "source": [ + "wbt.remove_off_terrain_objects(\"dsm.tif\", \"dem.tif\", filter=25, slope=15.0)" + ] + }, + { + "cell_type": "markdown", + "id": "122", + "metadata": {}, + "source": [ + "### Visualize DEM\n", + "\n", + "Visualize the bear-earth DEM on the map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "123", + "metadata": {}, + "outputs": [], + "source": [ + "m.add_raster(\"dem.tif\", palette=\"terrain\", layer_name=\"DEM\")\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "124", + "metadata": {}, + "source": [ + "### Create CHM\n", + "\n", + "\n", + "We can a **Canopy Height Model (CHM)** by subtracting the DEM from the DSM." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "125", + "metadata": {}, + "outputs": [], + "source": [ + "chm = wbt.subtract(\"dsm.tif\", \"dem.tif\", \"chm.tif\")" + ] + }, + { + "cell_type": "markdown", + "id": "126", + "metadata": {}, + "source": [ + "Visualize the CHM on the map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "127", + "metadata": {}, + "outputs": [], + "source": [ + "m.add_raster(\"chm.tif\", palette=\"gist_earth\", layer_name=\"CHM\")\n", + "m.add_layer_manager()\n", + "m" + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -928,7 +1392,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.8" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/book/geospatial/whitebox.md b/book/geospatial/whitebox.md index 47f4663..bee0796 100644 --- a/book/geospatial/whitebox.md +++ b/book/geospatial/whitebox.md @@ -4,25 +4,82 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.16.2 + jupytext_version: 1.16.4 kernelspec: - display_name: Python 3 + display_name: Python 3 (ipykernel) language: python name: python3 --- -# WhiteboxTools +# Whitebox [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/giswqs/geog-312/blob/main/book/geospatial/whitebox.ipynb) +++ +## Overview + +In this lecture, we will explore the use of [**WhiteboxTools**](https://github.com/jblindsay/whitebox-tools), a powerful open-source library for performing geospatial analysis. Specifically, we will focus on two key applications: **watershed analysis** and **LiDAR data analysis**. You will learn how to manipulate geospatial data using Python, conduct hydrological analysis, and derive digital elevation models (DEMs) and canopy height models (CHMs) from LiDAR data. + +This lecture is structured into two main sections: +1. **Watershed Analysis**: Using DEMs and hydrological tools to delineate watersheds, calculate flow accumulation, and extract stream networks. +2. **LiDAR Data Analysis**: Processing LiDAR point cloud data to derive DEMs, DSMs, and CHMs while removing outliers and improving data quality. + +By the end of this session, you will have gained hands-on experience with **WhiteboxTools** and **leafmap**, allowing you to perform a wide range of geospatial and hydrological analyses. + +## Learning Objectives + +By the end of this lecture, you will be able to: +- Install and configure **WhiteboxTools** and **leafmap** for geospatial analysis. +- Create interactive maps to visualize basemaps and geospatial datasets. +- Perform watershed analysis by delineating watersheds, flow directions, and stream networks. +- Manipulate and analyze Digital Elevation Models (DEMs) to conduct hydrological modeling. +- Process and analyze LiDAR data to generate **Digital Surface Models (DSMs)**, **Digital Elevation Models (DEMs)**, and **Canopy Height Models (CHMs)**. +- Integrate **WhiteboxTools** with Python workflows to automate geospatial analysis. + + +## Introduction + +Below is a brief introduction to Whitebox, a powerful open-source library for geospatial analysis. For more information, please refer to the whiteboxgeo website: https://www.whiteboxgeo.com. + +### What is Whitebox? + +Whitebox is geospatial data analysis software originally developed at the University of Guelph‘s Geomorphometry and Hydrogeomatics Research Group ([GHRG](https://jblindsay.github.io/ghrg/index.html)) directed by Dr. John Lindsay. Whitebox contains over 550 tools for processing many types of geospatial data. It has many great features such as its extensive use of parallel computing, it doesn’t need other libraries to be installed (e.g., GDAL), it can be used from scripting environments, and it easily plugs into other geographical information system (GIS) software. The Whitebox Toolset Extension provides even more power. + +### What can Whitebox do? + +Whitebox can be used to perform common GIS and remote sensing analysis tasks. Whitebox also contains advanced tooling for spatial hydrological analysis and LiDAR data processing. Whitebox is not a cartographic or spatial data visualization package; instead it’s meant to serve as an analytical back-end for other data visualization software, like QGIS and ArcGIS. + +### How is Whitebox different? + +Whitebox doesn’t compete with QGIS, ArcGIS/Pro, and ArcPy but rather it extends them. You can plug WhiteboxTools into QGIS and ArcGIS and it’ll provide hundreds of additional tools for analyzing all kinds of geospatial data. You can also call Whitebox functions from Python scripts using [Whitebox Workflows](https://www.whiteboxgeo.com/whitebox-workflows-for-python) (WbW). Combine WbW with ArcPy to more effectively automate your data analysis workflows and streamline your geoprocessing solutions. + +There are many tools in Whitebox that you won’t find elsewhere. You can think of Whitebox as a portable, cross-platform GIS analysis powerhouse, allowing you to extend your preferred GIS or to embed Whitebox capabilities into your automated scripted workflows. Oh, and it’s fast, really fast! + + +## Useful Resources for Whitebox + +- GitHub Repository: https://github.com/jblindsay/whitebox-tools +- WhiteboxGeo: https://www.whiteboxgeo.com +- User Manual: https://whiteboxgeo.com/manual/wbt_book/preface.html +- Whitebox Workflows for Python: https://www.whiteboxgeo.com/whitebox-workflows-for-python +- Whitebox Workflows for Python Pro: https://www.whiteboxgeo.com/whitebox-workflows-professional +- Python Frontend: https://github.com/opengeos/whitebox-python +- Jupyter Frontend: https://github.com/opengeos/whiteboxgui +- R Frontend: https://github.com/opengeos/whiteboxR +- ArcGIS Toolbox: https://github.com/opengeos/WhiteboxTools-ArcGIS +- QGIS Plugin: https://plugins.qgis.org/plugins/whitebox_workflows_for_qgis + ## Installation -Uncomment and run the following cell to install necessary packages for this notebook, including leafmap, geopandas, localtileserver, rio-cogeo, pynhd, py3dep. +To get started, we need to install the required packages, such as `leafmap`, and `whitebox`. Uncomment the following lines to install the packages. ```{code-cell} ipython3 -# %pip install "leafmap[raster]" geopandas +# %pip install "leafmap[raster]" "leafmap[lidar]" mapclassify +``` + +```{code-cell} ipython3 +# %pip install numpy==1.26.4 ``` ## Import libraries @@ -30,20 +87,20 @@ Uncomment and run the following cell to install necessary packages for this note ```{code-cell} ipython3 import os import leafmap +import numpy as np ``` -## Create interactive maps - -Specify the map center, zoom level, and height. +Some of the raster datasets generated by whitebox will be int32 type with a nodata value of -32768. To make it easier to visualize these datasets, we set the nodata value as an environment variable, which will be used by leafmap to set the nodata value when visualizing the raster datasets. ```{code-cell} ipython3 -m = leafmap.Map(center=[40, -100], zoom=4, height="600px") -m +os.environ["NODATA"] = "-32768" ``` -## Add basemaps +## Part 1: Watershed Analysis + +### Create Interactive Maps -Add OpenTopoMap, USGS 3DEP Elevation, and USGS Hydrography basemaps. +To perform watershed analysis, we first create an interactive map using **leafmap**. This step allows us to visualize different basemaps. ```{code-cell} ipython3 m = leafmap.Map() @@ -53,77 +110,62 @@ m.add_basemap("USGS Hydrography") m ``` -Add NLCD land cover map and legend. - -```{code-cell} ipython3 -m = leafmap.Map(center=[40, -100], zoom=4) -m.add_basemap("HYBRID") -m.add_basemap("NLCD 2019 CONUS Land Cover") -m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type") -m -``` +### Download Watershed Data -Add WMS layers. +Next, we download watershed data for the **Calapooia River basin** in Oregon. We'll use the latitude and longitude of a point in the basin to extract watershed boundary data. ```{code-cell} ipython3 -m = leafmap.Map(center=[40, -100], zoom=4) -m.add_basemap("Esri.WorldImagery") -url = "https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?" -m.add_wms_layer( - url, - layers="NLCD_2019_Land_Cover_L48", - name="NLCD 2019 CONUS Land Cover", - format="image/png", - transparent=True, -) -m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type") +lat = 44.361169 +lon = -122.821802 + +m = leafmap.Map(center=[lat, lon], zoom=10) +m.add_marker([lat, lon]) m ``` -## Get watershed data - -Let's download watershed data for the Calapooia River basin in Oregon. +Download the watershed data and visualize it: ```{code-cell} ipython3 -gdf = leafmap.get_nhd_basins(feature_ids=23763529, fsource="comid", simplified=False) +geometry = {"x": lon, "y": lat} ``` -Plot the watershed boundary on the map. - ```{code-cell} ipython3 -m = leafmap.Map() -m.add_gdf(gdf, layer_name="Catchment", info_mode=None) -m +gdf = leafmap.get_wbd(geometry, digit=10, return_geometry=True) +gdf.explore() ``` -Save the watershed boundary to a GeoJSON or shapefile. +Save the watershed boundary to a **GeoJSON** file: ```{code-cell} ipython3 -gdf.to_file("basin.geojson", driver="GeoJSON") +gdf.to_file("basin.geojson") ``` -```{code-cell} ipython3 -gdf.to_file("basin.shp") -``` - -## Download DEM +### Download and Display DEM -Download a digital elevation model (DEM) for the watershed from the USGS 3DEP Elevation service. Convert the DEM to a Cloud Optimized GeoTIFF (COG). +We download a **Digital Elevation Model (DEM)** from the USGS 3DEP Elevation service to analyze the terrain of the watershed. The DEM will be used to delineate watersheds, calculate flow accumulation, and extract stream networks. The `leafmap.get_3dep_dem()` function returns the DEM as an `xarray.DataArray` object. Optionally, you can save the DEM to a GeoTIFF file by setting the `output` parameter. ```{code-cell} ipython3 -leafmap.get_3dep_dem( - gdf, resolution=30, output="dem.tif", dst_crs="EPSG:3857", to_cog=True +array = leafmap.get_3dep_dem( + gdf, + resolution=30, + output="dem.tif", + dst_crs="EPSG:3857", + to_cog=True, + overwrite=True, ) +array ``` -Display the DEM on the map. +Visualize the DEM on the map: ```{code-cell} ipython3 -m.add_raster("dem.tif", palette="terrain", layer_name="DEM") +m.add_raster("dem.tif", palette="terrain", nodata=np.nan, layer_name="DEM") m ``` -## Get DEM metadata +### Get DEM metadata + +We can get the metadata of the DEM, such as the spatial resolution, bounding box, and coordinate reference system (CRS). ```{code-cell} ipython3 metadata = leafmap.image_metadata("dem.tif") @@ -132,18 +174,21 @@ metadata Get a summary statistics of the DEM. ++++ + +### Add colorbar + +Add a colorbar to the map to show the elevation values. Use the `image_min_max()` function to get the minimum and maximum values of the DEM. + ```{code-cell} ipython3 -metadata["bands"] +leafmap.image_min_max("dem.tif") ``` -## Add colorbar - ```{code-cell} ipython3 m.add_colormap(cmap="terrain", vmin="60", vmax=1500, label="Elevation (m)") -m ``` -## Initialize WhiteboxTools +### Initialize WhiteboxTools Initialize the WhiteboxTools class. @@ -157,37 +202,49 @@ Check the WhiteboxTools version. wbt.version() ``` -Display the WhiteboxTools interface. +Display the WhiteboxTools interface, which lists all available tools. You can use this interface to search for specific tools. You can also run any tool using the interface. However, we will use the Python API to run the tools. ```{code-cell} ipython3 leafmap.whiteboxgui() ``` -## Set working directory +### Set working directory + +Set the working directory to save the intermediate and output files. Set `wbt.version=False` to suppress the Whitebox processing log. ```{code-cell} ipython3 wbt.set_working_dir(os.getcwd()) wbt.verbose = False ``` -## Smooth DEM +### Smooth DEM All WhiteboxTools functions will return 0 if they are successful, and 1 if they are not. +Smoothing the DEM enhances the quality of subsequent hydrological analysis. + ```{code-cell} ipython3 wbt.feature_preserving_smoothing("dem.tif", "smoothed.tif", filter=9) ``` -Display the smoothed DEM and watershed boundary on the map. +Visualize the smoothed DEM and watershed boundary on the map. + ++++ + +Visualize the smoothed DEM: ```{code-cell} ipython3 m = leafmap.Map() -m.add_raster("smoothed.tif", palette="terrain", layer_name="Smoothed DEM") +m.add_basemap("Satellite") +m.add_raster("smoothed.tif", colormap="terrain", layer_name="Smoothed DEM") m.add_geojson("basin.geojson", layer_name="Watershed", info_mode=None) +m.add_basemap("USGS Hydrography", show=False) m ``` -## Create hillshade +### Create hillshade + +Create a hillshade from the smoothed DEM. ```{code-cell} ipython3 wbt.hillshade("smoothed.tif", "hillshade.tif", azimuth=315, altitude=35) @@ -200,22 +257,23 @@ m.add_raster("hillshade.tif", layer_name="Hillshade") m.layers[-1].opacity = 0.6 ``` -## Find no-flow cells +### Find no-flow cells -Find cells with undefined flow, i.e. no valid flow direction, based on the D8 flow direction algorithm +Find cells with undefined flow, i.e. no valid flow direction, based on the D8 flow direction algorithm. ```{code-cell} ipython3 wbt.find_no_flow_cells("smoothed.tif", "noflow.tif") ``` -Display the no-flow cells on the map. +Visualize the no-flow cells on the map. ```{code-cell} ipython3 m.add_raster("noflow.tif", layer_name="No Flow Cells") -m ``` -## Fill depressions +### Fill depressions + +First, we fill any depressions in the DEM to ensure proper flow simulation. ```{code-cell} ipython3 wbt.fill_depressions("smoothed.tif", "filled.tif") @@ -232,17 +290,22 @@ wbt.find_no_flow_cells("breached.tif", "noflow2.tif") ``` ```{code-cell} ipython3 +m.layers[-1].visible = False m.add_raster("noflow2.tif", layer_name="No Flow Cells after Breaching") m ``` -## Delineate flow direction +### Delineate flow direction + +Next, we delineate the flow direction based on the D8 algorithm. ```{code-cell} ipython3 wbt.d8_pointer("breached.tif", "flow_direction.tif") ``` -## Calculate flow accumulation +### Calculate flow accumulation + +Now, calculate **flow accumulation** to understand how water collects across the landscape. ```{code-cell} ipython3 wbt.d8_flow_accumulation("breached.tif", "flow_accum.tif") @@ -250,20 +313,22 @@ wbt.d8_flow_accumulation("breached.tif", "flow_accum.tif") ```{code-cell} ipython3 m.add_raster("flow_accum.tif", layer_name="Flow Accumulation") -m ``` -## Extract streams +### Extract streams + +Extract the stream network using the flow accumulation data. ```{code-cell} ipython3 wbt.extract_streams("flow_accum.tif", "streams.tif", threshold=5000) ``` ```{code-cell} ipython3 +m.layers[-1].visible = False m.add_raster("streams.tif", layer_name="Streams") ``` -## Calculate distance to outlet +### Calculate distance to outlet ```{code-cell} ipython3 wbt.distance_to_outlet( @@ -275,7 +340,7 @@ wbt.distance_to_outlet( m.add_raster("distance_to_outlet.tif", layer_name="Distance to Outlet") ``` -## Vectorize streams +### Vectorize streams ```{code-cell} ipython3 wbt.raster_streams_to_vector( @@ -283,7 +348,7 @@ wbt.raster_streams_to_vector( ) ``` -The raster_streams_to_vector tool has a bug. The output vector file is missing the coordinate system. Use leafmap.vector_set_crs() to set the coordinate system. +The raster_streams_to_vector tool has a bug. The output vector file is missing the coordinate system. Use `leafmap.vector_set_crs()` to set the coordinate system. ```{code-cell} ipython3 leafmap.vector_set_crs(source="streams.shp", output="streams.shp", crs="EPSG:3857") @@ -291,23 +356,26 @@ leafmap.vector_set_crs(source="streams.shp", output="streams.shp", crs="EPSG:385 ```{code-cell} ipython3 m.add_shp( - "streams.shp", layer_name="Streams Vector", style={"color": "#ff0000", "weight": 3} + "streams.shp", + layer_name="Streams Vector", + style={"color": "#ff0000", "weight": 3}, + info_mode=None, ) m ``` -## Delineate basins +You can turn on the USGS Hydrography basemap to visualize the stream network and compare it with the extracted stream network. -```{code-cell} ipython3 -wbt.basins("flow_direction.tif", "basins.tif") -``` ++++ + +### Delineate the longest flow path + +You can delineate the longest flow path in the watershed. ```{code-cell} ipython3 -m.add_raster("basins.tif", layer_name="Basins") +wbt.basins("flow_direction.tif", "basins.tif") ``` -## Delineate the longest flow path - ```{code-cell} ipython3 wbt.longest_flowpath( dem="breached.tif", basins="basins.tif", output="longest_flowpath.shp" @@ -331,17 +399,23 @@ m.add_shp( m ``` -## Generate a pour point +### Generate a pour point + +To delineate a watershed, you need to specify a pour point. You can use the outlet of the longest flow path as the pour point or specify a different point. Use the drawing tool to place a marker on the map to specify the pour point. If no marker is placed, the default pour point below will be used. ```{code-cell} ipython3 if m.user_roi is not None: m.save_draw_features("pour_point.shp", crs="EPSG:3857") else: - coords = [-122.613559, 44.284383] - leafmap.coords_to_vector(coords, output="pour_point.shp", crs="EPSG:3857") + lat = 44.284642 + lon = -122.611217 + leafmap.coords_to_vector([lon, lat], output="pour_point.shp", crs="EPSG:3857") + m.add_marker([lat, lon]) ``` -## Snap pour point to stream +### Snap pour point to stream + +Snap the pour point to the nearest stream. ```{code-cell} ipython3 wbt.snap_pour_points( @@ -349,31 +423,205 @@ wbt.snap_pour_points( ) ``` +Visualize the snapped pour point on the map. + ```{code-cell} ipython3 -m.add_shp("pour_point_snapped.shp", layer_name="Pour Point") +m.add_shp("pour_point_snapped.shp", layer_name="Pour Point", info_mode=False) ``` -## Delineate watershed +### Delineate watershed + +Delineate the watershed using a **pour point** and the flow direction data. ```{code-cell} ipython3 wbt.watershed("flow_direction.tif", "pour_point_snapped.shp", "watershed.tif") ``` +Visualize the watershed boundary on the map. + ```{code-cell} ipython3 m.add_raster("watershed.tif", layer_name="Watershed") m ``` -## Convert watershed raster to vector +### Convert watershed raster to vector + +You can convert the watershed raster to a vector file. ```{code-cell} ipython3 wbt.raster_to_vector_polygons("watershed.tif", "watershed.shp") ``` +Visualize the watershed boundary on the map. + ```{code-cell} ipython3 +m.layers[-1].visible = False m.add_shp( "watershed.shp", layer_name="Watershed Vector", style={"color": "#ffff00", "weight": 3}, + info_mode=False, ) ``` + +## Part 2: LiDAR Data Analysis + +In this section, we will process LiDAR data to derive **Digital Surface Models (DSMs)**, **Digital Elevation Models (DEMs)**, and **Canopy Height Models (CHMs)**. We will also remove outliers and improve the quality of the LiDAR data. + +### Set up whitebox + +First, we set up the WhiteboxTools and leafmap libraries. + +```{code-cell} ipython3 +import leafmap +``` + +```{code-cell} ipython3 +wbt = leafmap.WhiteboxTools() +wbt.set_working_dir(os.getcwd()) +wbt.verbose = False +``` + +### Download a sample dataset + +We download a sample LiDAR dataset for **Madison**. + +```{code-cell} ipython3 +url = "https://github.com/opengeos/datasets/releases/download/lidar/madison.zip" +filename = "madison.las" +``` + +```{code-cell} ipython3 +leafmap.download_file(url, "madison.zip", quiet=True) +``` + +### Read LAS/LAZ data + +Load and inspect the LiDAR data: + +```{code-cell} ipython3 +laz = leafmap.read_lidar(filename) +``` + +```{code-cell} ipython3 +laz +``` + +```{code-cell} ipython3 +str(laz.header.version) +``` + +### Upgrade file version + +Upgrade the LAS file version to 1.4. + +```{code-cell} ipython3 +las = leafmap.convert_lidar(laz, file_version="1.4") +``` + +```{code-cell} ipython3 +str(las.header.version) +``` + +### Write LAS data + +Save the LAS data to a new file. + +```{code-cell} ipython3 +leafmap.write_lidar(las, "madison.las") +``` + +### Histogram analysis + +Plot the histogram of the LiDAR data. The histogram shows the distribution of the LiDAR points based on their elevation values. The tool generates a histogram of the LiDAR data and saves it to an HTM file. You can open the HTM file in a web browser to view the histogram. + +```{code-cell} ipython3 +wbt.lidar_histogram("madison.las", "histogram.html") +``` + +### Visualize LiDAR data + +Run the `view_lidar()` function to visualize the LiDAR data in 3D. Note that the `view_lidar()` function may not work in some environments, such as Google Colab. + +```{code-cell} ipython3 +leafmap.view_lidar("madison.las") +``` + +### Remove outliers + +Remove outliers from the LiDAR dataset: + +```{code-cell} ipython3 +wbt.lidar_elevation_slice("madison.las", "madison_rm.las", minz=0, maxz=450) +``` + +### Visualize LiDAR data after removing outliers + +We can visualize the LiDAR data after removing the outliers. + +```{code-cell} ipython3 +leafmap.view_lidar("madison_rm.las", cmap="terrain") +``` + +### Create DSM + +Using the LiDAR data to create a **Digital Surface Model (DSM)**. + +```{code-cell} ipython3 +wbt.lidar_digital_surface_model( + "madison_rm.las", "dsm.tif", resolution=1.0, minz=0, maxz=450 +) +``` + +The DSM generated by whitebox is missing the coordinate system. Use `leafmap.raster_set_crs()` to set the coordinate system. + +```{code-cell} ipython3 +:jp-MarkdownHeadingCollapsed: true + +leafmap.add_crs("dsm.tif", epsg=2255) +``` + +### Visualize DSM + +Visualize the DSM on the map. + +```{code-cell} ipython3 +m = leafmap.Map() +m.add_basemap("Satellite") +m.add_raster("dsm.tif", colormap="terrain", layer_name="DSM") +m +``` + +### Create DEM + +We can create a bare-earth DEM from the DSM. The tool is typically applied to LiDAR DEMs which frequently contain numerous off-terrain objects (OTOs) such as buildings, trees and other vegetation, cars, fences and other anthropogenic objects. + +```{code-cell} ipython3 +wbt.remove_off_terrain_objects("dsm.tif", "dem.tif", filter=25, slope=15.0) +``` + +### Visualize DEM + +Visualize the bear-earth DEM on the map. + +```{code-cell} ipython3 +m.add_raster("dem.tif", palette="terrain", layer_name="DEM") +m +``` + +### Create CHM + + +We can a **Canopy Height Model (CHM)** by subtracting the DEM from the DSM. + +```{code-cell} ipython3 +chm = wbt.subtract("dsm.tif", "dem.tif", "chm.tif") +``` + +Visualize the CHM on the map. + +```{code-cell} ipython3 +m.add_raster("chm.tif", palette="gist_earth", layer_name="CHM") +m.add_layer_manager() +m +``` diff --git a/requirements.txt b/requirements.txt index 7c1c837..4fe3823 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,4 @@ ghp-import==2.1.* jupyter-book==1.0.* jupytext==1.16.* pmtiles==3.4.* +py3dep