From point cloud to raster DEM/DSM
E. Hardin, M.O. Kurum, H. Mitasova
Preprocessing lidar point cloud text files
# Raw LIDAR data often needs to be modified or altered to make it easier to handle.
# For example, the text files may too large and cannot be opened using common
# text editors like notepad or gedit but editing files can be done quickly
# and easily using awk or sed and some examples of how to do so are given here.
# Eric Pement compiled helpful awk and sed commands
#
# View only the first line of a file, for example, to identify a header or
# to determine the delimiter
sed q /path/filename
# View the first ten lines:
sed 10q /path/filename
# Remove the first line with awk or comment it out with sed
# (for example a header that is not needed)
awk 'NR>1{print $0}' /path/inputFile > /path/outputFile
sed s/E/#E/ /path/inputFile > /path/outputFile
# Change the delimiter, for example, from | to ,
awk -F"|" '{OFS=","}{print $1,$2,$3}' /pathTo/MyInputFile > /pathTo/MyOutputFile
# Merge many tiles, each with a header, and remove the headers.
# In this example, each tile is a text file (.txt) and the headers start with
# the word Easting as is the case of data obtained from NOAA.
# Within the directory with the tiles concatenate every .txt file into a text
# file called MyRawDataWHeaders and remove every line beginning with the word
# Easting regardless of what follows:
cat *.txt > MyRawDataWHeaders
sed '/Easting/d' MyRawDataWHeaders > MyRawData WoutHeaders
Binning
# Compute raster map where each raster cell shows the elevation range within that cell
# If the values in this map are greater than the LIDAR accuracy,
# then the DEM resolution should be finer
g.region res=10 -p
r.in.xyz input=/pathTo/MyRawLidarData out=MyRangeMap fs=, method=range
# Import vector points within the current computational region (the -r flag)
v.in.ascii -ztbr input=/pathTo/MyRawLidar output=MyImportedLidar format=point \
fs=, skip=0 x=1 y=2 z=3 cat=0
# Interpolation, smoothing and computation of slope and curvature
v.surf.rst -t MyImportedLidar elev=MyElevationMap slope=MySlopeMap pcu=MyProfCurvMap \
layer=0 segm=40 npmin=120 dmin=0.3 ten=1200 sm=0.5
# Create a mask to limit the area that is being interpolated to the study region, exclude
# water surface or interpolate smaller sections to prevent v.surf.rst running out of memory
r.digit
v.digit
r.mask or g.copy rast=mymask,MASK example here
# Masking based on contours, for example, shorelines
# Extract the shoreline at mean-high water level
# copy the contour before digitizing because changes
# (including mistakes) are applied to the map permanently
r.contour MyElevation out=MyShoreline_036m level=0.36 cut=300
g.copy vect=MyShoreline_036m,MyShorelineCopy_036m
# Use v.digit to delete spurious contours, to connect the ends of the line
# in order to define an area, and finally to add a centroid inside of the area
v.digit map=MyShorelineCopy_036m bgcmd="d.rast MyMap"
# Convert the line to a boundary
is this step needed - should we save it from v.digit as area?
v.type input= MyShorelineCopy_036m output=MyVectorArea type=line,boundary
# Convert the vector area to a raster
v.to.rast input=MyVectorArea output=MyRasterArea use=val type=point,line,area
# Multiply the elevation map by this raster map (mask)
# to get an elevation map without the ocean
r.mapcalc 'MyElevationWOutOcean=MyElevation*MyRasterArea'
# For large data sets, the module v.surf.rst may require more RAM, than available.
# and the processing needs to be done by overlapping tiles that are patched together
show here definition of the tiles for given data set to illustrate the overlap
g.region WholeStudyArea
v.in.ascii -ztbr input=/pathTo/MyRawLidar output=MyImportedLidar format=point \
fs=, skip=0 x=1 y=2 z=3 cat=0
g.region bottomHalfOfArea
v.surf.rst -t MyImportedLidar elev=MyElevationMap_Bottom layer=0 segm=40 npmin=120 \
dmin=0.3 ten=1200 sm=0.5
g.region TopHalfOfArea
v.surf.rst -t MyImportedLidar elev=MyElevationMap_Top layer=0 segm=40 npmin=120 \
dmin=0.3 ten=1200 sm=0.5
r.patch input=MyElevationMap_bottom,MyElevationMap_Top output=MyElevationMap_Whole
g.region rast=coast_2001_2m -p
d.erase
d.rast coast_2001_2m
Identification and reduction of systematic error