Workflow for generating DEM from Lidar tiles

The original dataset comprises Lidar tiles with classified ground and with XYZ ASCII format of this type :

13.83535621 42.12214528 401.940 32 2
13.83534713 42.12213870 401.870 28 2
13.83533818 42.12213222 401.670 33 2
13.83532879 42.12211645 401.950 30 1
13.83533763 42.12212258 401.630 28 2
13.83534682 42.12212933 401.830 29 2
13.83535625 42.12213609 402.250 25 1

2 processing methods are proposes. The first is more straightforward, it uses bash only. The second uses ArcGIS.

Method 1 : with bash, GDAL, LAStools

All the input tiles are compressed in zip format.

1- zip extraction :

for i in *.zip; do unzip $i -d .; done

The output is  xyz files, with classification as last field. 1 for unclassified, 2 for ground (see example above)

delete zip files with :
rm *.zip

2-extract only ground points from xyz files

awk '{if ($5==2) print $0}' input > output

for batch processing, see script “script_xyz2xyzground.sh”

3-convert xyz format to las format

las2las -i input -o output

4-interpolation in grid with points2grid

for i in *.las; do BASE=`basename $i .las`;points2grid -i $i -o OUTPUT/$BASE -r 0.00000707 --input_format las --output_format arc --mean --resolution 0.00001 --fill; done

5- convert asc format to tif

for i in *.asc; do BASE=`basename $i .asc`;gdal_translate -a_srs EPSG:4326 $i ${BASE}.tif;done

6- merge all tif tiles

gdal_merge.py -o MOSAIC/Marche_lidar_full_mosaic_4326.tif -co COMPRESS=DEFLATE -co BIGTIFF=YES -co TILED=YES -a_nodata 0 *.tif

Optional step : fill holes

for i in *.tif; do BASE=`basename $i ..tif`;gdal_fillnodata.py $i ${BASE}_filled.tif -co compress=deflate -co bigtiff=yes -co tiled=yes;done

Method 2 : with ArcGIS

Toolbox of several small scripts for processing Lidar point cloud and generating DTMs

The original datasets are point clouds with ground already classified

Steps :

1 – extract only ground from the classified point cloud
This is done with awk. See description at this site https://sigeo.cerege.fr/?p=199

ex. of command for processing a file :
awk ‘{ if ( $5!=1 ) print $0 }’ input.xyz > output.xyz
to batch it in bash :
for i in *.xyz; do awk ‘{ if ( $5!=1 ) print $0 }’ input.xyz > output.xyz; done

2- Convert ASCII XYZ file to LAS
This is done with ArcGis using the txt2las tool of LAStools, but it can also be done with libLAS or PDAL
The script in the ArcGIS toolbox is used in a batch mode

3- Convert LAS files to LAS datasets
we created a small python script based on Arcpy that we use add to a Toolbox and use it in batch mode
see “convert_las2lasd.py”

4- Convert Las datasets to TIN
we created a small python script based on Arcpy that we use add to a Toolbox and use it in batch mode
see “convert_lasd2tin.py”

5- Convert Tin to raster
we created a small python script based on Arcpy that we use add to a Toolbox and use it in batch mode
see “convert_tin2raster.py”

The scripts of steps 2 to 5 are used from a ArcGIS Toolbox