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