GIS/MEA 584:
Mapping and Analysis Using UAS

Analyzing and Visualizing surfaces from multiple return lidar data

Outline:
  • analyze and import lidar point cloud for Mid Pines area
  • compute bare earth and canopy surfaces and derived parameters
  • employ visualization techniques to highlight subtle terrain features
Data: Tools:
  • Your GRASS GIS 7 installation should also include libLAS library which is used by GRASS modules v.in.lidar and r.in.lidar (standalone GRASS GIS for MS Windows, OSGeo4W and Ubuntu packages contain libLAS).
  • Note: if v.in.lidar or r.in.lidar are not accessible, use r.in.xyz and v.in.ascii with text files.
  • libLAS installation should include command line tools lasinfo and las2txt.
  • Note that for merging las tiles (not necessary for this assignment) you need LASlib's lasmerge tool (part of LAStools).
  • Web-based plas.io (might not work with older browsers, most functionality is available in Chrome)

Quick data exploration

First view the LAS file in your web browser using plas.io. Then launch GRASS (all following commands should be pasted in the GUI command line). 
grass70

Change current working directory to the directory where you downloaded the LAS file using cd command and path or in case you work in command line in GRASS GUI just type cd and press enter and select directory using a dialog. 

cd ~/Downloads

Let's look at the metadata, particularly at the classes,
we will further work work with the classes 1 and 2.

lasinfo mid_pines_spm_2013.las

  Point Classifications
---------------------------------------------------------
	1340658 Unclassified (1) 
	2580704 Ground (2) 
	66 Low Point (noise) (7) 
	1960603 Reserved for ASPRS Definition (11)
Class 11 - for the explanation of Categories in this class refer to p. 10 in ASPRS LAS specification file

Get the geographic extent of the point cloud and then set it: 

r.in.lidar -sgo input=mid_pines_spm_2013.las
g.region n=220218 s=218694 e=637795 w=636271 -p

Density of points

First we set the resolution: 
g.region res=1 -pa

Compute the density of all points: 

r.in.lidar -o input=mid_pines_spm_2013.las output=mid_pines_dens_all method=n
d.rast mid_pines_dens_all
d.legend -f -s -d rast=mid_pines_dens_all at=5,50,7,10

Compute the density of ground points and set histogram equalized color table for both densities so that we can compare them: 

r.in.lidar -o input=mid_pines_spm_2013.las output=mid_pines_dens_ground class_filter=2 method=n
r.colors -e map=mid_pines_dens_all,mid_pines_dens_ground color=bcyr
d.rast mid_pines_dens_ground
d.legend -f -s -d rast=mid_pines_dens_ground at=5,50,7,10
Later, we can will examine the density also using v.outlier and imported vector points.

Raster binning and classification

Compute different surfaces by binning. Explore what different returns and classes show.
Create a raster map of classes: 1 - ground, 2 - vegetation, 3 - buildings.

Create DSM: 

r.in.lidar -o input=mid_pines_spm_2013.las output=mid_pines_all_max method=max resolution=3 class_filter=1,2

Create surface based on last return (represents canopy): 

r.in.lidar -o input=mid_pines_spm_2013.las output=mid_pines_last_mean resolution=3 method=mean return_filter=last class_filter=1,2

Create surface representing ground based on already classified points: 

r.in.lidar -o input=mid_pines_spm_2013.las output=mid_pines_ground_mean resolution=3 class_filter=2

Combine surfaces and create classes: 1 - ground, 2 - vegetation, 3 - buildings: 

r.mapcalc "classes = if( ! isnull(mid_pines_last_mean), 2, if( !isnull(mid_pines_ground_mean), 1, if( !isnull(mid_pines_all_max),3, null())))"

Set colors, copy and paste the following color rules into "or enter values directly" text field located in Define tab (option rules): 

r.colors map=classes 

1 220:220:180
2 0:180:0
3 150:0:0

High resolution DEM/DSM

Import as points (without attribute table, not necessary to build points): 
v.in.lidar -t -o input=mid_pines_spm_2013.las output=mid_pines_ground class_filter=2
v.in.lidar -t -o input=mid_pines_spm_2013.las output=mid_pines_surface class_filter=1,2 return_filter=first

Remove the vector layers from your Layer Manager if they were added (you can disable automatic adding of layers in the dialog for the module).

Check the point overall point density using v.outlier: 

v.outlier -e input=mid_pines_ground output=dummy outlier=dummy

Estimated point density: 1.111
Estimated mean distance between points: 0.9487

Interpolate with resolution 0.3 meter, also create slope and profile curvature map.
Set higher npmin to reduce artifacts from segmentation visible on slope and curvature maps (will be much slower!): 

g.region res=0.3
v.surf.rst input=mid_pines_ground elevation=mid_pines_ground_elev slope=mid_pines_ground_slope pcurv=mid_pines_ground_profcurv npmin=80 tension=20 smooth=1
v.surf.rst input=mid_pines_surface elevation=mid_pines_surface_elev slope=mid_pines_surface_slope pcurv=mid_pines_surface_profcurv npmin=80 tension=20 smooth=1
Visualize DEM and DSM in 3D view, use cross-sections.

Leave just raster map 'mid_pines_surface_elev' in Layer Manager, hide legend. Go to 3D view, set surface resolution 1 and color to map 'classes'.

Compute the difference of lidar DEM and GCP heights

First, create a copy of the map: 
g.copy vector=GCP_12_degrees,GCP_12_differences
Then, add the columns: 
v.db.addcolumn map=GCP_12_differences columns="dem_height DOUBLE, height_difference DOUBLE"
Query (sample) the raster elevation at the locations of GCPs: 
v.what.rast -i map=GCP_12_differences raster=mid_pines_ground_elev column=dem_height
Compute the difference and save it in the attribute column: 
v.db.update map=GCP_12_differences column=height_difference query_column="ASL - dem_height"
View the results: 
v.db.select map=GCP_12_differences columns=ASL,dem_height,height_difference
Compute extended univariate statistics 
v.univar -e GCP_12_differences column=height_difference
Alternatively use Attribute Table Manager to view the results.

Visualize point density in 3D (optional)

Convert LAS file into text file: 
las2txt -i mid_pines_spm_2013.las --parse xyzc -o mid_pines.txt

Set smaller region for creating the 3D raster: 

g.region n=219502 s=219348 w=637070 e=637276 b=110 t=135 res=2 res3=2 tbres=0.5 -p3

Use binning to create point density 3D raster: 

r3.in.xyz input=mid_pines.txt output=mid_pines_dens method=n separator=comma
You can also run r3.in.lidar, directly on the LAS file, but you have to use latest development version (7.3) of GRASS GIS. 
r3.in.lidar input=mid_pines_spm_2013.las n=mid_pines_dens

Set custom color table with r3.colors: 

0% 255:255:100
5% green
100% red

Visualize point density in 3D view using slices and isosurfaces:

  1. Leave only 3D raster mid_pines_dens in Layer Manager.
  2. Set computational region based on the 3D raster and zoom to it.
  3. Switch to 3D.
  4. On View page set Height to 50 and Z-exag to 2. Nothing is visible yet.
  5. On Data page, check Draw wire box, set resolution to 1.
  6. Add isosurface and set isosurface value 1 and press Enter.
  7. Check toggle normal direction and set different color of isosurface.
  8. Experiment with different isosurface levels (press enter to confirm the new value).
  9. Remove isosurface, change draw mode to slices, add slice and set its axis to Z.
  10. Manipulate with the slice.

3D raster profiles of 3D raster of proportional density

Note: latest version of GRASS GIS is needed (7.3).

Module r3.in.lidar can compute proportional count - count per 3D cell relative to the count per vertical column. 

r3.in.lidar input=mid_pines_spm_2013.las n=mid_pines_dens proportional_n=mid_pines_dens_proportional
We will visualize it using vertical 3D raster profiles. Download extension r3.profile: 
g.extension r3.profile
Compute profile between 2 points. 
r3.profile input=mid_pines_dens_proportional raster_output=mid_pines_dens_proportional_profile coordinates=637127,219478,637219,219469
Add raster mid_pines_dens_proportional_profile and zoom to it (lower left corner is set to 0, 0).