This is an unmaintained course material, please see current material at:
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
- Classified point cloud (merged from 4 tiles) in las format
- Secref and Mid Pines point clouds
- for additional point cloud tiles:
- see the overview vector file in GRASS,
- download preselected tiles from here,
- follow the instructions for importing lidar data to GRASS GIS.
- ASPRS LAS specification file with the list of ASPRS Standard LIDAR Point Classes on p.10
- 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).
- 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)
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
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
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
v.db.addcolumn map=GCP_12_differences columns="dem_height DOUBLE, height_difference DOUBLE"
v.what.rast -i map=GCP_12_differences raster=mid_pines_ground_elev column=dem_height
v.db.update map=GCP_12_differences column=height_difference query_column="ASL - dem_height"
v.db.select map=GCP_12_differences columns=ASL,dem_height,height_difference
v.univar -e GCP_12 column=height_difference
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 (TBD):
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
Visualize in 3D view using slices and isosurfaces.