• Home
• methods
• data
• GIS workflow
• Outputs
• references

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 

# IN THE FOLLOWING EXAMPLES, YOU WILL NEED TO CD TO THE DIRECTORY WHERE YOUR DATA IS
# OR CHANGE THE PATHS OF THE INPUT 
# Compute raster maps where each raster cell shows the elevation range within that cell
# Compute univariate statistics on the maps
# If the average range is greater than the LIDAR accuracy, 
# then the DEM resolution should be finer
resolutions=( '0.5' '1' '10' )
method=range
for res in ${resolutions[@]}
do
    g.region res=${res} n=250777.5 s=250623.5 w=913677.5 e=913882.5
    r.in.xyz input=JR_19990909 out=JR_rinxyz_${method}_${res}m \
	fs="," method=${method} --o
    r.univar JR_rinxyz_${method}_${res}m
done
# The average ranges from the output of the r.univar command should be:
# 0.0288m, 0.120m, and 2.613m for the three resolutions

# Find extents of data using r.in.xyz and -s flag
r.in.xyz -s input=JR_20080327 output=extentScan method=mean fs="," x=1 y=2

# Import vector points within the current computational region (the -r flag)
g.region res=0.5 n=252949 s=249152 w=911642 e=914771 
for date in 19961016 19971002 19980907 19990909 19990918 19991104 2001 20040925 20051126 20080327
do
    v.in.ascii -ztbr input=JR_${date} output=JR_${date} format=point \
	fs=, skip=0 x=1 y=2 z=3 cat=0 
done

# Interpolation, smoothing and computation of slope and curvature
dates=( 19961016 19971002 19980907 19990909 19990918 19991104 2001 20040925 20051126 20080327 )
smoothing=( 1200 1200 500 1000 1000 1000 1500 1500 2000 2000 )
i=0
while [ $i -lt ${#dates[@]} ]
do
    v.surf.rst -t JR_${dates[i]} elev=JR_${dates[i]}_05mrst slope=JR_${dates[i]}_05mrst_slp \
	pcu=JR_${dates[i]}_05mrst_pcu layer=0 segm=40 npmin=120 dmin=0.3 ten=${smoothing[i]} sm=0.5
    i=$(( $i+1 ))
done

# 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
north=252950
south=249152
increments=3
step=$(( ($north-$south)/$increments ))
s=$south
n=$(( $s+$step ))
dates=( 19961016 19971002 19980907 19990909 19990918 19991104 2001 20040925 20051126 20080327 )
smoothing=( 1200 1200 500 1000 1000 1000 1500 1500 2000 2000 )
for inc in `seq 1 1 $increments`
do
  g.region n=$n s=$s
  i=0
  while [ $i -lt ${#dates[@]} ]
  do
    v.surf.rst -t JR_${dates[i]} elev=JR_${dates[i]}_05mrst_tile${inc} slope=JR_${dates[i]}_05mrst_slp_tile${inc} \
       	    pcu=JR_${dates[i]}_05mrst_pcu_tile${inc} layer=0 segm=40 npmin=120 dmin=0.3 ten=${smoothing[i]} sm=0.5
    i=$(( $i+1 ))
  done
done
for date in ${dates[@]}
do
  r.patch input=JR_${dates[i]}_05mrst_tile1,JR_${dates[i]}_05mrst_tile2,JR_${dates[i]}_05mrst_tile3 \
	output=JR_${dates[i]}_05mrst
done

g.region rast=coast_2001_2m -p
d.erase
d.rast coast_2001_2m