{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Description: GRASS GIS offers, besides other things, numerous analytical tools for point clouds, terrain, and remote sensing. In this hands-on workshop we will explore the tools in GRASS GIS for processing point clouds obtained by lidar or through processing of UAV imagery. We will start with a brief and focused introduction into GRASS GIS graphical user interface (GUI) and we will continue with short introduction to GRASS GIS Python interface. Participants will then decide if they will use GUI, command line, Python, or online Jupyter Notebook for the rest of the workshop. We will explore the properties of the point cloud, interpolate surfaces, and perform advanced terrain analyses to detect landforms and artifacts. We will go through several terrain 2D and 3D visualization techniques to get more information out of the data and finish with vegetation analysis.\n", "Requirements: This workshop is accessible to beginners, but some basic knowledge of lidar processing or GIS is helpful for a smooth experience.\n", "Authors: [Vaclav Petras](/wiki/User:Wenzeslaus), [Anna Petrasova](/wiki/User:Annakrat), and Helena Mitasova from North Carolina State University\n", "Contributors: Robert S. Dzur and Doug Newcomb\n", "\n", "## Contents\n", "\n", "* [1 Preparation](#Preparation)\n", "\n", "* [1.1 Software](#Software)\n", "* [1.2 Addons](#Addons)\n", "* [1.3 Data](#Data)\n", "\n", "\n", "* [2 Basic introduction to graphical user interface](#Basic_introduction_to_graphical_user_interface)\n", "\n", "* [2.1 GRASS GIS Spatial Database](#GRASS_GIS_Spatial_Database)\n", "* [2.2 Importing data](#Importing_data)\n", "* [2.3 Computational region](#Computational_region)\n", "* [2.4 Modules](#Modules)\n", "\n", "\n", "* [3 Basic introduction to Python interface](#Basic_introduction_to_Python_interface)\n", "* [4 Decide if to use GUI, command line, Python, or online Jupyter Notebook](#Decide_if_to_use_GUI.2C_command_line.2C_Python.2C_or_online_Jupyter_Notebook)\n", "* [5 Binning of the point cloud](#Binning_of_the_point_cloud)\n", "* [6 Interpolation](#Interpolation)\n", "* [7 Terrain analysis](#Terrain_analysis)\n", "* [8 Vegetation analysis](#Vegetation_analysis)\n", "* [9 Bonus tasks](#Bonus_tasks)\n", "\n", "* [9.1 3D visualization of rasters](#3D_visualization_of_rasters)\n", "* [9.2 Analytical visualization of rasters in 3D](#Analytical_visualization_of_rasters_in_3D)\n", "* [9.3 Classify ground and non-ground points](#Classify_ground_and_non-ground_points)\n", "* [9.4 3D visualization of point clouds](#3D_visualization_of_point_clouds)\n", "* [9.5 Alternative DSM creation](#Alternative_DSM_creation)\n", "* [9.6 Explore layers of vegetation in a 3D raster](#Explore_layers_of_vegetation_in_a_3D_raster)\n", "\n", "\n", "* [10 Optimizations, troubleshooting and limitations](#Optimizations.2C_troubleshooting_and_limitations)\n", "* [11 Bleeding edge](#Bleeding_edge)\n", "* [12 See also](#See_also)\n", "\n", "\n", "\n", "## Preparation\n", "### Software\n", "GRASS GIS 7.2 compiled with libLAS is needed (e.g. [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) should work).\n", "OSGeo-Live\n", "All needed software is included in [OSGeo-Live](http://live.osgeo.org/).\n", " Ubuntu \n", "Install GRASS GIS from packages:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# execute manually the following or its equivalent:\n", "# sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable\n", "# execute manually the following or its equivalent:\n", "# sudo apt-get update\n", "# execute manually the following or its equivalent:\n", "# sudo apt-get install grass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Linux \n", "For other Linux distributions other then Ubuntu, please try to find GRASS GIS in their package managers.\n", "MS Windows\n", "Download the standalone GRASS GIS binaries from [grass.osgeo.org](https://grass.osgeo.org/).\n", "Mac OS\n", "Install GRASS GIS using [Homebrew](https://brew.sh/) [osgeo4mac](https://github.com/OSGeo/homebrew-osgeo4mac):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# execute manually the following or its equivalent:\n", "# brew tap osgeo/osgeo4mac\n", "# execute manually the following or its equivalent:\n", "# brew install numpy\n", "# execute manually the following or its equivalent:\n", "# brew install liblas --build-from-source\n", "# execute manually the following or its equivalent:\n", "# brew install grass7 --with-liblas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you have these already installed, you need to use `reinstall` instead of `install`.\n", "Note that on Mac OS in some versions the [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) module is not accessible, so you need to check this and use [r.in.ascii](https://grass.osgeo.org/grass72/manuals/r.in.ascii.html) in combination with libLAS or PDAL command line tools to achieve the same or, preferably, use OSGeo-Live.\n", "Note also that there is currently no recent version of installation on Mac OS where 3D view is accessible.\n", "\n", "### Addons\n", "You will need to install the following GRASS GIS Addons. This can be done through GUI, but for simplicity copy and paste and execute the following lines one by one in the command line:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# execute manually the following or its equivalent:\n", "# g.extension r.geomorphon\n", "# execute manually the following or its equivalent:\n", "# g.extension r.skyview\n", "# execute manually the following or its equivalent:\n", "# g.extension r.local.relief\n", "# execute manually the following or its equivalent:\n", "# g.extension r.shaded.pca\n", "# execute manually the following or its equivalent:\n", "# g.extension r.area\n", "# execute manually the following or its equivalent:\n", "# g.extension r.fill.gaps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For bonus tasks:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# execute manually the following or its equivalent:\n", "# g.extension v.lidar.mcc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Data\n", "We will need the two following ZIP files, downloaded and extracted:\n", "\n", "* [orthophoto](http://fatra.cnr.ncsu.edu/foss4g2017/nc_orthophoto_1m_spm.zip)\n", "* [point cloud](http://fatra.cnr.ncsu.edu/foss4g2017/nc_tile_0793_016_spm.zip)\n", "Later, for bonus task, download also this (much larger) file: \n", "\n", "* [UAV point cloud](http://fatra.cnr.ncsu.edu/foss4g2017/nc_uav_points_spm.zip)\n", "## Basic introduction to graphical user interface\n", "### GRASS GIS Spatial Database\n", "Here we provide an overview of [GRASS GIS](https://grass.osgeo.org). For this exercise it's not necessary to have a full understanding of how to use GRASS GIS. However, you will need to know how to place your data in the correct GRASS GIS database directory, as well as some basic GRASS functionality.\n", "GRASS uses specific database terminology and structure ([GRASS GIS Spatial Database](https://grass.osgeo.org/grass72/manuals/grass_database.html)) that are important to understand for working in GRASS GIS efficiently. You will create a new location and import the required data into that location. In the following we review important terminology and give step by step directions on how to download and place your data in the correct place.\n", "\n", "* A GRASS GIS Spatial Database (GRASS database) consists of directory with specific Locations (projects) where data (data layers/maps) are stored.\n", "* Location is a directory with data related to one geographic location or a project. All data within one Location has the same coordinate reference system.\n", "* Mapset is a collection of maps within Location, containing data related to a specific task, user or a smaller project.\n", "Start GRASS GIS, a start-up screen should appear.\n", "Unless you already have a directory called grassdata in your Documents directory (on MS Windows)\n", "or in your home directory (on Linux), create one. You can use the Browse button and the dialog in the GRASS GIS start up screen to do that.\n", "We will create a new location for our project with CRS (coordinate reference system) NC State Plane Meters with EPSG code 3358.\n", "Open Location Wizard with button New in the left part of the welcome screen. Select a name for the new location, select EPSG method and code 3358.\n", "When the wizard is finished, the new location will be listed on the start-up screen. Select the new location and mapset PERMANENT and press Start GRASS session.\n", "\n", "\n", "Note that a current working directory is a concept which is separate from the GRASS database, location and mapset discussed above. The current working directory is a directory where any program (no just GRASS GIS) writes and reads files unless a path to the file is provided. The current working directory can be changed from the GUI using Settings -> GRASS working environment -> Change working directory or from the Console using the cd command. This is advantageous when we are using command line and working with the same file again and again. This is often the case when working with lidar data. We can change the directory to the directory with the downloaded LAS file. IN case we don't change the directory, we need to provide full path to the file. Note that command line and GUI have their own settings of current working directory, so it needs to be changed separately for each of them." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# This is a quick introduction into Jupyter Notebook.\n", "# Python code can be excecuted like this:\n", "a = 6\n", "b = 7\n", "c = a * b\n", "print \"Answer is\", c\n", "# Python code can be mixed with command line code (Bash).\n", "# It is enough just to prefix the command line with an exclamation mark:\n", "!echo \"Answer is $c\"\n", "# Use Shift+Enter to execute this cell. The result is below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "import subprocess\n", "from IPython.display import Image\n", "\n", "# create GRASS GIS runtime environment\n", "gisbase = subprocess.check_output([\"grass\", \"--config\", \"path\"]).strip()\n", "os.environ['GISBASE'] = gisbase\n", "sys.path.append(os.path.join(gisbase, \"etc\", \"python\"))\n", "\n", "# do GRASS GIS imports\n", "import grass.script as gs\n", "import grass.script.setup as gsetup\n", "\n", "# set GRASS GIS session data\n", "rcfile = gsetup.init(gisbase, \"/home/jovyan/grassdata/\", \"workshop\", \"PERMANENT\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# default font displays\n", "os.environ['GRASS_FONT'] = 'sans'\n", "# overwrite existing maps\n", "os.environ['GRASS_OVERWRITE'] = '1'\n", "gs.set_raise_on_error(True)\n", "gs.set_capture_stderr(True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# set display modules to render into a file (named map.png by default)\n", "os.environ['GRASS_RENDER_IMMEDIATE'] = 'cairo'\n", "os.environ['GRASS_RENDER_FILE_READ'] = 'TRUE'\n", "os.environ['GRASS_LEGEND_FILE'] = 'legend.txt'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Importing data\n", "In this step we import the provided data into GRASS GIS. In menu File - Import raster data select Common formats import and in the dialog browse to find the orthophoto file, change the name to ortho.tif, and click button Import. All the imported layers should be added to GUI automatically, if not, add them manually. Point clouds will be imported later in a different way as part of the analysis.\n", "\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Import_raster_7.2.1.png)\n", "\t\t\t\n", "Import raster data: select the orthophoto file and change the name to ortho\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "The equivalent command is:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.import', input=\"nc_orthophoto_1m_spm.tif\", output=\"ortho\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computational region\n", "Before we use a module to compute a new raster map, we must properly set the computational region. All raster computations will be performed in the specified extent and with the given resolution.\n", "Computational region is an important raster concept in GRASS GIS. In GRASS a computational region can be set, subsetting larger extent data for quicker testing of analysis or analysis of specific regions based on administrative units. We provide a few points to keep in mind when using the computational region function:\n", "\n", "* defined by region extent and raster resolution\n", "* applies to all raster operations\n", "* persists between GRASS sessions, can be different for different mapsets\n", "* advantages: keeps your results consistent, avoid clipping, for computationally demanding tasks set region to smaller extent, check that your result is good and then set the computational region to the entire study area and rerun analysis\n", "* run `g.region -p` or in menu Settings - Region - Display region to see current region settings\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Computational_region_two_rasters.png)\n", "\t\t\t\n", "Computational region concept: A raster with large extent (blue) is displayed as well as another raster with smaller extent (green). The computational region (red) is now set to match the smaller raster, so all the computations are limited to the smaller raster extent even if the input is the larger raster. (Not shown on the image: Also the resolution, not only the extent, matches the resolution of the smaller raster.)\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:WxGUI_set_region.png)\n", "\t\t\t\n", "Simple ways to set computational region from GUI: On the left, set region to match raster map. On the right, select the highlighted option and then set region by drawing rectangle.\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Wxgui_computational_region_set_from_raster.png)\n", "\t\t\t\n", "Set computational region (extent and resolution) to match a raster (Layers tab in the Layer Manager)\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "The numeric values of computational region can be checked using:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After executing the command you will get something like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# execute manually the following or its equivalent:\n", "# north: 220750\n", "# execute manually the following or its equivalent:\n", "# south: 220000\n", "# execute manually the following or its equivalent:\n", "# west: 638300\n", "# execute manually the following or its equivalent:\n", "# east: 639000\n", "# execute manually the following or its equivalent:\n", "# nsres: 1\n", "# execute manually the following or its equivalent:\n", "# ewres: 1\n", "# execute manually the following or its equivalent:\n", "# rows: 750\n", "# execute manually the following or its equivalent:\n", "# cols: 700\n", "# execute manually the following or its equivalent:\n", "# cells: 525000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Computational region can be set using a raster map:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', raster=\"ortho\", flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Resolution can be set separately using the `res` parameter of the [g.region](https://grass.osgeo.org/grass72/manuals/g.region.html) module. The units are the units of the current location, in our case meters. This can be done in the Resolution tab of the [g.region](https://grass.osgeo.org/grass72/manuals/g.region.html) dialog or in the command line in the following way (using also the `-a` flag to print the new values):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', res=\"3\", flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The new resolution may be slightly modified in this case to fit into the extent which we are not changing. However, often we want the resolution to be the exact value we provide and we are fine with a slight modification of the extent. That's what `-a` flag is for. On the other hand, if alignment with cells of a specific raster is important, `align` parameter can be used to request the alignment to that raster (regardless of the the extent).\n", "The following example command will use the extent from the raster named `ortho`, use resolution `5` meters, modify the extent to align it to this 5 meter resolution, and print the values of this new computational region settings:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', raster=\"ortho\", res=\"5\", flags='apg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Modules\n", "GRASS GIS functionality is available through modules (tools, functions, commands). Modules respect the following naming conventions:\n", "\n", "\n", "\n", "\n", "\n", " Prefix \n", " Function \n", " Example\n", "\n", "\n", " r. \n", " raster processing \n", " [r.mapcalc](https://grass.osgeo.org/grass72/manuals/r.mapcalc.html): map algebra\n", "\n", "\n", " v. \n", " vector processing \n", " [v.surf.rst](https://grass.osgeo.org/grass72/manuals/v.surf.rst.html): surface interpolation\n", "\n", "\n", " i. \n", " imagery processing \n", " [i.segment](https://grass.osgeo.org/grass72/manuals/i.segment.html): image segmentation\n", "\n", "\n", " r3. \n", " 3D raster processing \n", " [r3.stats](https://grass.osgeo.org/grass72/manuals/r3.stats.html): 3D raster statistics\n", "\n", "\n", " t. \n", " temporal data processing \n", " [t.rast.aggregate](https://grass.osgeo.org/grass72/manuals/t.rast.aggregate.html): temporal aggregation\n", "\n", "\n", " g. \n", " general data management \n", " [g.remove](https://grass.osgeo.org/grass72/manuals/g.remove.html): removes maps\n", "\n", "\n", " d. \n", " display \n", "[d.rast](https://grass.osgeo.org/grass72/manuals/d.rast.html): display raster map\n", "\n", "These are the main groups of modules. There is few more for specific purposes. Note also that some modules have multiple dots in their names. This often suggests further grouping. For example, modules staring with v.net. deal with vector network analysis. The name of the module helps to understand its function, for example v.in.lidar starts with v so it deals with vector maps, the name follows with in which indicates that the module is for importing the data into GRASS GIS Spatial Database and finally lidar indicates that it deals with lidar point clouds.\n", "\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:R_in_lidar_dialog.png)\n", "\t\t\t\n", "r.in.lidar dialog with Output tab active and highlighted module name (blue), options and flags (red) and option values (green)\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Example_r.in.lidar_command_in_Bash.png)\n", "\t\t\t\n", "Example r.in.lidar command in Bash with highlighted module name (blue), options and flags (red) and option values (green)\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Example_r.in.lidar_command_in_Python.png)\n", "\t\t\t\n", "Example r.in.lidar command in Python with highlighted module name (blue), options and flags (red) and option values (green) and import (grey)\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Graphical_Modeler_with_r.in.lidar_and_terrain_analysis.png)\n", "\t\t\t\n", "GRASS GIS Graphical Modeler\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:WxGUI_console_completion.png)\n", "\t\t\t\n", "Console tab in the Layer Manager for running commands or opening module dialogs\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Grass_gis_cli_ubuntu_purple_r.in.lidar.png)\n", "\t\t\t\n", "Command line interface (CLI) in a system terminal\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Wxgui_module_parameters_r_neighbors.png)\n", "\t\t\t\n", "Layout of a module dialog\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:WxGUI_module_search.png)\n", "\t\t\t\n", "Search for a module in the Modules tab\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:G_search_modules_with_c_flag.png)\n", "\t\t\t\n", "Search for a module using advanced search with [g.search.modules](https://grass.osgeo.org/grass72/manuals/g.search.modules.html)\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "One of the advantages of GRASS GIS is the diversity and number of modules that let you analyze all manners of spatial and temporal data. GRASS GIS has over [500 different modules](https://grass.osgeo.org/grass72/manuals/full_index.html) in the core distribution and over [300 addon modules](https://grass.osgeo.org/grass72/manuals/addons/) that can be used to prepare and analyze data. The following table lists some of the main modules for point cloud analysis. \n", "\n", "\n", "\n", "\n", "\n", " Module \n", " Function \n", " Alternatives\n", "\n", "\n", " [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) \n", " binning into 2D raster, statistics \n", " [r.in.xyz](https://grass.osgeo.org/grass72/manuals/r.in.xyz.html), [v.vect.stats](https://grass.osgeo.org/grass72/manuals/v.vect.stats.html), [r.vect.stats](https://grass.osgeo.org/grass7/manuals/addons/r.vect.stats.html)\n", "\n", "\n", " [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html) \n", " import, decimation \n", " [v.in.ascii](https://grass.osgeo.org/grass72/manuals/v.in.ascii.html), [v.in.ogr](https://grass.osgeo.org/grass72/manuals/v.in.ogr.html), [v.import](https://grass.osgeo.org/grass72/manuals/v.import.html)\n", "\n", "\n", " [r3.in.lidar](https://grass.osgeo.org/grass72/manuals/r3.in.lidar.html) \n", " binning into 3D raster \n", " [r3.in.xyz](https://grass.osgeo.org/grass72/manuals/r3.in.xyz.html), [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html)\n", "\n", "\n", " [v.out.lidar](https://grass.osgeo.org/grass72/manuals/v.out.lidar.html) \n", " export of point cloud \n", " [v.out.ascii](https://grass.osgeo.org/grass72/manuals/v.out.ascii.html), [r.out.xyz](https://grass.osgeo.org/grass72/manuals/r.out.xyz.html)\n", "\n", "\n", " [v.surf.rst](https://grass.osgeo.org/grass72/manuals/v.surf.rst.html) \n", " interpolation surfaces from points \n", " [v.surf.bspline](https://grass.osgeo.org/grass72/manuals/v.surf.bspline.html), [v.surf.idw](https://grass.osgeo.org/grass72/manuals/v.surf.idw.html)\n", "\n", "\n", " [v.lidar.edgedetection](https://grass.osgeo.org/grass72/manuals/v.lidar.edgedetection.html) \n", " ground and object (edge) detection \n", " [v.lidar.mcc](https://grass.osgeo.org/grass7/manuals/addons/v.lidar.mcc.html), [v.outlier](https://grass.osgeo.org/grass72/manuals/v.outlier.html)\n", "\n", "\n", " [v.decimate](https://grass.osgeo.org/grass72/manuals/v.decimate.html) \n", " decimate (thin) a point cloud \n", " [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html), [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html)\n", "\n", "\n", " [r.slope.aspect](https://grass.osgeo.org/grass72/manuals/r.slope.aspect.html) \n", " topographic parameters \n", " [v.surf.rst](https://grass.osgeo.org/grass72/manuals/v.surf.rst.html), [r.param.scale](https://grass.osgeo.org/grass72/manuals/r.param.scale.html)\n", "\n", "\n", " [r.relief](https://grass.osgeo.org/grass72/manuals/r.relief.html) \n", " shaded relief computation \n", " [r.skyview](https://grass.osgeo.org/grass7/manuals/addons/r.skyview.html), [r.local.relief](https://grass.osgeo.org/grass7/manuals/addons/r.local.relief.html)\n", "\n", "\n", " [r.colors](https://grass.osgeo.org/grass72/manuals/r.colors.html) \n", " raster color table management \n", " [r.cpt2grass](https://grass.osgeo.org/grass7/manuals/addons/r.cpt2grass.html), [r.colors.matplotlib](https://grass.osgeo.org/grass7/manuals/addons/r.colors.matplotlib.html)\n", "\n", "\n", " [g.region](https://grass.osgeo.org/grass72/manuals/g.region.html) \n", " resolution and extent management \n", " [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html), GUI\n", "\n", "Modules and their descriptions with examples can be found the documentation. The documentation is included in the local installation and is also available online.\n", "\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Manual_pages_online_keywords_7.2.png)\n", "\t\t\t\n", "List of keywords (tags) in the online documentation\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:R_slope_aspect_3_manual.png)\n", "\t\t\t\n", "Manual page for a module is available also from the module dialog\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "## Basic introduction to Python interface\n", "The simplest way to execute the Python code which uses GRASS GIS packages is to use Simple Python Editor integrated in GRASS GIS (accessible from the toolbar or the Python tab in the Layer Manager). Another option is to use your favorite text editor and then run the script in GRASS GIS using the main menu File -> Launch script.\n", "\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Simple_python_editor_v_buffer.png)\n", "\t\t\t\n", "Simple Python Editor integrated in GRASS GIS\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:GRASS_GUI_Python_shell.png)\n", "\t\t\t\n", "Python tab with an interactive Python shell\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "We will use the Simple Python Editor to run the commands. You can open it from the Python tab.\n", "When you open Simple Python Editor, you find a short code snippet.\n", "It starts with importing GRASS GIS Python Scripting Library:\n", "\n", "```\n", "import grass.script as gscript\n", "```\n", "In the main function we call [g.region](https://grass.osgeo.org/grass72/manuals/g.region.html) to see the current computational region settings:\n", "\n", "```\n", "gscript.run_command(g.region, flags=p)\n", "```\n", "Note that the syntax is similar to command line syntax (`g.region -p`), only the flag is specified in a parameter. Now we can run the script by pressing the Run button in the toolbar. In Layer Manager we get the output of g.region.\n", "In this example, we set the computational extent and resolution to the raster layer ortho which would be done using `g.region raster=ortho` in the command line.\n", "To use the `run_command` to set the computational region, replace the previous [g.region](https://grass.osgeo.org/grass72/manuals/g.region.html) command with the following line:\n", "\n", "```\n", "gscript.run_command(g.region, raster=ortho)\n", "```\n", "The GRASS GIS Python Scripting Library provides functions to call GRASS modules within Python scripts as subprocesses. All functions are in a package called grass and the most common functions are in grass.script package which is often imported `import grass.script as gscript`. The most often used functions include:\n", "\n", "* [script.core.run_command()](https://grass.osgeo.org/grass71/manuals/libpython/script.html#script.core.run_command): used with modules which output raster or vector data and when text output is not expected\n", "* [script.core.read_command()](https://grass.osgeo.org/grass71/manuals/libpython/script.html#script.core.read_command): used when we are interested in text output which is returned as Python string\n", "* [script.core.parse_command()](https://grass.osgeo.org/grass71/manuals/libpython/script.html#script.core.parse_command): used with modules producing text output as key=value pair which is automatically parsed into a Python dictionary\n", "* [script.core.write_command()](https://grass.osgeo.org/grass71/manuals/libpython/script.html#script.core.write_command): for modules expecting text input from either standard input or file\n", "Here we use `parse_command` to obtain the statistics as a Python dictionary\n", "\n", "```\n", "region = gscript.parse_command(g.region, flags=g)\n", "print region[ewres], region[nsres]\n", "```\n", "The results printed are the raster resolutions in E-W and N-S directions.\n", "In the above examples, we were calling the [g.region](https://grass.osgeo.org/grass72/manuals/g.region.html) module. Typically, the scripts (and GRASS GIS modules) don't change the computational region and often they don't need to even read it. The computational region can be defined before running the script so that the script can be used with different computational region settings.\n", "The library also provides several convenient wrapper functions for often called modules, for example [script.raster.raster_info()](https://grass.osgeo.org/grass71/manuals/libpython/script.html#script.raster.raster_info) (wrapper for [r.info](https://grass.osgeo.org/grass72/manuals/r.info.html)), [script.core.list_grouped()](https://grass.osgeo.org/grass71/manuals/libpython/script.html#script.core.list_grouped) (one of the wrappers for [g.list](https://grass.osgeo.org/grass72/manuals/g.list.html)), and [script.core.region()](https://grass.osgeo.org/grass71/manuals/libpython/script.html#script.core.region) (wrapper for [g.region](https://grass.osgeo.org/grass72/manuals/g.region.html)).\n", "When we want to run the script again, we need to either remove the created data beforehand using [g.remove](https://grass.osgeo.org/grass72/manuals/g.remove.html) or we need to tell GRASS GIS to overwrite the existing data. This can be done adding `overwrite=True` as an additional argument to the function call for each module or we can do it globally using `os.environ['GRASS_OVERWRITE'] = '1'` (requires `import os`).\n", "Finally, you may have noticed the first line of the script which says `#!/usr/bin/env python`. This is what Linux, Mac OS, and similar systems use to determine which interpreter to use. If you get something line `[Errno 8] Exec format error`, this line is probably incorrect or missing.\n", "\n", "## Decide if to use GUI, command line, Python, or online Jupyter Notebook\n", "* GUI: can be combined anytime with command line\n", "* command line: can be any time combined with the GUI, most of the following instructions will be for command line (but can be easily transfered to GUI or Python), both system command line and Console tab in GUI will work well\n", "* Python: you will need to change the syntax to Python\n", "* Jupyter Notebook online: link will be distributed by the instructor\n", "* Jupyter Notebook locally: this is recommended only if you are using OSGeoLive or if you are on Linux and are familiar with Jupyter\n", "## Binning of the point cloud\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Binning_and_decimation_workflow_schema_for_point_clouds.png)\n", "\t\t\t\n", "Point cloud is either binned into a raster (e.g. [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html)) and then analyzed as raster or optionally decimated (e.g. using [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html)), converted to vector (still using [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html)) and then interpolated (e.g. [v.surf.rst](https://grass.osgeo.org/grass72/manuals/v.surf.rst.html)) into raster or analyzed as a vector (e.g. [v.vect.stats](https://grass.osgeo.org/grass72/manuals/v.vect.stats.html)). Data are in gray, processes in yellow.\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Binning_count_explanation.png)\n", "\t\t\t\n", "In basic case, binning of points into a 2D raster consists of counting the number of points falling into each cell. The resulting cell value is then count of points in that cell.\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Binning_mean_explanation.png)\n", "\t\t\t\n", "In general binning involves also values associated with the points and computes statistics on these values. Here the mean of Z coordinates of all the point in each cell is computed and stored in the raster. The cells without any points are NULL (NoData) shown in white here.\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "Fastest way to analyze basic properties of a point cloud is to use binning and create a raster map. We will now use [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) to create point count (point density) raster. At this point, we don't know the spatial extent of the point cloud, so we cannot set computation region ahead, but we can tell the [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) module to determine the extent first and use it for the output using the `-e` flag. We are using a coarse resolution to maximize the speed of the process. Additionally, the `-r` flag sets the computational region to match the output." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"count_10\", method=\"n\", resolution=\"10\", flags='en')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can see the distribution pattern, but let's examine also the numbers (in GUI using right click on the layer in Layer Manager and then Metadata or using [r.info](https://grass.osgeo.org/grass72/manuals/r.info.html) directly):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('r.info', map=\"count_10\", flags='g')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can do a quick check of some of the values using query tool in the Map Display.\n", "Since there is a lot of points per one cell, we can use a finer resolution:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"count_1\", method=\"n\", resolution=\"1\", flags='en')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Look at the distribution of the values using a histogram. Histogram is accessible from the context menu of a layer in Layer Manager, from Map Display toolbar Analyze map button or using the [d.histogram](https://grass.osgeo.org/grass72/manuals/d.histogram.html) module (`d.histogram map=count_1`).\n", "\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Lidar_point_density_with_coarse_resolution_and_ortho_in_the_background.png)\n", "\t\t\t\n", "Lidar point count per cell using [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) with 10 m cells (`resolution=10`) and `-e` flag\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:GRASS_GIS_Histogramming_Tool_d.histogram_-_count_of_point.png)\n", "\t\t\t\n", "GRASS GIS Histogramming Tool (based on [d.histogram](https://grass.osgeo.org/grass72/manuals/d.histogram.html))\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:GRASS_GIS_Histogramming_Tool_wxPython_-_count_of_point.png)\n", "\t\t\t\n", "GRASS GIS Histogramming Tool (based on wxPython)\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "Now we have appropriate number of points in most of the cells and there are no density artifacts around the edges, so we can use this raster as the base for extent and resolution we will use from now on." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', raster=\"count_1\", flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, this region has a lot of cells and some of them would be empty, esp. when we start filtering the point cloud. To account for that, we can modify just the resolution of the computational region:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', res=\"3\", flags='apg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use binning to obtain the digital surface model (DSM) (resolution and extent are taken from the computational region):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"binned_dsm\", method=\"max\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To understand what the map shows, compute statistics using [r.report](https://grass.osgeo.org/grass72/manuals/r.report.html) (Report raster statistics in the layer context menu):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print gs.read_command('r.report', map=\"binned_dsm\", units=\"c\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shows that there is an outlier (at around 600 m). Change the color table to histogram equalized (`-e` flag) to see the contrast in the rest of the map (using the `viridis` color table):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"binned_dsm\", color=\"elevation\", flags='e')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check the outliers also using minimum (which is `method=min`, we already used `method=max` for DSM):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"minimum\", method=\"min\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare this with range and mean of the values in the ground raster and decide what is the permissible range." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print gs.read_command('r.report', map=\"minimum\", units=\"c\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, to see the data, we can use histogram equalized color table:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"minimum\", color=\"elevation\", flags='e')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[](/wiki/File:R_in_lidar_explanation_of_zrange_option.png) [](/wiki/File:R_in_lidar_explanation_of_zrange_option.png)The `zrange` parameter filters out points which don't fall into the specified range of Z values\n", "Now when we know the outliers and the expected range of values (60-200 m seems to be a safe and but broad enough range), use the `zrange` parameter to filter out any possible outliers (don't be confused with the `intensity_range` parameter) when computing a new DSM:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"binned_dsm_limited\", method=\"max\", zrange=\"60,200\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpolation\n", "Now we will interpolate a digital surface model (DSM) and for that we can increase resolution to obtain as much detail as possible (we could use 0.5 m, i.e. `res=0.5`, for high detail or 2 m, i.e. `res=2 -a`, for speed):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', raster=\"count_1\", flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before interpolating, let's confirm that the spatial distribution of the points allows us to safely interpolate. We need to use the same filters as we will use for the DSM itself:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"count_dsm\", method=\"n\", return_filter=\"first\", zrange=\"60,200\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First check the numbers:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print gs.read_command('r.report', map=\"count_dsm\", units=\"h,c,p\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then check the spatial distribution. For that we will need to use histogram equalized color table (legend can be limited just to a specific range of values: `d.legend raster=count_interpolation range=0,5`)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"count_dsm\", color=\"viridis\", flags='e')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we import LAS file as vector points, we keep only first return points and limit the import vertically to avoid using outliers found in the previous step. Before running it, uncheck Add created map(s) into layer tree in the [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html) dialog if you are using GUI." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"first_returns\", return_filter=\"first\", zrange=\"60,200\", flags='bt')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* \n", "\t\t\t[](/wiki/File:V.in.lidar_dialog_do_not_add_into_layer_tree.png)\n", "\t\t\t\n", "Unchecked Add created map(s) into layer tree in v.in.lidar\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Lidar_point_density_with_fine_resolution_showing_swath_overlap.png)\n", "\t\t\t\n", "Point count at fine resolution with histogram equalized viridis color table and legend with a limited range\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "Then we interpolate:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.surf.rst', input=\"first_returns\", elevation=\"dsm\", tension=\"25\", smooth=\"1\", npmin=\"80\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we could visualize the resulting DSM by computing shaded relief with [r.relief](https://grass.osgeo.org/grass72/manuals/r.relief.html) or by using 3D view (see the following sections).\n", "\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Map_Display_with_DSM_and_legend_with_histogram.png)\n", "\t\t\t\n", "DSM with legend and histogram.\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "## Terrain analysis\n", "Set the computational region based on an existing raster map:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', raster=\"count_1\", flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check if the point density is sufficient (ground is class number 2, resolution and extent are taken from the computational region):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"count_ground\", method=\"n\", class_filter=\"2\", zrange=\"60,200\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import points (the -t flag disables creation of attribute table and the -b flag disables building of topology; uncheck Add created map(s) into layer tree):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"points_ground\", class_filter=\"2\", zrange=\"60,200\", flags='bt')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* \n", "\t\t\t[](/wiki/File:Counting_ground_points_per_cell_with_r.in.lidar.png)\n", "\t\t\t\n", "Map Display with legend created using [d.legend](https://grass.osgeo.org/grass72/manuals/d.legend.html), [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) dialog, ground point density pattern in the background\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:V.in.lidar_dialog_do_not_add_into_layer_tree.png)\n", "\t\t\t\n", "Unchecked Add created map(s) into layer tree in v.in.lidar\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "The interpolation will take some time, but we can set the computational region to a smaller area to save some time before we determine what are the best parameters for interpolation. This can be done in the GUI from the Map Display toolbar or here just quickly using the prepared coordinates:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', n=\"223751\", s=\"223418\", w=\"639542\", e=\"639899\", flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* \n", "\t\t\t[](/wiki/File:WxGUI_set_region.png)\n", "\t\t\t\n", "Simple ways to set computational region extent from GUI. On the left, set region to match raster map. On the right, select the highlighted option and then set region by drawing rectangle. Then zoom to computational region.\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Show_computational_extent_in_Map_Display.png)\n", "\t\t\t\n", "Show the current computational region extent in the Map Display\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "Now we interpolate the ground surface using regularized spline with tension (implemented in [v.surf.rst](https://grass.osgeo.org/grass72/manuals/v.surf.rst.html)) and at the same time we also derive slope, aspect and curvatures (the following code is one long line):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.surf.rst', input=\"points_ground\", tension=\"25\", smooth=\"1\", npmin=\"100\", elevation=\"terrain\", slope=\"slope\", aspect=\"aspect\", pcurvature=\"profile_curvature\", tcurvature=\"tangential_curvature\", mcurvatur=\"mean_curvature\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we examine the results, especially the curvatures curvatures show a pattern which may be caused by some problems with the point cloud collection. We decrease the tension which will cause the surface to hold less to the points and we increase the smoothing. Since the raster map called range already exists, use the overwrite flag, i.e. `--overwrite` in the command line or checkbox in GUI to replace the existing raster by the new one. We use the `--overwrite` flag shortened to just `--o`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.surf.rst', input=\"points_ground\", tension=\"20\", smooth=\"5\", npmin=\"100\", elevation=\"terrain\", slope=\"slope\", aspect=\"aspect\", pcurvature=\"profile_curvature\", tcurvature=\"tangential_curvature\", mcurvatur=\"mean_curvature\", overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we are satisfied with the result, we get back to our desired raster extent:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', raster=\"count_1\", flags='pg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally interpolate the ground surface in the full extent:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.surf.rst', input=\"points_ground\", tension=\"20\", smooth=\"5\", npmin=\"100\", elevation=\"terrain\", slope=\"slope\", aspect=\"aspect\", pcurvature=\"profile_curvature\", tcurvature=\"tangential_curvature\", mcurvatur=\"mean_curvature\", overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute shaded relief:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.relief', input=\"terrain\", output=\"relief\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now combine the shaded relief raster and the elevation raster. This can be done in several ways. In GUI by changing opacity of one of them. A better result can be obtained using the [r.shade](https://grass.osgeo.org/grass72/manuals/r.shade.html) module which combines the two rasters and creates a new one. Finally, this can be done on the fly without creating any raster using the [d.shade](https://grass.osgeo.org/grass72/manuals/d.shade.html) module. The module can be used from GUI through toolbar or using the Console tab:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('d.shade', shade=\"relief\", color=\"terrain\")\n", "Image(filename=\"map.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, instead of using [r.relief](https://grass.osgeo.org/grass72/manuals/r.relief.html), we will use [r.skyview](https://grass.osgeo.org/grass7/manuals/addons/r.skyview.html):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.skyview', input=\"terrain\", output=\"skyview\", ndir=\"8\", colorized_output=\"terrain_skyview\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Combine the terrain and skyview also on the fly:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('d.shade', shade=\"skyview\", color=\"terrain\")\n", "Image(filename=\"map.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Analytical visualization based on shaded relief:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.shaded.pca', input=\"terrain\", output=\"pca_shade\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Local relief model:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.local.relief', input=\"terrain\", output=\"lrm\", shaded_output=\"shaded_lrm\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we use automatic detection of landforms (using 50 m search window):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.geomorphon', elevation=\"terrain\", forms=\"landforms\", search=\"50\", flags='m')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* \n", "\t\t\t[](/wiki/File:Different_terrain_analyses_and_visualizations_in_multiple_Map_Displays.png)\n", "\t\t\t\n", "Different terrain analyses and visualizations in multiple Map Displays\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "## Vegetation analysis\n", "[](/wiki/File:R_in_lidar_explanation_of_zrange_option.png) [](/wiki/File:R_in_lidar_explanation_of_zrange_option.png)The `zrange` parameter filters out points which don't fall into the specified range of Z values\n", "Although for many vegetation-related application, coarser resolution is more appropriate because more points are needed for the statistics, we will use just:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('g.region', raster=\"count_1\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use all the points (not using the classification of the points), but with the Z filter to get the range of heights in each cell:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"range\", method=\"range\", zrange=\"60,200\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use all the points, but with the Z filter to get the range of heights in each cell:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.mapcalc(\"vegetation_by_range = if(range > 2, 1, null())\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or lower res to avoid gaps (and possible smoothing during their filling).\n", "\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Change_color_table_interactively_with_menu_and_labels.png)\n", "\t\t\t\n", "Change color table for a raster map from layer context menu using Set color table interactively." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"height_above_ground\", method=\"max\", base_raster=\"terrain\", zrange=\"0,100\", flags='d')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case we `if(height_above_ground > 2, height_above_ground, null())`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.mapcalc(\"above_2m = if(height_above_ground > 2, 1, null())\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use [r.grow](https://grass.osgeo.org/grass72/manuals/r.grow.html) to extend the patches and fill in the holes:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.grow', input=\"above_2m\", output=\"vegetation_grow\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider the result to represent the vegetated areas. We clump (connect) the individual cells into patches using [r.clump](https://grass.osgeo.org/grass72/manuals/r.clump.html)):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.clump', input=\"vegetation_grow\", output=\"vegetation_clump\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some of the patches are very small. Using [r.area](https://grass.osgeo.org/grass7/manuals/addons/r.area.html) we remove all patches smaller than the given threshold:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.area', input=\"vegetation_clump\", output=\"vegetation_by_height\", lesser=\"100\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now convert these areas to vector:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.to.vect', input=\"vegetation_by_height\", output=\"vegetation_by_height\", type=\"area\", flags='s')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far we were using elevation of the points, now we will use intensity.\n", "The intensity is used by [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) when `-j` (or `-i`) flag is provided:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"intensity\", zrange=\"60,200\", flags='j')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With this high resolution, there are some cells without any points, so we use [r.fill.gaps](https://grass.osgeo.org/grass7/manuals/addons/r.fill.gaps.html) to fill these cells ([r.fill.gaps](https://grass.osgeo.org/grass7/manuals/addons/r.fill.gaps.html) also smooths the raster as part of the gap-filling process)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.fill.gaps', input=\"intensity\", output=\"intensity_filled\", uncertainty=\"uncertainty\", distance=\"3\", mode=\"wmean\", power=\"2.0\", cells=\"8\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Grey scale color table is more appropriate for intensity:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"intensity_filled\", color=\"grey\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are some areas with very high intensity, so to visually explore the other areas, we use histogram equalized color table:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"intensity_filled\", color=\"grey\", flags='e')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's use [r.geomorphon](https://grass.osgeo.org/grass7/manuals/addons/r.geomorphon.html) again, but now with DSM and different settings. The following settings shows building roof shapes:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.geomorphon', elevation=\"dsm\", forms=\"dsm_forms\", search=\"7\", skip=\"4\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Different settings, especially higher flatness threshold (`flat` parameter) shows forested areas as combination of footslopes, slopes, and shoulders. Individual trees are represented as shoulders." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.geomorphon', elevation=\"dsm\", forms=\"dsm_forms\", search=\"12\", skip=\"8\", flat=\"10\", overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lowering the skip parameter decreases generalization and brings in summits which often represents tree tops:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.geomorphon', elevation=\"dsm\", forms=\"dsm_forms\", search=\"12\", skip=\"2\", flat=\"10\", overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we extract summits using [r.mapcalc](https://grass.osgeo.org/grass72/manuals/r.mapcalc.html) raster algebra expression `if(dsm_forms==2, 1, null())`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.mapcalc(\"trees = if(dsm_forms==2, 1, null())\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bonus tasks\n", "### 3D visualization of rasters\n", "We can explore our study area in 3D view (use image on the right if clarification is needed for below steps):\n", "\n", "* Add raster dsm and uncheck or remove any other layers. Note that unchecking (or removing) any other layers is important because any layer loaded in Layers is interpreted as a surface in 3D view.\n", "* Switch to 3D view (in the right corner of Map Display).\n", "* Adjust the view (perspective, height).\n", "* In Data tab, set Fine mode resolution to 1 and set ortho as the color of the surface (the orthophoto is draped over the DSM).\n", "* Go back to View tab and explore the different view directions using the green puck.\n", "* Go to Appearance tab and change the light conditions (lower the light height, change directions).\n", "* Try also using the landforms raster as the color of the surface.\n", "* When finished, switch back to 2D view.\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:ICC_workshop_3Dview_ortho.png)\n", "\t\t\t\n", "3D visualization of DSM with orthophoto draped over: arrows pointing to the name of the surface, resolution used for its display, and raster used as its color.\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Geomorphons_in_3D_view_Appearance_tab.png)\n", "\t\t\t\n", "Landforms over the DSM in 3D: The light source position is changed in the Appearance tab.\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "### Analytical visualization of rasters in 3D\n", "We can explore our study area in 3D view (use image on the right if clarification is needed for below steps):\n", "\n", "* Add rasters terrain and dsm and uncheck or remove any other layers. (Any layer in Layers will be interpreted as surface in 3D view.)\n", "* Switch to 3D view (in the right corner of Map Display).\n", "* Adjust the view (perspective, height). Set z-exaggeration to 1.\n", "* In Data tab, in Surface, set Fine mode resolution to 1 for both rasters.\n", "* Set a different color for each surface.\n", "* For the dsm, set the position in Z direction to 1. This is a position relative the actual position of the raster. This small offset will help is see relation in between the terrain and the DSM surfaces.\n", "* Go to Analysis tab and activate the first cutting plane. Set Shading to bottom color to get the bottom color to color the space between the terrain and DSM.\n", "* Start sliding the sliders for X, Y and Rotation.\n", "* Go back to View tab in case you want to change the view.\n", "* When finished, switch back to 2D view.\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Nviz_data_tab_dsm_set_relative_positon.png)\n", "\t\t\t\n", "In the Data tab of 3D view, set the same fine mode resolution for both surfaces. Choose different colors. Set (relative) position of the DSM to 1 m (above its actual position).\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Nviz_cutting_plane_analysis_tab.png)\n", "\t\t\t\n", "In the Analysis tab of 3D view, activate cutting plane, set shading to bottom color and start sliding the sliders for X, Y and rotation.\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Nviz_cutting_plane_dem_dsm.png)\n", "\t\t\t\n", "To see longer part of the transect, stretch the window horizontally.\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "### Classify ground and non-ground points\n", "UAV point clouds usually require classification of ground points (bare earth) when we want to create ground surface.\n", "Import all the points using the [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html) command in the previous section, unless you already have them. Since the metadata about projection are wrong, we use the `-o` flag to skip the projection consistency check." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_uav_points_spm.las\", output=\"uav_density_05\", method=\"n\", resolution=\"0.5\", flags='eno')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set a small computational region to make all the computations faster (skip this to operate in the whole area):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('g.region', n=\"219415\", s=\"219370\", w=\"636981\", e=\"637039\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import the points using the [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html) module but limit the extent just to computational region (`-r` flag), do not build topology (`-b` flag), do not create attribute table (`-t` flag), and do not assign categories (ids) to points (`-c` flag). There is more points than we need for interpolating the the resolution 0.5 m, so we will decimate (thin) the point cloud during import removing 75% of the points using `preserve=4` (which uses count-based decimation which assumes spatially uniform distribution of the points)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.in.lidar', input=\"nc_uav_points_spm.las\", output=\"uav_points\", preserve=\"4\", flags='tcbro')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then use [v.lidar.mcc](https://grass.osgeo.org/grass7/manuals/addons/v.lidar.mcc.html) to classify ground and non-ground points:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.lidar.mcc', input=\"points_all\", ground=\"mcc_ground\", nonground=\"mcc_nonground\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Interpolate the ground surface:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.surf.rst', input=\"mcc_ground\", tension=\"20\", smooth=\"5\", npmin=\"100\", elevation=\"mcc_ground\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you are using the UAV data, you can extract the RGB values from the points. First, import the points again, but this time with all the attributes:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.in.lidar', input=\"nc_uav_points_spm.las\", output=\"uav_points\", preserve=\"4\", flags='ro')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then extract the RGB values (stored in attribute columns) into separate channels (rasters):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.to.rast', input=\"uav_points\", output=\"red\", use=\"attr\", attr_col=\"red\")\n", "gs.run_command('v.to.rast', input=\"uav_points\", output=\"green\", use=\"attr\", attr_col=\"green\")\n", "gs.run_command('v.to.rast', input=\"uav_points\", output=\"blue\", use=\"attr\", attr_col=\"blue\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then set the color table to grey:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"red\", color=\"grey\")\n", "gs.run_command('r.colors', map=\"green\", color=\"grey\")\n", "gs.run_command('r.colors', map=\"blue\", color=\"grey\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These channels could be used for some remote sensing tasks. Now we will just visualize them using [d.rgb](https://grass.osgeo.org/grass72/manuals/d.rgb.html):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('d.rgb', red=\"red\", green=\"green\", blue=\"blue\")\n", "Image(filename=\"map.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you are using the lidar data for this part, you can compare the new DEM with the previous dataset:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.mapcalc(\"mcc_lidar_differece = ground - mcc_ground\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the color table:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"mcc_lidar_differece\", color=\"difference\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3D visualization of point clouds\n", "Unless you have already done that, familiarize yourself with the 3D view by doing the exercise above. Then import all points but this time do not add them to the Layer Manager (there is a check box at the bottom of the module dialog)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"points_all\", zrange=\"60,200\", flags='obt')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can explore point clouds in 3D view:\n", "\n", "* Uncheck or remove any other layers. (Any layer in Layer Manager is interpreted as surface in 3D view.)\n", "* Switch to 3D view (in the right corner of Map Display).\n", "* When 3D is active, Layer Manager has a new icon for default settings of the 3D which we need to open and change.\n", "* In the settings, change settings for default point shape from sphere to X.\n", "* Then add the points vector.\n", "* You may need to adjust the size of the points if they are all merged together.\n", "* Adjust the view (perspective, height).\n", "* Go back to View tab and explore the different view directions using the green puck.\n", "* Go to Appearance tab and change the light conditions (lower the light height, change directions).\n", "* When finished, switch back to 2D view.\n", "\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:WxGUI_nviz_with_point_cloud_ground_and_non-ground_data_tab.png)\n", "\t\t\t\n", "UAV point cloud classified to group points (orange) and non-ground points cloud (green).\n", "\n", "\t\t\t\n", "\t\t\n", "\t\t\n", "* \n", "\t\t\t[](/wiki/File:Selection_409.png)\n", "\t\t\t\n", "Detail showing trees and outliers.\n", "\n", "\t\t\t\n", "\t\t\n", "\n", "### Alternative DSM creation\n", "First returns not always represent top of the canopy, but they can simply represent first hit somewhere in the canopy. To account for this we create a raster with maximum for each cell instead of importing the raw points:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.in.lidar', input=\"nc_tile_0793_016_spm.las\", output=\"maximum\", method=\"maximum\", return_filter=\"first\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The raster has many missing cells, that's why we want to interpolate, but at the same time, we also reduced the number of points by replacing all points in one cells just by this one cell. We can consider the cell to represent one point in the middle of the cell and use [r.to.vect](https://grass.osgeo.org/grass72/manuals/r.to.vect.html) to create points from these cells:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.to.vect', input=\"maximum\", output=\"points_for_dsm\", type=\"point\", flags='b')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have a decimated (thinned) point cloud, decimated using binning into raster. We interpolate these points to get the DSM:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('v.surf.rst', input=\"first_returns\", elevation=\"dsm\", tension=\"25\", smooth=\"1\", npmin=\"80\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Explore layers of vegetation in a 3D raster\n", "Set the top and bottom of the computational region to match the height of vegetation (we are interested in vegetation from 0 m to 30 m above ground)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.parse_command('g.region', b=\"0\", t=\"30\", flags='p3g')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly to 2D, we can perform the binning also in three dimensions using 3D rasters. This is implemented in a module called [r3.in.lidar](https://grass.osgeo.org/grass72/manuals/r3.in.lidar.html):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r3.in.lidar', input=\"nc_tile_0793_016_spm.las\", n=\"count\", sum=\"intensity_sum\", mean=\"intensity\", proportional_n=\"prop_count\", proportional_sum=\"prop_intensity\", base_raster=\"terrain\", flags='d')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convert the horizontal layers (slices) of 3D raster into 2D rasters:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r3.to.rast', input=\"prop_count\", output=\"slices_prop_count\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Open a second Map Display and add all the created rasters. Then set the a consistent color table to all of them." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gs.run_command('r.colors', map=\"slices_prop_count_001,slices_prop_count_002,...\", color=\"viridis\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now open the Animation Tool and add all to 2D rasters there and use slider to explore the layers.\n", "\n", "## Optimizations, troubleshooting and limitations\n", "Speed optimizations:\n", "\n", "* Rasterize early.\n", "* For many use cases, there is less cells than points. Additionally, rasters will be faster, e.g. for visualization, because raster has a natural spatial index. Finally, many algorithms simply use rasters anyway.\n", "* If you cannot rasterize, see if you can decimate (thin) the point cloud (using [v.in.lidar](https://grass.osgeo.org/grass72/manuals/v.in.lidar.html), [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html) + [r.to.vect](https://grass.osgeo.org/grass72/manuals/r.to.vect.html), [v.decimate](https://grass.osgeo.org/grass72/manuals/v.decimate.html)).\n", "* Binning (e.g. [r.in.lidar](https://grass.osgeo.org/grass72/manuals/r.in.lidar.html), [r3.in.lidar](https://grass.osgeo.org/grass72/manuals/r3.in.lidar.html)) is much faster than interpolation and can carry out part of the analysis.\n", "* Faster count-based decimation performs often the same on a given point cloud for terrain as slower grid-based decimation (Petras et al. 2016).\n", "* r.in.lidar\n", "* choose computation region extent and resolution ahead\n", "* have enough memory to avoid using percent option (high memory usage versus high I/O)\n", "* v.in.lidar\n", "* -r limit import to computation region extent\n", "* -t do not create attribute table\n", "* -b do not build topology (applicable to other modules as well)\n", "* -c store only coordinates, no categories or IDs\n", "Memory requirements:\n", "\n", "* for r.in.lidar\n", "* depend on the side of output and type of analysis\n", "* can be reduced by percent option\n", "* ERROR: G_malloc: unable to allocate ... bytes of memory means that there is not enough memory for the output raster, use percent option, coarser resolution, or smaller extent\n", "* on Linux available memory for process is RAM + SWAP partition, but when SWAP will be used, processing will be much slower\n", "* for v.in.lidar\n", "* low when not building topology (-b flag)\n", "* to build topology but keeping low memory, use export GRASS_VECTOR_LOWMEM=1 (os.environ['GRASS_VECTOR_LOWMEM'] = '1' in Python)\n", "Limits:\n", "\n", "* Vector features with topology are limited to about 2 billion features per vector map (2^31 - 1)\n", "* Points without topology are limited only by disk space and what the modules are able to process later. (Theoretically, the count limited only by the 64bit architecture which would be 16 exbibytes per file, but the actual value depends on the file system.)\n", "* For large vectors with attributes, PostgreSQL backend is recommended for attributes ([v.db.connect](https://grass.osgeo.org/grass72/manuals/v.db.connect.html)).\n", "* There is more limits for the 32bit versions for operations which require memory. (Since 2016 there is 64bit version even for MS Windows. Even the 32bit versions have large file support (LFS) and can work with files exceeding the 32bit limitations.)\n", "* Reading and writing to disk (I/O) usually limits speed. (Can be made faster for rasters since 7.2 using different compression algorithms set using the GRASS_COMPRESSOR variable.)\n", "* There is limit for number of open files set by the operating system.\n", "* The limit is often 1024. On Linux, it can be changed using ulimit.\n", "* Individual modules may have their own limitations or behaviors with large data based on the algorithm they are using. However, in general modules are made to deal with large datasets. For example,\n", "* [r.watershed](https://grass.osgeo.org/grass72/manuals/r.watershed.html) can process 90,000 x 100,000 (9 billion cells) in 77.2 hours (2.93GHz), and\n", "* [v.surf.rst](https://grass.osgeo.org/grass72/manuals/v.surf.rst.html) can process 1 million points and interpolate 202,000 cells in 13 minutes.\n", "* Read the documentation as it provides descriptions of limitations for each module and ways to deal with them.\n", "* Write to [grass-user mailing list](https://lists.osgeo.org/listinfo/grass-user) or [GRASS GIS bug tracker](http://trac.osgeo.org/grass/) in case you hit limits. For example, if you get negative number as the number of points or cells, open a ticket.\n", "## Bleeding edge\n", "[](/wiki/File:Raster3d_example_small_data.png) [](/wiki/File:Raster3d_example_small_data.png)Vegetation fragmentation computed by [r3.forestfrag](https://grass.osgeo.org/grass7/manuals/addons/r3.forestfrag.html) from lidar point cloud represented as 3D raster\n", "* Processing point clouds as 3D rasters\n", "* r3.forestfrag\n", "* r3.count.categories\n", "* r3.scatterplot\n", "* Point profiles\n", "* v.profile.points\n", "* PDAL integration\n", "* prototypes: v.in.pdal, r.in.pdal\n", "* goal: provide access to selected PDAL algorithms (filters)\n", "* r.in.kinect\n", "* Has some features from r.in.lidar, v.in.lidar, and PCL (Point Cloud Library), but the point cloud is continuously updated from a Kinect scanner using OpenKinect libfreenect2 (used in [Tangible Landscape](http://tangible-landscape.github.io/)).\n", "* Combining elevation data from different sources (e.g. global and lidar DEM or lidar and UAV DEM)\n", "* [r.patch.smooth](https://github.com/petrasovaa/r.patch.smooth): Seamless fusion of high-resolution DEMs from two different sources\n", "* [r.mblend](https://grass.osgeo.org/grass7/manuals/addons/r.mblend.html): Blending rasters of different spatial resolution\n", "* show large point clouds in Map Display (2D)\n", "* prototype: d.points\n", "* Store return and class information as category\n", "* experimental implementation in v.in.lidar\n", "* Decimations\n", "* implemented: v.in.lidar, r.in.lidar\n", "* v.decimate (with experimental features)\n", "Contributions are certainly welcome. You can discuss them at the [grass-dev mailing list](https://lists.osgeo.org/listinfo/grass-dev).\n", "\n", "\n", "## See also\n", "* [GRASS GIS at FOSS4G Boston 2017](/wiki/GRASS_GIS_at_FOSS4G_Boston_2017)\n", "* [From GRASS GIS novice to power user (workshop at FOSS4G Boston 2017)](/wiki/From_GRASS_GIS_novice_to_power_user_(workshop_at_FOSS4G_Boston_2017))\n", "* [Analytical data visualizations at ICC 2017](/wiki/Analytical_data_visualizations_at_ICC_2017)\n", "* [Unleash the power of GRASS GIS at US-IALE 2017](/wiki/Unleash_the_power_of_GRASS_GIS_at_US-IALE_2017)\n", "* [Creating animation from FUTURES output in GRASS GIS](/wiki/Creating_animation_from_FUTURES_output_in_GRASS_GIS)\n", "* [Introduction to GRASS GIS with terrain analysis examples](/wiki/Introduction_to_GRASS_GIS_with_terrain_analysis_examples)\n", "* [Lidar Analysis of Vegetation Structure](/wiki/Lidar_Analysis_of_Vegetation_Structure)\n", "* [LIDAR](/wiki/LIDAR)\n", "* [GRASS GIS and Python](/wiki/GRASS_and_Python)\n", "* [GRASS GIS and R](/wiki/R_statistics)\n", "* [GRASS Location Wizard](/wiki/GRASS_Location_Wizard)\n", "External links:\n", "\n", "* [Processing lidar and general point cloud data in GRASS GIS](http://wenzeslaus.github.io/grass-lidar-talks/) (presentation slides)\n", "* [Seamless fusion of high-resolution DEMs from multiple sources](https://petrasovaa.github.io/dem-fusion-talk) (presentation slides)\n", "* [libLAS](https://www.liblas.org)\n", "* [PDAL](http://pdal.io/)\n", "References:\n", "\n", "* Brovelli M. A., Cannata M., Longoni U.M. 2004. LIDAR Data Filtering and DTM Interpolation Within GRASS. Transactions in GIS, April 2004, vol. 8, iss. 2, pp. 155-174(20), Blackwell Publishing Ltd. ([full text at ResearchGate](https://www.researchgate.net/publication/220606097_LiDAR_data_filtering_and_DTM_interpolation_within_GRASS))\n", "* Evans, J. S., Hudak, A. T. 2007. A Multiscale Curvature Algorithm for Classifying Discrete Return LiDAR in Forested Environments. IEEE Transactions on Geoscience and Remote Sensing 45(4): 1029 - 1038. ([full text](http://www.fs.fed.us/rm/pubs_other/rmrs_2007_evans_j001.pdf))\n", "* Jasiewicz, J., Stepinski, T. 2013. Geomorphons - a pattern recognition approach to classification and mapping of landforms, Geomorphology. vol. 182, 147-156. [DOI: 10.1016/j.geomorph.2012.11.005](http://dx.doi.org/10.1016/j.geomorph.2012.11.005)\n", "* Leito, J.P., Prodanovic, D., Maksimovic, C. 2016. Improving merge methods for grid-based digital elevation models. Computers & Geosciences, Volume 88, March 2016, Pages 115-131, ISSN 0098-3004. [DOI: 10.1016/j.cageo.2016.01.001](http://doi.org/10.1016/j.cageo.2016.01.001). \n", "* Mitasova, H., Mitas, L. and Harmon, R.S., 2005, Simultaneous spline approximation and topographic analysis for lidar elevation data in open source GIS, IEEE GRSL 2 (4), 375- 379. ([full text](http://www4.ncsu.edu/~hmitaso/gmslab/papers/IEEEGRSL2005.pdf))\n", "* Petras, V., Newcomb, D. J., Mitasova, H. 2017. Generalized 3D fragmentation index derived from lidar point clouds. In: Open Geospatial Data, Software and Standards. [DOI: 10.1186/s40965-017-0021-8](http://dx.doi.org/10.1186/s40965-017-0021-8)\n", "* Petras, V., Petrasova, A., Jeziorska, J., Mitasova, H. 2016. Processing UAV and lidar point clouds in GRASS GIS. ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences XLI-B7, 945952, 2016 ([full text at ResearchGate](https://www.researchgate.net/publication/304340172_Processing_UAV_and_lidar_point_clouds_in_GRASS_GIS))\n", "* Petrasova, A., Mitasova, H., Petras, V., Jeziorska, J. 2017. Fusion of high-resolution DEMs for water flow modeling. In: Open Geospatial Data, Software and Standards. [DOI: 10.1186/s40965-017-0019-2](http://dx.doi.org/10.1186/s40965-017-0019-2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# end the GRASS session\n", "os.remove(rcfile)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" } }, "nbformat": 4, "nbformat_minor": 0 }