From d5cf87efac105cb528e651febb398913f821aadf Mon Sep 17 00:00:00 2001 From: Alen Date: Fri, 23 Aug 2024 15:07:27 +0200 Subject: [PATCH 1/6] lidar tutorial --- jupyter_lidar_inspection.ipynb | 918 +++++++++++++++++++++++++++++++++ 1 file changed, 918 insertions(+) create mode 100644 jupyter_lidar_inspection.ipynb diff --git a/jupyter_lidar_inspection.ipynb b/jupyter_lidar_inspection.ipynb new file mode 100644 index 0000000..589db88 --- /dev/null +++ b/jupyter_lidar_inspection.ipynb @@ -0,0 +1,918 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1cd6b205", + "metadata": {}, + "source": [ + "# LiDAR point clouds with GRASS GIS & Python in Jupyter Notebooks" + ] + }, + { + "cell_type": "markdown", + "id": "2cc05b0d", + "metadata": {}, + "source": [ + "This tutorial is meant to show you how to check LiDAR point clouds with GRASS in python. When we get new airborne data we have to see what are we dealing with. Different sensors, planes, height of acqusiton, season, plant cover, weather and in general different situations are impacting LiDAR data. Before producing elevation models and before all the fun (simulatons and maps) we need to inspect the data.\n", + "\n", + "Let's get started!\n", + "\n", + "This tutorial can be run locally. You need to have GRASS GIS 8.4+ and Jupyter installed. The processes demands some specific libraries. Be sure to have numpy, pandas, geopandas, pdal, python-pdal and tiledbb with pybabylonjs. \n", + "\n", + "The first thing we need to do for any of the cases we'll see further on, is to import GRASS GIS python packages. In order to do so, we need to add GRASS GIS python package to PATH. Let's see how we do that." + ] + }, + { + "cell_type": "markdown", + "id": "1b88ac81", + "metadata": {}, + "source": [ + "## Setup" + ] + }, + { + "cell_type": "markdown", + "id": "26a46471-616d-4ff6-bcea-ff2f56bd52a6", + "metadata": {}, + "source": [ + "Install micromamba and setup a new invironment. It works aso with conda and mamba." + ] + }, + { + "cell_type": "markdown", + "id": "e590a2de", + "metadata": {}, + "source": [ + "`\"${SHELL}\" <(curl -L micro.mamba.pm/install.sh)`
\n", + "`micromamba create -n grass_lidar numpy geopandas pdal python-pdal laspy tiledb-py pyarrow ipywidgets==7.7.2 jupyterlab==3.4.5 wxpython`" + ] + }, + { + "cell_type": "markdown", + "id": "0b6626bd-537b-423d-ae7c-533592fc6faf", + "metadata": {}, + "source": [ + "Import Python standard library and IPython packages we need." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f7a0d3d", + "metadata": {}, + "outputs": [], + "source": [ + "# pip install pybabylonjs\n", + "import os\n", + "import subprocess\n", + "import sys\n", + "import pdal\n", + "#pip install pdal python-pdal\n", + "import tiledb\n", + "#pip install tiledb-py pyarrow pandas pdal python-pdal\n", + "import laspy\n", + "#pip install lazrs\n", + "import pandas as pd\n", + "import geopandas as gpd\n", + "import numpy as np\n", + "import pybabylonjs" + ] + }, + { + "cell_type": "markdown", + "id": "0d020acc-fb86-4a2f-ac39-9f36fcba47f8", + "metadata": {}, + "source": [ + "Check where GRASS GIS python packages are and add them to PATH." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "569269bc", + "metadata": {}, + "outputs": [], + "source": [ + "sys.path.append(\n", + " subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "7e2346f7-293a-44e4-b201-6d34c99368fd", + "metadata": {}, + "source": [ + "Import the GRASS GIS libraries we need." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aba908f3", + "metadata": {}, + "outputs": [], + "source": [ + "import grass.script as gs\n", + "import grass.jupyter as gj\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "markdown", + "id": "0a70d98a-812f-49a1-9f58-4b36784a2aa9", + "metadata": {}, + "source": [ + "We can check the available commands with this command" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a54d7060", + "metadata": {}, + "outputs": [], + "source": [ + "gs.get_commands()" + ] + }, + { + "cell_type": "markdown", + "id": "cc67099a", + "metadata": {}, + "source": [ + "## Project Initialization and Import" + ] + }, + { + "cell_type": "markdown", + "id": "884dbab0-cbb5-4a0e-a269-eb9b4722b65a", + "metadata": {}, + "source": [ + "Create a temporary folder where to place our GRASS project." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3665f24", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "tempdir = tempfile.TemporaryDirectory()\n", + "# Create a project in the temporary directory\n", + "gs.create_project(path=tempdir.name, name=\"lidar_inspection\", epsg=\"3794\", overwrite=True)\n", + "# Start GRASS in the recently created project\n", + "session = gj.init(Path(tempdir.name,\"lidar_inspection\"))" + ] + }, + { + "cell_type": "markdown", + "id": "7d05141c-46ee-4946-806b-6d6b921fa6df", + "metadata": {}, + "source": [ + "Or create/connect a project in an existing GRASS database." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eaefd740", + "metadata": {}, + "outputs": [], + "source": [ + "gisdbase = '/your/path/grassdata'\n", + "gs.create_project(path=gisdbase, name=\"lidar_inspection\", epsg=\"3794\", overwrite=False)\n", + "session = gj.init(Path(gisdbase,\"lidar_inspection\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f67c8b8f", + "metadata": {}, + "outputs": [], + "source": [ + "# Let's check the environment settings: with os.environ[] we can see all of them\n", + "os.environ['GRASS_OVERWRITE'] # This environent variable shows us if we can overwrite our outputs. " + ] + }, + { + "cell_type": "markdown", + "id": "71e59b8d", + "metadata": {}, + "source": [ + "- Question: What can wee see from the output above? Is it the same as we are used in GRASS?" + ] + }, + { + "cell_type": "markdown", + "id": "e3b2f12d-e1f5-402e-b1a5-e72b94e7c017", + "metadata": {}, + "source": [ + "We can also use PyGRASS. For example, like this we get the documentation in Python." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "64a049ae", + "metadata": {}, + "outputs": [], + "source": [ + "from grass.pygrass.modules import Module\n", + "r_in_pdal = Module(\"r.in.pdal\")\n", + "print(r_in_pdal.__doc__) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4dc3bb3", + "metadata": {}, + "outputs": [], + "source": [ + "# Set the location of our point cloud\n", + "point_cloud = \"/your/path/point_cloud.laz\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51cadb94", + "metadata": {}, + "outputs": [], + "source": [ + "# Check our point cloud\n", + "print(gs.read_command(\"r.in.pdal\", \n", + " input=point_cloud,\n", + " flags=\"p\"))" + ] + }, + { + "cell_type": "markdown", + "id": "b622c5e3-a32e-4a85-9db8-7a29e1dec3ef", + "metadata": {}, + "source": [ + "Know your data format!\n", + "https://www.asprs.org/wp-content/uploads/2019/07/LAS_1_4_r15.pdf" + ] + }, + { + "cell_type": "markdown", + "id": "ab824297", + "metadata": {}, + "source": [ + "- Question 1: What is the difference between ReturnNumber and NumberOfReturns?\n", + "- Question 2: What is the classification flag value have points classified as ground?\n", + "- Question 3: What is the difference between the Point Format 7 and 8 and what can dimensions Red, Green, Blue, Infrared contain?" + ] + }, + { + "cell_type": "markdown", + "id": "ddb495cb-162f-4f41-9333-20d9727ae697", + "metadata": {}, + "source": [ + "Another way to do it using PyGRASS." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d66466a", + "metadata": {}, + "outputs": [], + "source": [ + "from grass.pygrass.modules import Module\n", + "Module(\"r.in.pdal\", \n", + " input=point_cloud,\n", + " flags=\"p\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a9e2ce0a", + "metadata": {}, + "outputs": [], + "source": [ + "# See computational region\n", + "gs.run_command(\"g.region\",flags=\"p\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61b147de", + "metadata": {}, + "outputs": [], + "source": [ + "# Set the resolution\n", + "gs.run_command(\"g.region\",res=1) " + ] + }, + { + "cell_type": "markdown", + "id": "ffbbbd92", + "metadata": {}, + "source": [ + "We are interested to see how different filters, gridding (resolution) and methods affects our raszerization\n", + "Our point cloud has no CRS defined but it fits our project projection, se we overrired the projection chech with `-o`" + ] + }, + { + "cell_type": "markdown", + "id": "c394f189-30fe-4bc3-b9a6-90570e1d9733", + "metadata": {}, + "source": [ + "Import our point cloud as a raster DSM, so we can check it a bit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4078e3f3", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "# With -e we use the extent of the input for the raster extent \n", + "# With -n we set the computation region to match the new raster map\n", + "\n", + "gs.run_command(\"r.in.pdal\",\n", + " input=point_cloud,\n", + " output='point_cloud_1m',\n", + " resolution=1,\n", + " flags=\"eno\")\n", + "\n", + "gs.run_command(\"r.in.pdal\",\n", + " input=point_cloud,\n", + " output='point_cloud_1m_last',\n", + " resolution=1,\n", + " return_filter='last', # Let's see what points are last\n", + " flags=\"eo\")\n", + "\n", + "gs.run_command(\"r.in.pdal\",\n", + " input=point_cloud,\n", + " output='point_cloud_1m_ground',\n", + " resolution=1,\n", + " class_filter=2, # Let's see what points are classified as ground\n", + " flags=\"eo\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fcfc296f", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a list of the rasters we produced\n", + "point_cloud_check=gs.list_strings(type=['raster'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93346e2b", + "metadata": {}, + "outputs": [], + "source": [ + "series = gj.SeriesMap(height = 700)\n", + "series.add_rasters(point_cloud_check)\n", + "series.d_grid(size=250,color='orange') # With a grid will be easier to focus later on checking specific parts\n", + "series.show()\n", + "\"series.save(\"/your/path/image.gif\") # You can export that into a GIF, add it to your presentation and be cool" + ] + }, + { + "cell_type": "markdown", + "id": "d3b999c1", + "metadata": {}, + "source": [ + "Set a smaller computational region for checking how would we like the data processed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b572a576", + "metadata": {}, + "outputs": [], + "source": [ + "gs.run_command(\"g.region\",save='my_region',flags=\"s\") # Let's save our region in the GRASS way\n", + "region_dict = gs.parse_command(\"g.region\",flags=\"gu\") # Tthe parse_command parses the region parameters into a dictionary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90b07d6b", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a copy of the original region dictionary\n", + "smaller_region_dict = region_dict.copy()\n", + "\n", + "# Update specific keys in the new dictionary\n", + "# Set the values of the extent you want to examine\n", + "updates = {\n", + " 'n': '79000',\n", + " 's': '78750',\n", + " 'w': '478250',\n", + " 'e': '478500'\n", + "}\n", + "# Apply the updates to the new dictionary\n", + "smaller_region_dict.update(updates)" + ] + }, + { + "cell_type": "markdown", + "id": "d7839235-f132-49bf-bee0-fc2490f07f09", + "metadata": {}, + "source": [ + "Let's update momentarily our region so we can visualize a smaller portion of it" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe7dc980", + "metadata": {}, + "outputs": [], + "source": [ + "n=int(smaller_region_dict['n'])\n", + "s=int(smaller_region_dict['s'])\n", + "w=int(smaller_region_dict['w'])\n", + "e=int(smaller_region_dict['e'])\n", + "univar_stat=gs.parse_command('r.univar',map='point_cloud_1m', flags='g')\n", + "z_min=float(univar_stat['min'])\n", + "z_max=float(univar_stat['max'])\n", + "\n", + "gs.run_command(\"g.region\",\n", + " save='smaller_region',\n", + " n=n,\n", + " s=s,\n", + " w=w,\n", + " e=e,\n", + " flags=\"d\")\n", + "gs.run_command(\"g.region\",region='smaller_region')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5b70bf8", + "metadata": {}, + "outputs": [], + "source": [ + "gs.run_command(\"g.region\",\n", + " save='smaller_region',\n", + " n=int(smaller_region_dict['n']),\n", + " s=int(smaller_region_dict['s']),\n", + " w=int(smaller_region_dict['w']),\n", + " e=int(smaller_region_dict['e']),\n", + " flags=\"d\")\n", + "gs.run_command(\"g.region\",region='smaller_region')" + ] + }, + { + "cell_type": "markdown", + "id": "6e356497-1662-443c-9de4-689e139a8c10", + "metadata": {}, + "source": [ + "## 3D Data Visualization" + ] + }, + { + "cell_type": "markdown", + "id": "755e62d4", + "metadata": {}, + "source": [ + "Let's do a 3D check of the region: tiledb will manage that our computer survives" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d4815a0", + "metadata": {}, + "outputs": [], + "source": [ + "data = os.path.expanduser(\"/your/path/point_cloud.laz\")\n", + "array_name = os.path.expanduser(\"/your/path/tiledb_pc\")\n", + "filename = \"/your/path/point_cloud.laz\"\"" + ] + }, + { + "cell_type": "markdown", + "id": "114bc70b-a904-4c2f-8331-e585eda9b8e4", + "metadata": {}, + "source": [ + "Let's create a PDAL pipeline to import our LAZ into a tiledb database." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d84f2b86", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "pipeline = pdal.Reader.las(filename=data).pipeline()\n", + "pipeline |= pdal.Writer.tiledb(filename=array_name, x_tile_size=200, y_tile_size=200, z_tile_size=100)\n", + "pipeline.execute() # once we create it, we can just load it next time\n", + "# The output states the number of points" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aec2215e", + "metadata": {}, + "outputs": [], + "source": [ + "from pybabylonjs import Show as show\n", + "A = tiledb.open(array_name)\n", + "df = A.query(attrs=('Red', 'Green','Blue')).df[w:e, s:n, z_min:z_max]\n", + "# df = A.query(attrs=('Infrared', 'Green','Red')).df[478000+500:478000+700, 77999+500:77999+700, 406.14:593.856000]\n", + "data = {\n", + " 'X': df['X'],\n", + " 'Y': df['Y'],\n", + " 'Z': df['Z'],\n", + " 'Red': df['Red'], # or df['Infrared'], df['Green'], df['Red']\n", + " 'Green': df['Green'],\n", + " 'Blue': df['Blue']\n", + "}\n", + "\n", + "show.point_cloud(source=\"dict\",\n", + " data=data,\n", + " point_size = 4,\n", + " rgb_max=65535,\n", + " width = 1200,\n", + " height = 700)" + ] + }, + { + "cell_type": "markdown", + "id": "c9aa1f54", + "metadata": {}, + "source": [ + "- Question 1: What is the difference between a DSM and a DTM?\n", + "- Question 2: What resolution, method and filters would you use to generate the DTM and what to generate a DSM with vegetation classes?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b38d9f1", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "# Let's see what check how do different resolutions look\n", + "resolutions = np.arange(0.50, 3.25, 0.25) # We set steps of 0.25m from 0.5 to 3m\n", + "\n", + "# Create a list the will later be populated with outputs with different resolustions\n", + "rasters_different_res = []\n", + " \n", + "# Iterate over the list of resolutions\n", + "for resolution in resolutions:\n", + " # Format the name of the output\n", + " resolution_str = f\"{resolution:.2f}\".replace('.', '')\n", + "\n", + " # Create the output name\n", + " output_name = f\"point_cloud_{resolution_str}\"\n", + " \n", + " # Run the command with the current resolution and output name\n", + " gs.run_command(\"r.in.pdal\",\n", + " input=point_cloud,\n", + " output=output_name,\n", + " resolution=resolution,\n", + " class_filter=2, # Filter points classified as ground\n", + " flags=\"eo\")\n", + " \n", + " # Add the output name to the list of rasters\n", + " rasters_different_res.append(output_name)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eeeaf121", + "metadata": {}, + "outputs": [], + "source": [ + "series = gj.SeriesMap(height = 700)\n", + "series.add_rasters(rasters_different_res)\n", + "series.show()" + ] + }, + { + "cell_type": "markdown", + "id": "6664c8fd", + "metadata": {}, + "source": [ + "What gridding would you chouse to create a DTM?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "09e3874e", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "# We use them all! Let patch them hierarchically and then interpolate.\n", + "gs.run_command(\"g.region\",res=0.5)\n", + "gs.run_command(\"r.patch\",input=rasters_different_res,output='dtm_patched',nprocs=8,mem=2000)\n", + "gs.run_command(\"r.relief\",input='dtm_patched',output='dtm_patched_hs',altitude=30,azimuth=315)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c2234ba", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "# List to store the names of the output shaded relief rasters\n", + "rasters_different_res_hs = []\n", + "\n", + "# Iterate over each raster and generate the shaded relief\n", + "for raster in rasters_different_res:\n", + " input_raster = raster\n", + " output_name = f\"{input_raster}_hs\"\n", + "\n", + " # Run r.relief\n", + " gs.run_command(\n", + " \"r.relief\",\n", + " input=input_raster,\n", + " output=output_name,\n", + " altitude=30,\n", + " azimuth=315\n", + " )\n", + "\n", + " # Add the output name to the list of rasters\n", + " rasters_different_res_hs.append(output_name)\n", + "\n", + "rasters_different_res_hs.append('dtm_patched_hs')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8e04327", + "metadata": {}, + "outputs": [], + "source": [ + "series = gj.SeriesMap(height = 700)\n", + "series.add_rasters(rasters_different_res_hs)\n", + "series.show()\n", + "# You can see here that better resolution doesn't mean more detail." + ] + }, + { + "cell_type": "markdown", + "id": "4ccb7f8b", + "metadata": {}, + "source": [ + "Is that composite better or worse than an interpolated, raw one for your future data analysis?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bfce78e", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "gs.run_command('r.resamp.bspline',\n", + " input='point_cloud_100', # I decided that this one has a good resolution:detail payoff\n", + " output='point_cloud_100_intp',\n", + " method='bicubic')\n", + "gs.run_command(\"r.relief\",input='point_cloud_100_intp',output='point_cloud_100_intp_hs',altitude=30,azimuth=315)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08833664", + "metadata": {}, + "outputs": [], + "source": [ + "series = gj.SeriesMap(height = 700)\n", + "series.add_rasters(['dtm_patched_hs','point_cloud_100_intp_hs'])\n", + "series.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6693fdfc", + "metadata": {}, + "outputs": [], + "source": [ + "m = gj.InteractiveMap(width=700, height=700, tiles=None)\n", + "m.add_raster('dtm_patched_hs')\n", + "m.add_raster('point_cloud_100_intp_hs')\n", + "m.add_layer_control()\n", + "m.show()" + ] + }, + { + "cell_type": "markdown", + "id": "9e2b55e8", + "metadata": {}, + "source": [ + "Let's first interpolate the missing values in the patched DTM and then dig a bit more." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97706b1c", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "gs.run_command('r.resamp.bspline', \n", + " input='dtm_patched',\n", + " output='dtm_patched_int',\n", + " method='bicubic')" + ] + }, + { + "cell_type": "markdown", + "id": "24cdeade", + "metadata": {}, + "source": [ + "Guess, what: GRASS has also a ton of addons and we can install them easily also here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0d9ea7d", + "metadata": {}, + "outputs": [], + "source": [ + "# We can see differences on microtopography:\n", + "# Let's denoise and see how does local noise affect the elevation mode\n", + "\n", + "# Let's install the extension\n", + "gs.run_command(\"g.extension\",\n", + " extension='r.denoise',\n", + " operation='add')" + ] + }, + { + "cell_type": "markdown", + "id": "fd959956", + "metadata": {}, + "source": [ + "Create a temporary directory or if you created one in the beginning just use that." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7e5871a", + "metadata": {}, + "outputs": [], + "source": [ + "# Let's compile the requirements of r.denoise\n", + "import tempfile\n", + "temp_dir = tempfile.mkdtemp() # If you created a new one in the beginning just use that\n", + "repo_url = \"https://github.com/exuberant/mdenoise.git\" # Life changes, repos change: check it\n", + "subprocess.run([\"git\", \"clone\", repo_url, temp_dir], check=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "94787d40", + "metadata": {}, + "outputs": [], + "source": [ + "# Change directory to the cloned repo\n", + "os.chdir(temp_dir)\n", + "\n", + "# Compile the code (if needed)\n", + "compile_cmd = \"g++ -o mdenoise mdenoise.cpp triangle.c\"\n", + "subprocess.run(compile_cmd, shell=True, check=True)\n", + "\n", + "print(\"Yep, no errors.\")" + ] + }, + { + "cell_type": "markdown", + "id": "f27e9c41", + "metadata": {}, + "source": [ + "Let's denoise and see how does local noise affect the elevation mode." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7fcf4a5c", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "thresholds = np.arange(0.85, 1.0, 0.02)\n", + "\n", + "# Create the list the will later be populated with outputs with different resolustions\n", + "dtm_different_denoise = []\n", + " \n", + "# Iterate over the list of resolutions\n", + "for threshold in thresholds:\n", + " # Format the name of the output\n", + " threshold_str = f\"{threshold:.2f}\".replace('.', '')\n", + "\n", + " # Create the output name\n", + " output_name = f\"dtm_denoise_{threshold_str}\"\n", + " \n", + " # Run the command with the current resolution and output name\n", + " gs.run_command(\"r.denoise\",\n", + " input='dtm_patched_int',\n", + " output=output_name,\n", + " iterations=5,\n", + " threshold=threshold)\n", + " \n", + " # Add the output name to the list of rasters\n", + " dtm_different_denoise.append(output_name)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef66ab4e", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "# List to store the names of the output shaded relief rasters\n", + "dtm_different_denoise_hs = []\n", + "\n", + "# Iterate over each raster and generate the shaded relief\n", + "for dtm in dtm_different_denoise:\n", + " input_raster = dtm\n", + " output_name = f\"{input_raster}_hs\"\n", + "\n", + " # Run r.relief\n", + " gs.run_command(\n", + " \"r.relief\",\n", + " input=input_raster,\n", + " output=output_name,\n", + " altitude=30,\n", + " azimuth=315\n", + " )\n", + "\n", + " # Add the output name to the list of rasters\n", + " dtm_different_denoise_hs.append(output_name)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c91e6349", + "metadata": {}, + "outputs": [], + "source": [ + "series = gj.SeriesMap(height = 700)\n", + "series.add_rasters(dtm_different_denoise_hs)\n", + "series.show()\n", + "# Pick the one that you need!" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.12.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 3cc793e06a1e90fb4d3ab8ce7fb99d47c7c9f9fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alen=20Mangafi=C4=87?= Date: Thu, 24 Oct 2024 11:06:47 +0200 Subject: [PATCH 2/6] Update jupyter_lidar_inspection.ipynb Co-authored-by: Veronica Andreo --- jupyter_lidar_inspection.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/jupyter_lidar_inspection.ipynb b/jupyter_lidar_inspection.ipynb index 589db88..e7205e3 100644 --- a/jupyter_lidar_inspection.ipynb +++ b/jupyter_lidar_inspection.ipynb @@ -35,7 +35,7 @@ "id": "26a46471-616d-4ff6-bcea-ff2f56bd52a6", "metadata": {}, "source": [ - "Install micromamba and setup a new invironment. It works aso with conda and mamba." + "Install micromamba and setup a new environment. It works also with conda and mamba." ] }, { From 18efde570cacc4af384085b8f76b4d7e200fcd40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alen=20Mangafi=C4=87?= Date: Thu, 24 Oct 2024 11:06:59 +0200 Subject: [PATCH 3/6] Update jupyter_lidar_inspection.ipynb Co-authored-by: Veronica Andreo --- jupyter_lidar_inspection.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/jupyter_lidar_inspection.ipynb b/jupyter_lidar_inspection.ipynb index e7205e3..c6111b0 100644 --- a/jupyter_lidar_inspection.ipynb +++ b/jupyter_lidar_inspection.ipynb @@ -123,7 +123,7 @@ "id": "0a70d98a-812f-49a1-9f58-4b36784a2aa9", "metadata": {}, "source": [ - "We can check the available commands with this command" + "We can check the available tools with this command" ] }, { From 52056f132a3d224932b9285171af1479c590d992 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alen=20Mangafi=C4=87?= Date: Thu, 24 Oct 2024 11:07:06 +0200 Subject: [PATCH 4/6] Update jupyter_lidar_inspection.ipynb Co-authored-by: Veronica Andreo --- jupyter_lidar_inspection.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/jupyter_lidar_inspection.ipynb b/jupyter_lidar_inspection.ipynb index c6111b0..038efe0 100644 --- a/jupyter_lidar_inspection.ipynb +++ b/jupyter_lidar_inspection.ipynb @@ -172,7 +172,7 @@ "id": "7d05141c-46ee-4946-806b-6d6b921fa6df", "metadata": {}, "source": [ - "Or create/connect a project in an existing GRASS database." + "Or optionally create/connect a project in an existing GRASS database." ] }, { From ae96926427361ac3c2d870b9062330c9d402278e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alen=20Mangafi=C4=87?= Date: Thu, 24 Oct 2024 11:08:38 +0200 Subject: [PATCH 5/6] Update jupyter_lidar_inspection.ipynb Co-authored-by: Veronica Andreo --- jupyter_lidar_inspection.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/jupyter_lidar_inspection.ipynb b/jupyter_lidar_inspection.ipynb index 038efe0..5195bcd 100644 --- a/jupyter_lidar_inspection.ipynb +++ b/jupyter_lidar_inspection.ipynb @@ -627,7 +627,7 @@ "outputs": [], "source": [ "%%time\n", - "# We use them all! Let patch them hierarchically and then interpolate.\n", + "# We'll use them all! Let's patch them hierarchically and then interpolate.\n", "gs.run_command(\"g.region\",res=0.5)\n", "gs.run_command(\"r.patch\",input=rasters_different_res,output='dtm_patched',nprocs=8,mem=2000)\n", "gs.run_command(\"r.relief\",input='dtm_patched',output='dtm_patched_hs',altitude=30,azimuth=315)" From cb5aac5b27b92f679a16de99f75bb94af7fa4760 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alen=20Mangafi=C4=87?= Date: Thu, 24 Oct 2024 11:08:56 +0200 Subject: [PATCH 6/6] Update jupyter_lidar_inspection.ipynb Co-authored-by: Veronica Andreo --- jupyter_lidar_inspection.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/jupyter_lidar_inspection.ipynb b/jupyter_lidar_inspection.ipynb index 5195bcd..9b94e31 100644 --- a/jupyter_lidar_inspection.ipynb +++ b/jupyter_lidar_inspection.ipynb @@ -816,7 +816,7 @@ "id": "f27e9c41", "metadata": {}, "source": [ - "Let's denoise and see how does local noise affect the elevation mode." + "Let's denoise and see how local noise affects the elevation model." ] }, {