GRASS-Wiki - User contributions [en]

LIDAR

2010-01-08T18:10:57Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ removed windows issues. likely to be wrong, have spoken with Hamish about this. Could have been due to bug #809

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[MB-System|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + [[MB-System|MBsys]] + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).
* NB for {{cmd|v.outlier}}, {{cmd|v.lidar.edgedetection}} and {{cmd|v.surf.bspline}}, one spline steps equates to 1m. It is recommended as a starting point that the choice of spline step is roughly 3 or 4 times the planimetric resolution (potential grid resolution) of your data. Experiment from there to obtain better results.

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
# and note use of bicubic for DSM and bilinear for DTM here and below
v.surf.bspline -c elev_lid_dsm sie=100 sin=100 type=bicubic
v.surf.bspline -c elev_lid_dtm sie=100 sin=100 type=bilinear

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1 type=bicubic
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001 type=bilinear

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2010-01-07T10:18:38Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ very minor edit

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[MB-System|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + [[MB-System|MBsys]] + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).
* NB for {{cmd|v.outlier}}, {{cmd|v.lidar.edgedetection}} and {{cmd|v.surf.bspline}}, one spline steps equates to 1m. It is recommended as a starting point that the choice of spline step is roughly 3 or 4 times the planimetric resolution (potential grid resolution) of your data. Experiment from there to obtain better results.

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
# and note use of bicubic for DSM and bilinear for DTM here and below
v.surf.bspline -c elev_lid_dsm sie=100 sin=100 type=bicubic
v.surf.bspline -c elev_lid_dtm sie=100 sin=100 type=bilinear

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1 type=bicubic
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001 type=bilinear

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

'''Windows users (6.4.x only?):'''
* When using the DBF driver (default) it is necessary to manually create the DBF table prior to importing points with {{cmd|v.in.ascii}}:
# Manually connect the DBF table
db.connect driver=dbf database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf'

# Display the connection (bug: has to be a separate command and run or the connection does not work)
db.connect -p

* Bug?: {{cmd|v.in.ascii}} is the only module that requires the mapset to be specified with the input file like with {{cmd|v.lidar.growing}} above.

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T16:52:27Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ added to windows users section

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).
* NB for {{cmd|v.outlier}}, {{cmd|v.lidar.edgedetection}} and {{cmd|v.surf.bspline}}, one spline steps equates to 1m. It is recommended as a starting point that the choice of spline step is roughly 3 or 4 times the planimetric resolution (potential grid resolution) of your data. Experiment from there to obtain better results.

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
# and note use of bicubic for DSM and bilinear for DTM here and below
v.surf.bspline -c elev_lid_dsm sie=100 sin=100 type=bicubic
v.surf.bspline -c elev_lid_dtm sie=100 sin=100 type=bilinear

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1 type=bicubic
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001 type=bilinear

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

'''Windows users (6.4.x only?):'''
* When using the DBF driver (default) it is necessary to manually create the DBF table prior to importing points with {{cmd|v.in.ascii}}:
# Manually connect the DBF table
db.connect driver=dbf database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf

# Display the connection (bug: has to be a separate command and run or the connection does not work)
db.connect -p

* Bug?: {{cmd|v.in.ascii}} is the only module that requires the mapset to be specified with the input file like with {{cmd|v.lidar.growing}} above.

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T16:50:51Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ added to windows users section

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).
* NB for {{cmd|v.outlier}}, {{cmd|v.lidar.edgedetection}} and {{cmd|v.surf.bspline}}, one spline steps equates to 1m. It is recommended as a starting point that the choice of spline step is roughly 3 or 4 times the planimetric resolution (potential grid resolution) of your data. Experiment from there to obtain better results.

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
# and note use of bicubic for DSM and bilinear for DTM here and below
v.surf.bspline -c elev_lid_dsm sie=100 sin=100 type=bicubic
v.surf.bspline -c elev_lid_dtm sie=100 sin=100 type=bilinear

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1 type=bicubic
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001 type=bilinear

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

'''Windows users (6.4.x only?):'''
* When using the DBF driver (default) it is necessary to manually create the DBF table prior to importing points with {{cmd|v.in.ascii}}:
# Manually connect the DBF table
db.connect driver=dbf database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf

# Display the connection (has to be a separate command and run or the connection does not work)
db.connect -p

* Bug?: {{cmd|v.in.ascii}} is the only module that requires the mapset to be specified with the input file like with {{cmd|v.lidar.growing}} above.

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T15:56:31Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ spelling error of mine :)

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).
* NB for {{cmd|v.outlier}}, {{cmd|v.lidar.edgedetection}} and {{cmd|v.surf.bspline}}, one spline steps equates to 1m. It is recommended as a starting point that the choice of spline step is roughly 3 or 4 times the planimetric resolution (potential grid resolution) of your data. Experiment from there to obtain better results.

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
# and note use of bicubic for DSM and bilinear for DTM here and below
v.surf.bspline -c elev_lid_dsm sie=100 sin=100 type=bicubic
v.surf.bspline -c elev_lid_dtm sie=100 sin=100 type=bilinear

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1 type=bicubic
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001 type=bilinear

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

'''Windows users (6.4.x only?):'''
* When using the DBF driver (default) it is necessary to manually create the DBF table prior to importing points with {{cmd|v.in.ascii}}:
# Manually connect the DBF table
db.connect driver=dbf database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf

# Display the connection (has to be a separate command and run or the connection does not work)
db.connect -p

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T15:55:11Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ addition of bicubic and bilinear commands in v.surf.spline section

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).
* NB for {{cmd|v.outlier}}, {{cmd|v.lidar.edgedetection}} and {{cmd|v.surf.bspline}}, one spline steps equates to 1m. It is recommend as a starting point that the choice of spline step is roughly 3 or 4 times the planimetric resolution (potential grid resolution) of your data. Experiment from there to obtain better results.

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
# and note use of bicubic for DSM and bilinear for DTM here and below
v.surf.bspline -c elev_lid_dsm sie=100 sin=100 type=bicubic
v.surf.bspline -c elev_lid_dtm sie=100 sin=100 type=bilinear

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1 type=bicubic
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001 type=bilinear

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

'''Windows users (6.4.x only?):'''
* When using the DBF driver (default) it is necessary to manually create the DBF table prior to importing points with {{cmd|v.in.ascii}}:
# Manually connect the DBF table
db.connect driver=dbf database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf

# Display the connection (has to be a separate command and run or the connection does not work)
db.connect -p

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T15:44:37Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ added NB about spline steps = 1m

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).
* NB for {{cmd|v.outlier}}, {{cmd|v.lidar.edgedetection}} and {{cmd|v.surf.bspline}}, one spline steps equates to 1m. It is recommend as a starting point that the choice of spline step is roughly 3 or 4 times the planimetric resolution (potential grid resolution) of your data. Experiment from there to obtain better results.

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
v.surf.bspline -c elev_lid_dsm sie=100 sin=100
v.surf.bspline -c elev_lid_dtm sie=100 sin=100

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

'''Windows users (6.4.x only?):'''
* When using the DBF driver (default) it is necessary to manually create the DBF table prior to importing points with {{cmd|v.in.ascii}}:
# Manually connect the DBF table
db.connect driver=dbf database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf

# Display the connection (has to be a separate command and run or the connection does not work)
db.connect -p

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T15:36:42Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ added windows user points at bottom of section.

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB. currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
v.surf.bspline -c elev_lid_dsm sie=100 sin=100
v.surf.bspline -c elev_lid_dtm sie=100 sin=100

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

'''Windows users (6.4.x only?):'''
* When using the DBF driver (default) it is necessary to manually create the DBF table prior to importing points with {{cmd|v.in.ascii}}:
# Manually connect the DBF table
db.connect driver=dbf database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf

# Display the connection (has to be a separate command and run or the connection does not work)
db.connect -p

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T13:14:05Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ v.outlier error warning

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data (NB. currently a selection of </=4 for soe/son in {{cmd|v.outlier}} results in an error message).
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
v.surf.bspline -c elev_lid_dsm sie=100 sin=100
v.surf.bspline -c elev_lid_dtm sie=100 sin=100

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-05T13:09:45Z

⚠️John Tate: /* DEM/DSM separation the more complex way */ deleted previous edit

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data.
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
v.surf.bspline -c elev_lid_dsm sie=100 sin=100
v.surf.bspline -c elev_lid_dtm sie=100 sin=100

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]

LIDAR

2009-11-03T20:30:46Z

⚠️John Tate: added using -z flag to lidar point cloud line in general procedure for DEM/DSM complex way

<h1> LIDAR and Multi-beam Swath bathymetry data </h1>

== Modules ==

=== Import ===

* {{cmd|r.in.xyz}} - Create a raster map from an assemblage of many coordinates using univariate statistics. ([http://hamish.bowman.googlepages.com/grassfiles#xyz example])

* {{cmd|v.in.ascii}} - Import data from an ASCII file to GRASS vector format.
: ''Due to memory overhead vector point imports will be limited to a few million data points unless topology and database creation is skipped with the '''-bt''' flags''. It may also be useful to clip the import file to only accept points falling within the current region by using the '''-r''' flag. See {{cmd|g.region}} for details on specifying the region bounds.

=== Analysis ===

* {{cmd|v.outlier}} - Removes outliers from vector point data.

* {{cmd|v.lidar.edgedetection}} - Detect the object's edges from a LIDAR data set.

* {{cmd|v.lidar.growing}} - Building contour determination and Region Growing algorithm for determining the building inside.

* {{cmd|v.lidar.correction}} - Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.

=== Surface generation ===

* {{cmd|v.surf.rst}} - Spatial approximation and topographic analysis using regularized spline with tension.

* {{cmd|v.surf.idw}} - Surface interpolation from vector point data by Inverse Distance Squared Weighting.

* {{cmd|v.surf.bspline}} - Surface interpolation from vector point data by bicubic or bilineal interpolation with Tykhonov regularization.

* {{cmd|r.fillnulls}} - Fills no-data areas in raster maps using v.surf.rst splines interpolation.

== Swath Bathymetry Tools ==

''see also the [[Marine_Science#Multibeam_sonar_processing]] wiki page''

* The [http://david.p.finlayson.googlepages.com/swathwidth v.swathwidith] module by David Finlayson for planning surveys. (development page)

* An example of [http://bambi.otago.ac.nz/hamish/grass/gdal/sidescan_warp.html post-processing scanned paper sidescan swaths] using thin plate spline warping with [http://www.gdal.org/ GDAL's] "<tt>gdalwarp -tps</tt>" function. (debugging page)

* [[Marine_Science#Import_into_GRASS|GRASS integration]] with [http://www.ldeo.columbia.edu/res/pi/MB-System/ MB-System] (GPL) software for processing Multibeam and Sidescan Sonar data. GRASS + MBsys + [[GMT]] make a nice scriptable trio.

== LIDAR Tools ==

* {{cmd|r.terraflow}} - computation of flow direction, flow accumulation and other basic topographic terrain indices from massive raster digital elevation models (DEM). From the Duke University [http://terrain.cs.duke.edu/ STREAM] project.

* [http://mpa.itc.it/markus/grass61/demos/rlake/ Flood simulation] using {{cmd|r.lake}}. Includes fancy NVIZ visualization of Trento, Italy, by Markus Neteler.

* [http://www.cs.unc.edu/~isenburg/software/ LAStools] are a set of simple command line tools (including source code) for converting to/from ASCII, viewing, comparing, and compressing LIDAR data.

* [http://www.liblas.org libLAS] ASPRS LiDAR data translation tools

* See also [[LIDAR#Analysis|LIDAR Analysis]]

== Micro-tutorial for LAS data==

=== Preparation ===

==== Conversion of text files to LAS ====

[http://www.liblas.org/ libLAS] supports the following column types:
x - x coordinate
y - y coordinate
z - z coordinate
a - scan angle
i - intensity
n - number of returns for given pulse (1..n)
r - number of this return (1..r)
c - classification
u - user data (but only 1 byte)
p - point source ID
e - edge of flight line
d - direction of scan flag
t - GPS time
s - skip column

Sample text data such as:
returntime,pulse,east,north,height,intensity,stripe
549778.907200,1,673999.940,5099680.080,507.425,20.0,45105
...

can be converted to LAS format like this:
returntime - t
pulse - r
east - x
north - y
height - z
intensity - i
stripe - s

First parse (sanity check):
txt2las -parse trxyzis data.asc

Then convert:
txt2las -parse trxyzis -i data.asc -o data.las

While you may not need to do this conversion for GRASS import, the resulting files are way smaller than the uncompressed text files; additionally, they are in a defined format.

=== LAS Import ===

* Use the [http://www.liblas.org libLAS] utilities to import LAS data into GRASS.
* Data stored in text files can generally be imported directly into GRASS.

* The data used in the rest of this micro-tutorial can be found at [http://www.appliedimagery.com/Serpent%20Mound%20Model%20LAS%20Data.las appliedimagery.com].

==== Import LAS as raster DEM ====

Check bounds and SRS:
$ lasinfo "Serpent_Mound_Model_LAS_Data.las"

[...]
Min X Y Z 289020.900000 4320942.610000 166.780000
Max X Y Z 290106.020000 4323641.570000 215.480000
Spatial Reference +proj=utm +zone=17 +ellps=WGS84 +units=m

After creating a suitable UTM zone 17 location (EPSG:32617)
set the region according to the information from lasinfo at 1m resolution, rounding grid outwards to align to whole meters:
GRASS> {{cmd|g.region}} n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap

Finally, import with {{cmd|r.in.xyz}} with data piped directly from the <tt>las2txt</tt> program and set a nice equalized color table:
BASEMAP="Serpent_Mound_Model_LAS"

las2txt --stdout "${BASEMAP}_Data.las" | \
r.in.xyz in=- out=${BASEMAP}_Data fs=space method=mean

{{cmd|r.colors}} ${BASEMAP}_Data color=bcyr -e

The above example uses the default z-elevation level as the 3rd term, but by using the <tt>las2txt --parse</tt> command other fields (such as intensity) can be imported instead via "<tt>--parse xyi</tt>".

==== Import LAS as vector points ====

Region setting (establishing the grid) is not needed for vector features so we can go directly to the import step. To deal with millions of input points {{cmd|v.in.ascii}} should be run with the options to skip creation of an attribute database and building topology as these can consume large amounts of memory. Note that vector maps without topology built are somewhat limited in their ability to be processed. Most LIDAR specific modules have been adapted to not require
topology. Even so, after initial cleaning steps it is often more efficient to work with huge datasets in GRASS as raster data.

las2txt --stdout "${BASEMAP}_Data.las" | \
v.in.ascii -tbz z=3 out="${BASEMAP}_pts" fs=space

If topology was built, you can use {{cmd|d.vect}}'s -z flag to colorize by elevation value. Without topology you can still colorize, but you need to use color rules based on absolute elevations, not percentage of scale.

# display colorized points for data with built topology
d.vect map=lidar_pts size=1 -z zcolor=elevation

* After import as points you can then use the v.lidar tools to clean the data: {{cmd|v.lidar.correction}}, {{cmd|v.lidar.edgedetection}}, {{cmd|v.lidar.growing}}, and {{cmd|v.outlier}}.

* Vector points can be interpolated into raster DEMs with a number of modules, including {{cmd|v.surf.rst}}, {{cmd|v.surf.bspline}}, and {{cmd|v.surf.idw}}.

=== Clean imported raster DEM ===

(fill holes)

# convert to vector points
r.to.vect -z feature=point in=${BASEMAP}_Data out=${BASEMAP}_pt

# interpolate using a regularized spline fit
# this is very slow, but produces very high quality output
v.surf.rst layer=0 in=${BASEMAP}_pt elev=${BASEMAP}.rst

# create 5m buffer area around original data points
r.buffer in=${BASEMAP}_Data out=${BASEMAP}.5m_buff dist=5

# crop interpolated DEM to only include areas nearby actual data
r.mapcalc "${BASEMAP}.filled = \
if( isnull(${BASEMAP}.5m_buff), null(), ${BASEMAP}.rst)"

# set colors to something nice
r.colors ${BASEMAP}.filled color=bcyr -e

Depending on your needs, {{cmd|r.fillnulls}}, {{cmd|v.surf.bspline}}, {{cmd|v.surf.idw}}, {{cmd|r.surf.idw}}, {{cmd|r.surf.idw2}}, or [[GRASS_AddOns#r.surf.nnbathy|r.surf.nnbathy]] may be faster than the {{cmd|v.surf.rst}} method.

=== Visualize raster DEM in 3D ===

nviz ${BASEMAP}.filled

* Set z-exag to 2.0
* In Visualize → Raster Surfaces set the fine (final) resolution to 1, and coarse (preview) resolution to 5.
* Set the height to 500.0, the perspective to 15.0, and drag the view-puck to the North-West and reasonably zoomed in.
: You should now be able to see the serpent:

[[Image:LAS_serpent_nviz.jpg|thumb|center|400px|The Great Serpent Mound, Adams County, Ohio, USA]]

It is possible to show vector points in 3D, but millions of them may make the program run slow. Topology is required ({{cmd|v.build}}). Tick the "3D" box in the Visualize → Vector points dialog.

=== LAS Export ===

==== Export to LAS ====

This is the reverse of the import step, but using {{cmd|v.out.ascii}} or {{cmd|r.out.xyz}} with ''txt2las''.
Because the v.out.ascii module exports category number which we are not interested in, we cut it away with the UNIX ''cut'' utility.

v.out.ascii ${BASEMAP}_pts fs=space | cut -f1-3 -d' ' \
> ${BASEMAP}_export.txt

txt2las --parse xyz -i ${BASEMAP}_export.txt

== Micro-tutorial for LIDAR data analysis ==

=== Simple analysis of single return data ===

Data source: North Carolina sample data set (Lidar data from 2001, near Raleigh, North Carolina). Download from
http://www.grassbook.org/ncexternal/
-> File: BE3720079200WC20020829m.txt 3.6M (lidar bare earth points [spm])

Scan input file for spatial extent. The -g flag shows the result in a convenient copy-paste format for {{cmd|g.region}}:
r.in.xyz BE3720079200WC20020829m.txt out=dummy -s -g

We use this output to set region, predefine 2m raster cells, and polish the odd coordinates with -a:
g.region n=222504.439000 s=219456.442000 e=640081.274000 w=637033.274000 res=2 -a -p

==== QUESTION 1: Are these Lidar data sufficiently dense? ====

We compute a raster map representing number of points per cell
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn2m method=n

Look at the resulting raster map:
d.mon x0
d.font Vera
d.rast.leg lid_792_binn2m

Report point distribution in percent:
r.report lid_792_binn2m unit=p

Reduce the resolution to 6m to get at least one point per cell:
g.region res=6 -a -p
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binn6m method=n
d.rast.leg lid_792_binn6m
# ... a few holes remain but that's probably acceptable.

Compute raster maps representing mean elevation for each cell:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binmean6m \
meth=mean
d.rast.leg lid_792_binmean6m

Compute range and variation:
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binrange6m \
meth=range
d.rast.leg lid_792_binrange6m
r.in.xyz BE3720079200WC20020829m.txt out=lid_792_binvar6m \
meth=coeff_var
d.rast.leg lid_792_binvar6m

Overlay other GIS maps to map of mean elevation:
d.rast.leg lid_792_binmean6m
d.vect streets_wake
d.vect lakes type=boundary
d.vect streams

Now we continue to work in a zoomed spatial subset of the area.

We import only points in the rural area without building topology and using z-coordinate for elevation:
g.region rural_1m -p
v.in.ascii -ztbr BE3720079200WC20020829m.txt \
out=elev_lidrural_bepts z=3

Clear monitor show black-white orthophoto:
d.erase -f
d.rast ortho_2001_t792_1m

We now use already prepared Lidar vector point maps:
* elev_lid792_bepts: Rural area (in tile 792), bare earth lidar point cloud
* elev_lidrural_mrpts: Rural area multiple return lidar point cloud

Look at ground data:
d.vect elev_lidrural_bepts siz=2 col=red

Look at "surface" data:
d.vect elev_lidrural_mrpts siz=1 col=green

Visualize 3D Lidar multi-return points:
nviz elev_lid792_1m point=elev_lidrural_mrpts
# Visualize -> Vector Points -> 3D points
# -> Icon size 2.25
# --> GRASS book p. 253, fig. 3.15 (see also screenshot [http://www.grassbook.org/about_menu3rd.php here])

==== DEM/DSM separation the simple way by selection of Lidar returns ====

Find out where we have multiple returns:
d.rast ortho_2001_t792_1m
d.vect elev_lidrural_mrpts where="return=1" col=red siz=2
d.vect elev_lidrural_mrpts where="return=2" col=green siz=3
d.vect elev_lidrural_mrpts where="return=3" col=blue
d.vect elev_lidrural_mrpts where="return=4" col=yellow

DTM: extract last return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_first where="Return < 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_first

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DTM
nviz elev_lidrural_mrpts_DTM col=ortho_2001_t792_1m

DSM: extract first return(s):
v.extract elev_lidrural_mrpts out=elev_lidrural_mrpts_last where="Return > 2"
nviz elev_lid792_1m point=elev_lidrural_mrpts_last

Interpolate to map and look at it:
v.surf.rst elev_lidrural_mrpts_first layer=0 elev=elev_lidrural_mrpts_DSM
nviz elev_lidrural_mrpts_DSM col=ortho_2001_t792_1m

==== DEM/DSM separation the more complex way ====

''TODO: verify order of first and last returns in below text''

General procedure:
* Lidar point clouds (first and last return) are imported with {{cmd|v.in.ascii}} (-b flag to not build the topology, and -z flag to create 3D vector map).
* Outlier detection is done with {{cmd|v.outlier}} on both first and last return data.
* Then, with {{cmd|v.lidar.edgedetection}}, edges are detected from last return data.
* The building are generated by {{cmd|v.lidar.growing}} from detected edges.
* The resulting data are post-processed with {{cmd|v.lidar.correction}}.
* Finally, the DTM and DSM are generated with {{cmd|v.surf.bspline}} (DTM: uses the 'v.lidar.correction' output; DSM: uses last return output from outlier detection).

export GRASS_OVERWRITE=1
v.extract elev_lidrural_mrpts out=elev_lidfirst_pts \
where="return = 1"
v.extract elev_lidrural_mrpts out=elev_lidlast_pts \
where="return >= 2"
d.vect elev_lidfirst_pts col=red
d.vect elev_lidlast_pts col=green

Outlier detection and separation into two maps
# 1st return
v.outlier elev_lidfirst_pts output=elev_lidfirst_clean \
outlier=elev_lidfirst_outl
d.erase
d.vect elev_lidfirst_clean siz=2
d.vect elev_lidfirst_outl col=red

2nd return
v.outlier elev_lidlast_pts output=elev_lidlast_clean \
outlier=elev_lidlast_outl
d.erase
d.vect elev_lidlast_clean siz=2
d.vect elev_lidlast_outl col=red
# -> no outliers visible

Run an edge detection on cleaned last return:
v.lidar.edgedetection elev_lidlast_clean \
out=elev_lidlast_edges

Buildings/vegetation are generated from detected edges (bug: you may need to specify the mapset):
v.lidar.growing elev_lidlast_edges@lidar out=elev_lidlast_grow \
first=elev_lidfirst_clean

Compare:
d.vect elev_lidlast_pts col=blue
d.vect elev_lidlast_grow col=green

Correction (this is applied twice):
v.lidar.correction elev_lidlast_grow out=elev_lidlast_corr1 \
terrain=elev_lidlast_terr1
v.lidar.correction elev_lidlast_corr1 out=elev_lid_dsm \
terrain=elev_lid_dtm

DEM and DSM are generated:
# Estimation of lambda_i parameter with cross validation (watch for RMS!)
v.surf.bspline -c elev_lid_dsm sie=100 sin=100
v.surf.bspline -c elev_lid_dtm sie=100 sin=100

# From the cross-validation, we select lambda with minimal RMS error:
# generate raster surfaces at 1m resolution
v.surf.bspline elev_lid_dsm raster=lidar_dsm lambda=0.1
v.surf.bspline elev_lid_dtm raster=lidar_dtm lambda=0.0001

d.rast lidar_dsm
d.rast lidar_dtm

nviz lidar_dsm,lidar_dtm \
col=ortho_2001_t792_1m,ortho_2001_t792_1m

...with the position slider you can visually separate DSM and DEM, increase z slider (Visualize -> Raster Surface -> Position).

== Sample LIDAR data ==

=== Widely used in GRASS tutorials ===

* [http://www.grassbook.org/data_menu2nd.php Jockey's Ridge, NC, LIDAR dataset]
* [http://www.grassbook.org/data_menu3rd.php North Carolina OSGeo Edu data set] (includes multi-return LIDAR data)

=== Other ===

* United States Antarctic Resource Center: LIDAR High-resolution DEM Final DATA Downloads http://usarc.usgs.gov/lidar_dload.shtml

* LIDAR ALSM Research, Arizona State University Ative Tectonics, Research Group http://lidar.asu.edu/research.html and http://www.geongrid.org/science/lidar.html

* USGS Center for LIDAR Information Coordination and Knowledge (aka CLICK) - USGS LiDAR point cloud distribution site http://lidar.cr.usgs.gov

* Washington State Geospatial Data Archive, Mount Saint Helens - Lidar Data https://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/

* Puget Sound Lidar Consortium, public-domain high-resolution topography for western Washington http://pugetsoundlidar.ess.washington.edu/index.htm

* NOAA Topographic Change Mapping LIDAR Data Retrieval Tool (LDART) NOAA Coastal Services Center http://maps.csc.noaa.gov/TCM/

* Landmap, LIDAR Data from the Environment Agency http://www.landmap.ac.uk/lidar/lidar.html

* Northern California LIDAR data http://quake.usgs.gov/research/geology/lidar/ and http://core2.gsfc.nasa.gov/lidar/terrapoint/

* IDAHO GEOSPATIAL , Bare Earth LIDAR DEM Download - UTM http://inside.uidaho.edu/geodata/LiDAR/LiDARBareEarthDEM_DownloadUTM.htm

* EarthScope Spatial Data Explorer - A java application for querying, browsing, and acquiring data from the EarthScope Spatial Data Repository. Currently includes a number of LiDAR datasets. http://www.earthscope.org/data/lidar.php

* LIDAR Galley https://peoiewswebinfo.monmouth.army.mil/JPSD/LIDAR_Gallery.htm

* South Tyrol - Download of DTMs (Homepage in German or Italian) http://www.provinz.bz.it/raumordnung/grundkarten/utm/default_d.htm

== Links ==

* [http://liblas.org/ libLAS] - LAS 1.0/1.1 ASPRS LiDAR data translation toolset

[[Category:Documentation]]