diff --git a/MANIFEST.in b/MANIFEST.in index 7a1e512..bbfadbe 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,2 +1 @@ include LICENSE README.md HISTORY.rst AUTHORS.rst -include dem_30m.png tpi_500M.png \ No newline at end of file diff --git a/README.ipynb b/README.ipynb new file mode 100644 index 0000000..203b675 --- /dev/null +++ b/README.ipynb @@ -0,0 +1,231 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# topo-descriptors\n", + "\n", + "A python library to compute DEM-based topographical descriptors.\n", + "\n", + "## Usage\n", + "\n", + "Let's install `topo-descriptors` with few additional packages that will help us\n", + "to run a simple example (remember to use a virtual environment):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install topo-descriptors elevation rioxarray matplotlib --quiet" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The [elevation](https://github.com/bopen/elevation) package is an python library that\n", + "provides an easy access to global elevation data. Here we are going to clip the SRTM 30m\n", + "DEM around the Basodino region in southern Switzerland, around 46.4N 8.5E:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "!eio clip -o Basodino-30m-DEM.tif --bounds 8.2 46.30 8.6 46.55" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import logging\n", + "\n", + "logger = logging.getLogger()\n", + "handler = logging.StreamHandler()\n", + "formatter = logging.Formatter(\"%(asctime)s %(name)-12s %(levelname)-8s %(message)s\")\n", + "handler.setFormatter(formatter)\n", + "logger.addHandler(handler)\n", + "logger.setLevel(logging.INFO)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now in python we can use the xarray interface to rasterio to easily import the\n", + "`Basodino-30m-DEM.tif` file generated above:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "\n", + "dem = xr.open_rasterio(\"Basodino-30m-DEM.tif\")\n", + "dem = dem.isel(band=0, drop=True)\n", + "dem.plot(robust=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from topo_descriptors import topo, helpers\n", + "\n", + "scale_meters = 500\n", + "scale_pixel, __ = helpers.scale_to_pixel(scale_meters, dem)\n", + "topo.tpi(dem, scale_pixel).plot(vmin=-100, vmax=100, cmap=\"bwr\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Sx is used to describe the horizon in a given direction and spatial scale.\n", + "In the example below we compute the Sx for a 0° azimuth (i.e., looking North)\n", + "and a radius of 500 meters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sx_500m = topo.sx(dem, azimuth=0, radius=500)\n", + "xr.DataArray(sx_500m, coords=dem.coords).plot.imshow()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Other topographical descriptors are available, such as slope, aspect, derivatives,\n", + "and more. As an example, below we show how to compute a range of descriptors for two\n", + "distinc spatial scales (200 and 2000 meters)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "output_dir = Path(\"out/\")\n", + "output_dir.mkdir(exist_ok=True)\n", + "\n", + "scales_meters = [200, 2000]\n", + "domain = {\"x\": slice(8.25, 8.55), \"y\": slice(46.50, 46.35)}\n", + "\n", + "topo.compute_gradient(dem, scales_meters, sig_ratios=1, crop=domain, outdir=output_dir)\n", + "topo.compute_std(dem, scales_meters, crop=domain, outdir=output_dir)\n", + "topo.compute_tpi(dem, scales_meters, crop=domain, outdir=output_dir)\n", + "topo.compute_sx(dem, azimuth=0, radius=scales_meters[0], crop=domain, outdir=output_dir)\n", + "topo.compute_sx(dem, azimuth=0, radius=scales_meters[1], crop=domain, outdir=output_dir)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Above, the output was written directly to disk, while in the cell below we show how \n", + "to easly import the results and visualize them using xarray." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_mfdataset(str(output_dir / \"topo_*.nc\"))\n", + "min_max = ds.quantile(q=[0.05, 0.95])\n", + "ds = (ds - min_max.isel(quantile=0)) / (\n", + " min_max.isel(quantile=1) - min_max.isel(quantile=0)\n", + ")\n", + "ds.to_array().plot.imshow(\n", + " col=\"variable\",\n", + " col_wrap=len(scales_meters),\n", + " robust=True,\n", + " add_colorbar=False,\n", + " vmin=0,\n", + " vmax=1,\n", + ")\n", + "ds.close()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Build the README\n", + "\n", + "To use this Jupyter Notebook to compile the markdown's version for GitHub, first install\n", + "the conda environment using the `environment.yaml` file:\n", + "\n", + "```shell\n", + "conda env create -f environment.yaml\n", + "conda activate topo\n", + "```\n", + "\n", + "Then generate the `README.md` by running:\n", + "\n", + "```shell\n", + "jupyter nbconvert --execute --to markdown README.ipynb\n", + "```\n", + "\n", + "The associated figures are saved in the `README_files/` folder." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.10.8 ('topo')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:25:29) [Clang 14.0.6 ]" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "deb7fae88368ab96eb09766f06009a343f9a7d51b6b1f227ac5998e9f38709ff" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/README.md b/README.md index 845179a..a84082b 100644 --- a/README.md +++ b/README.md @@ -2,42 +2,76 @@ A python library to compute DEM-based topographical descriptors. - ## Usage +## Usage Let's install `topo-descriptors` with few additional packages that will help us -to run a simple example: +to run a simple example (remember to use a virtual environment): + +```python +%pip install topo-descriptors elevation rioxarray matplotlib --quiet ``` -pip install topo-descriptors -pip install elevation -pip install rasterio -``` -[elevation](https://github.com/bopen/elevation) is an python library that provides -an easy access to global elevation data. Here we are going to clip the SRTM 30m + Note: you may need to restart the kernel to use updated packages. + + +The [elevation](https://github.com/bopen/elevation) package is an python library that +provides an easy access to global elevation data. Here we are going to clip the SRTM 30m DEM around the Basodino region in southern Switzerland, around 46.4N 8.5E: + +```python +!eio clip -o Basodino-30m-DEM.tif --bounds 8.2 46.30 8.6 46.55 ``` -eio clip -o Basodino-30m-DEM.tif --bounds 8.2 46.30 8.6 46.55 + + make: Nothing to be done for `download'. + make: Nothing to be done for `all'. + cp SRTM1.vrt SRTM1.0e5622d0845a4ad9a0f3e44ad90db19d.vrt + gdal_translate -q -co TILED=YES -co COMPRESS=DEFLATE -co ZLEVEL=9 -co PREDICTOR=2 -projwin 8.2 46.55 8.6 46.3 SRTM1.0e5622d0845a4ad9a0f3e44ad90db19d.vrt /Users/daniele/src/topo-descriptors/Basodino-30m-DEM.tif + rm -f SRTM1.0e5622d0845a4ad9a0f3e44ad90db19d.vrt + + + +```python +import logging + +logger = logging.getLogger() +handler = logging.StreamHandler() +formatter = logging.Formatter("%(asctime)s %(name)-12s %(levelname)-8s %(message)s") +handler.setFormatter(formatter) +logger.addHandler(handler) +logger.setLevel(logging.INFO) ``` -Now in python we can use the xarray interface to rasterio to easily import our -`Basodino-30m-DEM.tif` file: +Now in python we can use the xarray interface to rasterio to easily import the +`Basodino-30m-DEM.tif` file generated above: + ```python import xarray as xr dem = xr.open_rasterio("Basodino-30m-DEM.tif") dem = dem.isel(band=0, drop=True) -dem.plot() +dem.plot(robust=True) ``` -![png](dem_30m.png) + + /var/folders/v1/d7jjg8d52fl77y27qbv75bw80000gp/T/ipykernel_9632/102662558.py:3: DeprecationWarning: open_rasterio is Deprecated in favor of rioxarray. For information about transitioning, see: https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html + dem = xr.open_rasterio("Basodino-30m-DEM.tif") + + + + + + + + + + + +![png](README_files/README_6_2.png) + -With the DEM data imported as a xarray DataArray, we can use topo-descriptors to -compute two established topographical descriptors: the topographical position index -(TPI) and the Sx. The TPI describes the elevation of a given point relatively -to its neighbourhood, here defined on a 500 meter scale: ```python from topo_descriptors import topo, helpers @@ -46,17 +80,149 @@ scale_meters = 500 scale_pixel, __ = helpers.scale_to_pixel(scale_meters, dem) topo.tpi(dem, scale_pixel).plot(vmin=-100, vmax=100, cmap="bwr") ``` -![png](tpi_500M.png) + + 2023-01-09 19:11:08,498 yaconfigobject INFO Loading /Users/daniele/src/topo-descriptors/topo_descriptors/config/topo_descriptors.conf. + 2023-01-09 19:11:08,498 yaconfigobject INFO Loading configuration file: /Users/daniele/src/topo-descriptors/topo_descriptors/config/topo_descriptors.conf + 2023-01-09 19:11:09,032 topo_descriptors.helpers INFO Computed in 0:00:00 (HH:mm:ss) + + + + + + + + + + + +![png](README_files/README_7_2.png) + + The Sx is used to describe the horizon in a given direction and spatial scale. In the example below we compute the Sx for a 0° azimuth (i.e., looking North) and a radius of 500 meters. + ```python sx_500m = topo.sx(dem, azimuth=0, radius=500) xr.DataArray(sx_500m, coords=dem.coords).plot.imshow() ``` -![png](sx_0azimuth_500radius.png) + + OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead. + 2023-01-09 19:11:13,162 topo_descriptors.helpers INFO Computed in 0:00:03 (HH:mm:ss) + + + + + + + + + + + +![png](README_files/README_9_2.png) + + Other topographical descriptors are available, such as slope, aspect, derivatives, -and more. \ No newline at end of file +and more. As an example, below we show how to compute a range of descriptors for two +distinc spatial scales (200 and 2000 meters). + + +```python +from pathlib import Path + +output_dir = Path("out/") +output_dir.mkdir(exist_ok=True) + +scales_meters = [200, 2000] +domain = {"x": slice(8.25, 8.55), "y": slice(46.50, 46.35)} + +topo.compute_gradient(dem, scales_meters, sig_ratios=1, crop=domain, outdir=output_dir) +topo.compute_std(dem, scales_meters, crop=domain, outdir=output_dir) +topo.compute_tpi(dem, scales_meters, crop=domain, outdir=output_dir) +topo.compute_sx(dem, azimuth=0, radius=scales_meters[0], crop=domain, outdir=output_dir) +topo.compute_sx(dem, azimuth=0, radius=scales_meters[1], crop=domain, outdir=output_dir) +``` + + 2023-01-09 19:11:13,301 topo_descriptors.topo INFO ***Starting gradients computation for scales [200, 2000] meters*** + 2023-01-09 19:11:13,432 topo_descriptors.topo INFO Computing scale 200 meters with sigma ratio 1 ... + 2023-01-09 19:11:13,476 topo_descriptors.helpers INFO Computed in 0:00:00 (HH:mm:ss) + 2023-01-09 19:11:13,487 topo_descriptors.helpers INFO saved: out/topo_WE_DERIVATIVE_200M_SIGRATIO1.nc + 2023-01-09 19:11:13,491 topo_descriptors.helpers INFO saved: out/topo_SN_DERIVATIVE_200M_SIGRATIO1.nc + 2023-01-09 19:11:13,495 topo_descriptors.helpers INFO saved: out/topo_SLOPE_200M_SIGRATIO1.nc + 2023-01-09 19:11:13,499 topo_descriptors.helpers INFO saved: out/topo_ASPECT_200M_SIGRATIO1.nc + 2023-01-09 19:11:13,499 topo_descriptors.topo INFO Computing scale 2000 meters with sigma ratio 1 ... + 2023-01-09 19:11:13,657 topo_descriptors.helpers INFO Computed in 0:00:00 (HH:mm:ss) + 2023-01-09 19:11:13,662 topo_descriptors.helpers INFO saved: out/topo_WE_DERIVATIVE_2000M_SIGRATIO1.nc + 2023-01-09 19:11:13,665 topo_descriptors.helpers INFO saved: out/topo_SN_DERIVATIVE_2000M_SIGRATIO1.nc + 2023-01-09 19:11:13,669 topo_descriptors.helpers INFO saved: out/topo_SLOPE_2000M_SIGRATIO1.nc + 2023-01-09 19:11:13,673 topo_descriptors.helpers INFO saved: out/topo_ASPECT_2000M_SIGRATIO1.nc + 2023-01-09 19:11:13,675 topo_descriptors.topo INFO ***Starting STD computation for scales [200, 2000] meters*** + 2023-01-09 19:11:13,808 topo_descriptors.topo INFO Computing scale 200 meters with smoothing factor None ... + 2023-01-09 19:11:13,860 topo_descriptors.helpers INFO Computed in 0:00:00 (HH:mm:ss) + 2023-01-09 19:11:13,865 topo_descriptors.helpers INFO saved: out/topo_STD_200M.nc + 2023-01-09 19:11:13,865 topo_descriptors.topo INFO Computing scale 2000 meters with smoothing factor None ... + 2023-01-09 19:11:13,923 topo_descriptors.helpers INFO Computed in 0:00:00 (HH:mm:ss) + 2023-01-09 19:11:13,928 topo_descriptors.helpers INFO saved: out/topo_STD_2000M.nc + 2023-01-09 19:11:13,929 topo_descriptors.topo INFO ***Starting TPI computation for scales [200, 2000] meters*** + 2023-01-09 19:11:14,061 topo_descriptors.topo INFO Computing scale 200 meters with smoothing factor None ... + 2023-01-09 19:11:14,086 topo_descriptors.helpers INFO Computed in 0:00:00 (HH:mm:ss) + 2023-01-09 19:11:14,090 topo_descriptors.helpers INFO saved: out/topo_TPI_200M.nc + 2023-01-09 19:11:14,090 topo_descriptors.topo INFO Computing scale 2000 meters with smoothing factor None ... + 2023-01-09 19:11:14,116 topo_descriptors.helpers INFO Computed in 0:00:00 (HH:mm:ss) + 2023-01-09 19:11:14,120 topo_descriptors.helpers INFO saved: out/topo_TPI_2000M.nc + 2023-01-09 19:11:14,121 topo_descriptors.topo INFO ***Starting Sx computation for azimuth 0 meters and radius 200*** + 2023-01-09 19:11:16,259 topo_descriptors.helpers INFO Computed in 0:00:02 (HH:mm:ss) + 2023-01-09 19:11:16,263 topo_descriptors.helpers INFO saved: out/topo_SX_RADIUS200_AZIMUTH0.nc + 2023-01-09 19:11:16,263 topo_descriptors.topo INFO ***Starting Sx computation for azimuth 0 meters and radius 2000*** + 2023-01-09 19:11:18,183 topo_descriptors.helpers INFO Computed in 0:00:01 (HH:mm:ss) + 2023-01-09 19:11:18,187 topo_descriptors.helpers INFO saved: out/topo_SX_RADIUS2000_AZIMUTH0.nc + + +Above, the output was written directly to disk, while in the cell below we show how +to easly import the results and visualize them using xarray. + + +```python +ds = xr.open_mfdataset(str(output_dir / "topo_*.nc")) +min_max = ds.quantile(q=[0.05, 0.95]) +ds = (ds - min_max.isel(quantile=0)) / ( + min_max.isel(quantile=1) - min_max.isel(quantile=0) +) +ds.to_array().plot.imshow( + col="variable", + col_wrap=len(scales_meters), + robust=True, + add_colorbar=False, + vmin=0, + vmax=1, +) +ds.close() +``` + + + +![png](README_files/README_13_0.png) + + + +## Build the README + +To use this Jupyter Notebook to compile the markdown's version for GitHub, first install +the conda environment using the `environment.yaml` file: + +```shell +conda env create -f environment.yaml +conda activate topo +``` + +Then generate the `README.md` by running: + +```shell +jupyter nbconvert --execute --to markdown README.ipynb +``` + +The associated figures are saved in the `README_files/` folder. diff --git a/README_files/README_13_0.png b/README_files/README_13_0.png new file mode 100644 index 0000000..ed0f570 Binary files /dev/null and b/README_files/README_13_0.png differ diff --git a/README_files/README_6_2.png b/README_files/README_6_2.png new file mode 100644 index 0000000..3880cdd Binary files /dev/null and b/README_files/README_6_2.png differ diff --git a/README_files/README_7_2.png b/README_files/README_7_2.png new file mode 100644 index 0000000..07d68c2 Binary files /dev/null and b/README_files/README_7_2.png differ diff --git a/README_files/README_9_2.png b/README_files/README_9_2.png new file mode 100644 index 0000000..909a07f Binary files /dev/null and b/README_files/README_9_2.png differ diff --git a/dem_30m.png b/dem_30m.png deleted file mode 100755 index f101ea4..0000000 Binary files a/dem_30m.png and /dev/null differ diff --git a/environment.yaml b/environment.yaml new file mode 100644 index 0000000..f7a67ac --- /dev/null +++ b/environment.yaml @@ -0,0 +1,8 @@ +name: topo +channels: +- conda-forge +- defaults +dependencies: + - python=3.10 + - gdal + - jupyter \ No newline at end of file diff --git a/scripts/compute_topo_descriptors.py b/scripts/compute_topo_descriptors.py index c135825..d512290 100755 --- a/scripts/compute_topo_descriptors.py +++ b/scripts/compute_topo_descriptors.py @@ -47,7 +47,9 @@ ) # TPI with prior smoothing - tp.compute_tpi(dem_da, scales_meters, smth_factors=1, ind_nans=ind_nans, crop=domain) + tp.compute_tpi( + dem_da, scales_meters, smth_factors=1, ind_nans=ind_nans, crop=domain + ) # Gradients with symmetric kernels tp.compute_gradient( diff --git a/sx_0azimuth_500radius.png b/sx_0azimuth_500radius.png deleted file mode 100644 index a66676d..0000000 Binary files a/sx_0azimuth_500radius.png and /dev/null differ diff --git a/test/test_helpers.py b/test/test_helpers.py index 47abfac..8f291ee 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -2,6 +2,7 @@ import topo_descriptors.helpers as hlp + def test_round_up_to_odd(): inputs = np.arange(0.1, 10, 0.7) outputs = hlp.round_up_to_odd(inputs) diff --git a/topo_descriptors/helpers.py b/topo_descriptors/helpers.py index edae040..f3905ea 100755 --- a/topo_descriptors/helpers.py +++ b/topo_descriptors/helpers.py @@ -63,6 +63,7 @@ def to_netcdf(array, coords, name, crop=None, outdir="."): da = xr.DataArray(array, coords=coords, name=name).sel(crop) filename = f"topo_{name}.nc" da.to_dataset().to_netcdf(outdir / filename) + logger.info(f"saved: {outdir / filename}") def scale_to_pixel(scales, dem_da): @@ -90,7 +91,7 @@ def scale_to_pixel(scales, dem_da): x_coords, y_coords = dem_da["x"].values, dem_da["y"].values epsg_code = dem_da.attrs["crs"].lower() if "epsg:4326" in epsg_code: - logger.warning( + logger.debug( f"Reprojecting coordinates from WGS84 to UTM to obtain units of meters" ) x_coords, y_coords = np.meshgrid(x_coords, y_coords) @@ -101,6 +102,7 @@ def scale_to_pixel(scales, dem_da): x_res = np.gradient(x_coords, axis=n_dims - 1) y_res = np.gradient(y_coords, axis=0) mean_res = np.mean(np.abs([x_res.mean(), y_res.mean()])) + logger.debug(f"Estimated resolution: {mean_res:.0f} meters.") return round_up_to_odd(np.array(scales) / mean_res), {"x": x_res, "y": y_res} diff --git a/topo_descriptors/topo.py b/topo_descriptors/topo.py index c68d258..11204ef 100755 --- a/topo_descriptors/topo.py +++ b/topo_descriptors/topo.py @@ -104,8 +104,8 @@ def tpi(dem, size, sigma=None): if isinstance(dem.data, da.Array): conv = da.map_overlap(conv_fn, dem.data, depth=size * 2, boundary="none") - elif isinstance(dem.data, np.ndarray): - conv = conv_fn(dem.values) + else: + conv = conv_fn(dem) return dem - conv / np.sum(kernel) @@ -136,7 +136,7 @@ def circular_kernel(size): else: xx, yy = np.mgrid[:size, :size] circle = (xx - middle) ** 2 + (yy - middle) ** 2 - kernel = np.asarray(circle <= (middle ** 2), dtype=np.float32) + kernel = np.asarray(circle <= (middle**2), dtype=np.float32) return kernel @@ -228,7 +228,7 @@ def std(dem, size, sigma=None): sum_dem = signal.convolve(dem, kernel, mode="same") sum_squared_dem = signal.convolve(squared_dem, kernel, mode="same") - variance = (sum_squared_dem - sum_dem ** 2 / kernel_sum) / (kernel_sum - 1) + variance = (sum_squared_dem - sum_dem**2 / kernel_sum) / (kernel_sum - 1) variance = np.clip(variance, 0, None) # avoid small negative values return np.sqrt(variance) @@ -585,7 +585,7 @@ def gradient(dem, sigma, res_meters, sig_ratio=1): _normalize_dxy(dx, dy, res_meters) - slope = np.sqrt(dx ** 2, dy ** 2) + slope = np.sqrt(dx**2, dy**2) aspect = ( 180 + np.degrees(np.arctan2(dx, dy)) ) % 360 # north faces = 0°, east faces = 90° @@ -909,7 +909,7 @@ def _sx_rolling(dem, distance, blines, height): def _sx_name(radius, azimuth): """Return name for the array in output of the Sx function""" - add = f"_RADIUS{int(radius[0])}-{int(radius)}_AZIMUTH{int(azimuth)}" + add = f"_RADIUS{int(radius)}_AZIMUTH{int(azimuth)}" return f"SX{add}" diff --git a/tpi_500M.png b/tpi_500M.png deleted file mode 100755 index 397ea05..0000000 Binary files a/tpi_500M.png and /dev/null differ