Processing lidar and UAV point clouds in GRASS GIS (workshop at FOSS4G Boston 2017): Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(→‎Modules: https://grasswiki.osgeo.org/wiki/Unleash_the_power_of_GRASS_GIS_at_US-IALE_2017)
(→‎Modules: images and updated table)
Line 108: Line 108:
Image:Import_raster_7.2.1.png|Import raster data
Image:Import_raster_7.2.1.png|Import raster data
</gallery></center>
</gallery></center>
=== Computational region ===
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.
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:
* defined by region extent and raster resolution
* applies to all raster operations
* persists between GRASS sessions, can be different for different mapsets
* 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
* run <code>g.region -p</code> or in menu ''Settings'' - ''Region'' - ''Display region'' to see current region settings
<center><gallery perrow=3 widths=380 heights=190>
Image:Computational_region_two_rasters.png|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.)
Image:WxGUI set region.png|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.
Image:Wxgui_computational_region_set_from_raster.png|Set computational region (extent and resolution) to match a raster (''Layers'' tab in the ''Layer Manager'')
</gallery></center>
The numeric values of computational region can be checked using:
g.region -p
After executing the command you will get something like this:
north:      220750
south:      220000
west:      638300
east:      639000
nsres:      1
ewres:      1
rows:      750
cols:      700
cells:      525000
Computational region can be set using a raster map:
g.region raster=ortho -p
Resolution can be set separately using the <code>res</code> parameter of the {{cmd|g.region}} module. The units are the units of the current location, in our case meters. This can be done in the ''Resolution'' tab of the {{cmd|g.region}} dialog or in the command line in the following way (using also the <code>-a</code> flag to print the new values):
g.region res=3 -p
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 <code>-a</code> flag is for. On the other hand, if alignment with cells of a specific raster is important, <code>align</code> parameter can be used to request the alignment to that raster (regardless of the the extent).
The following example command will use the extent from the raster named <code>ortho</code>, use resolution <code>5</code> meters, modify the extent to align it to this 5 meter resolution, and print the values of this new computational region settings:
g.region vector=ortho res=5 -a -p


=== Computational region ===
=== Computational region ===
Line 163: Line 210:


<center>
<center>
  {| class="wikitable"
{| class="wikitable"
|-
|-
! Prefix !! Function !! Example
! Prefix !! Function !! Example
Line 169: Line 216:
| r.* || raster processing || ''{{cmd|r.mapcalc}}'': map algebra
| r.* || raster processing || ''{{cmd|r.mapcalc}}'': map algebra
|-
|-
| v.* || vector processing  || ''{{cmd|v.clean}}'': topological cleaning
| v.* || vector processing  || ''{{cmd|v.surf.rst}}'': surface interpolation
|-
| i.* || imagery processing  || ''{{cmd|i.segment}}'': object recognition
|-
|-
| db.* || database management || ''{{cmd|db.select}}'': select values from table
| i.* || imagery processing  || ''{{cmd|i.segment}}'': image segmentation
|-
|-
| r3.* || 3D raster processing  || ''{{cmd|r3.stats}}'': 3D raster statistics
| r3.* || 3D raster processing  || ''{{cmd|r3.stats}}'': 3D raster statistics
Line 179: Line 224:
| t.* || temporal data processing  || ''{{cmd|t.rast.aggregate}}'': temporal aggregation
| t.* || temporal data processing  || ''{{cmd|t.rast.aggregate}}'': temporal aggregation
|-
|-
| g.* || general data management || ''{{cmd|g.rename}}'': renames map
| g.* || general data management || ''{{cmd|g.remove}}'': removes maps
|-
|-
| d.* || display ||''{{cmd|d.rast}}'': display raster map
| d.* || display ||''{{cmd|d.rast}}'': display raster map
|}</center>
|}</center>


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.
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.


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.
<center><gallery perrow=4 widths=300 heights=190>
Image:R_in_lidar_dialog.png|r.in.lidar dialog with Output tab active and highlighted module name (blue), options and flags (red) and option values (green)
Image:Example r.in.lidar command in Bash.png|Example r.in.lidar command in Bash with highlighted module name (blue), options and flags (red) and option values (green)
Image:Example r.in.lidar command in Python.png|Example r.in.lidar command in Python with highlighted module name (blue), options and flags (red) and option values (green) and import (grey)
Image:Graphical Modeler with r.in.lidar and terrain analysis.png|GRASS GIS ''Graphical Modeler''
</gallery></center>


== Basic introduction to Python interface ==
== Basic introduction to Python interface ==

Revision as of 17:54, 19 July 2017

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.

Requirements: This workshop is accessible to beginners, but some basic knowledge of lidar processing or GIS is helpful for a smooth experience.

Authors: Vaclav Petras, Anna Petrasova, and Helena Mitasova from North Carolina State University

Contributors: Robert S. Dzur and Doug Newcomb

Testers:

Preparation

Software

GRASS GIS 7.2 compiled with libLAS is needed (e.g. r.in.lidar should work).

OSGeo-Live

All needed software is included.

Ubuntu

Install GRASS GIS from packages:

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install grass

Linux

For other Linux distributions other then Ubuntu, please try to find GRASS GIS in their package managers.

MS Windows

Download the standalone GRASS GIS binaries from grass.osgeo.org.

Mac OS

Install GRASS GIS using Homebrew osgeo4mac:

brew tap osgeo/osgeo4mac
brew install grass7

Note that currently (summer 2017) r.in.lidar is often not accessible on Mac OS, use r.in.ascii in combination with libLAS or PDAL command line tools to achieve the same. Note also that the 3D view may not be accessible.

Addons

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:

g.extension r.geomorphon
g.extension r.skyview
g.extension r.local.relief
g.extension r.shaded.pca
g.extension r.area
g.extension r.terrain.texture
g.extension r.fill.gaps

For bonus tasks:

g.extension v.lidar.mcc

Data

Basic introduction to graphical user interface

GRASS GIS Spatial Database

Here we provide an overview of GRASS GIS. 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.

GRASS uses specific database terminology and structure (GRASS GIS Spatial Database) 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.

  • A GRASS GIS Spatial Database (GRASS database) consists of directory with specific Locations (projects) where data (data layers/maps) are stored.
  • 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.
  • Mapset is a collection of maps within Location, containing data related to a specific task, user or a smaller project.

Start GRASS GIS, a start-up screen should appear. Unless you already have a directory called grassdata in your Documents directory (on MS Windows) 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.

We will create a new location for our project with CRS (coordinate reference system) NC State Plane Meters with EPSG code 3358. 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. 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.

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.

Importing data

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 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.

Computational region

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.

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:

  • defined by region extent and raster resolution
  • applies to all raster operations
  • persists between GRASS sessions, can be different for different mapsets
  • 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
  • run g.region -p or in menu Settings - Region - Display region to see current region settings

The numeric values of computational region can be checked using:

g.region -p

After executing the command you will get something like this:

north:      220750
south:      220000
west:       638300
east:       639000
nsres:      1
ewres:      1
rows:       750
cols:       700
cells:      525000

Computational region can be set using a raster map:

g.region raster=ortho -p

Resolution can be set separately using the res parameter of the g.region 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 dialog or in the command line in the following way (using also the -a flag to print the new values):

g.region res=3 -p

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).

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:

g.region vector=ortho res=5 -a -p

Computational region

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.

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:

  • defined by region extent and raster resolution
  • applies to all raster operations
  • persists between GRASS sessions, can be different for different mapsets
  • 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
  • run g.region -p or in menu Settings - Region - Display region to see current region settings

The numeric values of computational region can be checked using:

g.region -p

After executing the command you will get something like this:

north:      220750
south:      220000
west:       638300
east:       639000
nsres:      1
ewres:      1
rows:       750
cols:       700
cells:      525000

Computational region can be set using a raster map:

g.region raster=ortho -p

Resolution can be set separately using the res parameter of the g.region 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 dialog or in the command line in the following way (using also the -a flag to print the new values):

g.region res=3 -p

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).

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:

g.region vector=ortho res=5 -a -p

Modules

One of the advantages of GRASS is the diversity and number of modules that let you analyze all manner of spatial and temporal. GRASS GIS has over 500 different modules in the core distribution and over 230 addon modules that can be used to prepare and analyze data layers.

GRASS functionality is available through modules (tools, functions). Modules respect the following naming conventions:

Prefix Function Example
r.* raster processing r.mapcalc: map algebra
v.* vector processing v.surf.rst: surface interpolation
i.* imagery processing i.segment: image segmentation
r3.* 3D raster processing r3.stats: 3D raster statistics
t.* temporal data processing t.rast.aggregate: temporal aggregation
g.* general data management g.remove: removes maps
d.* display d.rast: display raster map

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.

Basic introduction to Python interface

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.

The GRASS GIS 7 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:

We will use the Simple Python Editor to run the commands. You can open it from the Python tab.

When you open Simple Python Editor, you find a short code snippet. It starts with importing GRASS GIS Python Scripting Library:

import grass.script as gscript

In the main function we call g.region to see the current computational region settings:

gscript.run_command('g.region', flags='p')

Note that the syntax is similar to bash 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.

Before running any GRASS raster modules, you need to set the computational region. In this example, we set the computational extent and resolution to the raster layer ortho. Replace the previous g.region command with the following line:

gscript.run_command('g.region', raster='ortho')

The run_command() function is the most commonly used one. We will use it to compute viewshed using r.viewshed. Add the following line after g.region command:

gscript.run_command('r.viewshed', input='dsm', output='python_viewshed', coordinates=(640645, 224035), overwrite=True)

Parameter overwrite is needed if we rerun the script and rewrite the raster.

Now let's look at how big the viewshed is by using r.univar. Here we use parse_command to obtain the statistics as a Python dictionary

univar = gscript.parse_command('r.univar', map='python_viewshed', flags='g')
print univar['n']

The printed result is the number of cells of the viewshed.

GRASS GIS Python Scripting Library also provides several wrapper functions for often called modules. List of convenient wrapper functions with examples includes:

Here are two commands (to be executed in the Console) often used when working with scripts. First is setting the computational region. We can do that in a script, but it is better and more general to do it before executing the script (so that the script can be used with different computational region settings):

g.region raster=dsm

The second command is handy when we want to run the script again and again. In that case, we first need to remove the created raster maps, for example:

g.remove type=raster pattern="viewshed*"

The above command actually won't remove the maps, but it will inform you which it will remove if you provide the -f flag.

Decide if to use GUI, command line, Python, or online Jupyter Notebook

  • GUI: can be combined anytime with command line
  • 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
  • Python: you will need to change the syntax to Python
  • Jupyter Notebook online: link will be distributed by the instructor
  • Jupyter Notebook locally: this is recommended only if you are using OSGeoLive or if you are on Linux and are familiar with Jupyter

Binning of the point cloud

r.in.lidar input=.../nc_tile_0793_016_spm.las output=count_10 method=n -e resolution=10
r.info map=count_10
r.in.lidar input=.../nc_tile_0793_016_spm.las output=count_1 method=n -e resolution=1
g.region raster=density_1 -p
g.region -p res=3 -a
r.in.lidar --overwrite input=.../nc_tile_0793_016_spm.las output=ground class_filter=2
r.colors map=ground color=elevation
r.relief input=ground output=relief
d.shade shade=relief color=ground
r.in.lidar --overwrite input=.../nc_tile_0793_016_spm.las output=ground class_filter=2
r.in.lidar --overwrite input=.../nc_tile_0793_016_spm.las output=range method=range zrange=0,200

Compute density of point cloud:

 g.region  n=224034 s=223266 e=640086 w=639318 res=1
 r.in.lidar -o input=tile_0793_016_spm.las method=n output=tile_density_1m resolution=1
 r.in.lidar -o input=tile_0793_016_spm.las method=n output=tile_density_6m resolution=6

We can use binning to compute the range of elevation values and other statistics:

 r.in.lidar -o input=tile_0793_016_spm.las method=range output=tile_range_6m resolution=6
 # the map and the following statistics show there is an outlier at around 600 m
 r.report map=tile_range_6m units=c
 # change the color table to histogram equalized to see the contrast
 r.colors map=tile_range_6m color=viridis -e

We can filter data by class or return:

 r.in.lidar -o input=tile_0793_016_spm.las method=mean output=tile_ground_6m resolution=6 class_filter=2
 r.colors map=tile_ground_6m color=elevation

Interpolation

Increase resolution:

g.region -pa res=0.5

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. If using GUI, uncheck add to before running it.

 v.in.lidar -obt input=tile_0793_016_spm.las  output=tile_points return_filter=first zrange=0,300

Then we interpolate:

 g.region -a -p vector=tile_points res=2
 v.surf.rst input=tile_points elevation=tile_dsm tension=25 smooth=1 npmin=100

Now you can visualize the resulting DSM using 3D view (first uncheck all other layers and then switch to 3D) or by computing shaded relief with r.relief.

Terrain analysis

g.region raster=
v.in.lidar -obt input=.../nc_tile_0793_016_spm.las  output=points_ground class_filter=2 zrange=0,300
g.region n=223751 s=223418 w=639542 e=639899
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
g.region raster=
time 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 --o
r.geomorphon -m --overwrite elevation=terrain forms=landforms search=50 skip=0 dist=0
r.skyview input=terrain output=skyview ndir=8 colorized_output=terrain_skyview
d.shade shade=skyview color=terrain
r.shaded.pca input=terrain output=pca_shade
r.local.relief input=terrain output=lrm shaded_output=shaded_lrm
r.terrain.texture elevation=terrain thres=0 pitdensity=pit_density

Vegetation analysis

The zrange parameter filters out points which don't fall into the specified range of Z values
g.region raster= res=
r.in.lidar input=.../nc_tile_0793_016_spm.las output=range method=range zrange=0,200
r.mapcalc "height_above_ground_2m = if(height_above_ground > 2, height_above_ground, null())"
r.mapcalc "presence = if(height_above_ground > 2, 1, null())"
r.grow --overwrite input=presence output=vegetation
r.clump --overwrite input=vegetation output=veg_clump
r.area input=veg_clump output=veg_size lesser=100
r.to.vect -s --overwrite input=veg_size output=vegetation type=area
g.region res=1
r.in.lidar -e -i input=.../nc_tile_0793_016_spm.las output=intensity zrange=0,300
r.fill.gaps --overwrite input=intensity output=intensity_filled uncertainty=uncertainty distance=3 mode=wmean power=2.0 cells=8
r.colors map=intensity_filled color=grey
r.colors -e map=intensity_filled color=grey

Or lower res to avoid gaps (and possible smoothing during their filling).

r.in.lidar -e -d --overwrite input=.../nc_tile_0793_016_spm.las output=hag method=max base_raster=terrain resolution=3


r3.in.lidar -d 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
r3.to.rast input=prop_count output=slices_prop_count
r3.colors

Second Map Display and Animation Tool

Bonus tasks

3D visualization

3D visualization of DSM with orthophoto draped over

We can explore our study area in 3D view (use image on the right if clarification is needed for below steps):

  1. Add raster dsm and uncheck or remove any other layers. Any layer in Layer Manager is interpreted as surface in 3D view.
  2. Switch to 3D view (in the right corner of Map Display).
  3. Adjust the view (perspective, height).
  4. 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).
  5. Go back to View tab and explore the different view directions using the green puck.
  6. Go again to Data tab and change color to viewshed raster computed in the previous steps.
  7. Go to Appearance tab and change the light conditions (lower the light height, change directions).
  8. When finished, switch back to 2D view.

Classify ground and non-ground points

v.lidar.mcc input=points_all ground=mcc_ground nonground=mcc_nonground

Optimizations, troubleshooting and limitations

Speed optimizations:

  • rasterize early (often less cells than points, natural spatial index, many algorithms use it)
  • if you can rasterize, decimate
  • binning is much faster than interpolation
  • faster count-based decimation performs the same on a given point cloud for terrain as slower grid-based decimation (Petras et al. 2016)
  • r.in.lidar
    • choose computation region extent and resolution ahead
    • have enough memory to avoid using percent option (high memory usage versus high I/O)
  • v.in.lidar
    • -r limit import to computation region extent
    • -t do not create attribute table
    • -b do not build topology (applicable to other modules as well)
    • -c store only coordinates, no categories or IDs

Memory requirements:

  • r.in.lidar
    • depend on the side of output and type of analysis
    • can be reduced by percent option
    • ERROR: G_malloc: unable to allocate ... bytes of memory
    • on Linux available memory for process is RAM + SWAP partition
  • v.in.lidar
    • low when not building topology
    • for topology and low memory use export GRASS_VECTOR_LOWMEM=1 (os.environ['GRASS_VECTOR_LOWMEM'] = '1' in Python)

Limits:

  • vector features with topology: count limit is about 2 billion features per vector map 2^31 - 1
  • points without topology: count theoretically limited only by the 64bit architecture (may be 16 exbibytes per file but depends on the file system)
  • for large vectors with attributes PostgreSQL backend is recommended
  • more limits for 32bit versions for operations which require memory
  • read/write to disk often limits speed (faster for rasters in 7.1)
  • number of open files limitation (system limit)
    • often 1024, change using ulimit on Linux
  • r.watershed 90,000 x 100,000 (9 billion cells) in 77.2 hours (2.93GHz)
  • v.surf.rst: minutes = 1.5e-11 * points^2 (test: 13 min for 1 million points and 201880 cells)
  • read the documentation
  • write to grass-user mailing list or GRASS GIS Tracker in case you hit limits

Bleeding edge

r3.forestfrag

PDAL integration

v.in.pdal, r.in.pdal, ...

r.in.kinect

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).

See also

External links and references:

  • 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
  • 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, 945–952, 2016 (full text at ResearchGate)
  • Processing lidar and general point cloud data in GRASS GIS (slides)
  • 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