Strip Alignment and DTM generation through the estimation of the "median ground" terrain points. A methodology for DJI Zenmuse L1

:: Strip Alignment and Digital Terrain Model (DTM) generation through the estimation of the "median ground" terrain points ::

A methodology to process data from LiDAR systems with Livox MID-40, Livox Horizon or Livox AVIA laser heads https://www.livoxtech.com/, of course including the DJI Zenmuse L1: https://www.dji.com/zenmuse-l1.

This publication is based on the articles published on March 4, and April 8, 2021 by Dr. Martin Isenburg on LAStools company website https://lastools.org/, in which he shared the methodology for a small project that we developed at Samara, Costa Rica, using data from an Unmanned Laser Scanning manufactured by the company Geosun https://www.geosunlidar.com/ the GS-MID40. So I adapted the methodology to process DJI Zenmuse L1 data and I added strip alignment routine. Martin taught me the importance of sharing knowledge, so following his work philosophy, I am publishing this workflow.

Photo: Dr. Martin Isenburg and Nelson Mattie working at Samara, Costa Rica.

:: The project ::

The project is composed of a flight made with the drone DJI Matrice 300 RTK and Zenmuse L1.

Image 1. Unmanned Laser Scanning DJI Zenmuse L1.

:: Quality checks with LAStools ::

lasinfo ^

-i C:\ZenmuseL1\*.laz ^

-histo intensity 16 ^

-histo gps_time 10 ^

-histo z 1 ^

-odir C:\ZenmuseL1\01_quality -odix_info -otxt ^

-cores 12

Image 2. Misalignments between strips were detected

:: Quality checks and Strip alignment using Bayesmap Strip Alignment 2.2 software https://bayesmap.com/ ::

Sometimes, LiDAR point clouds generated with the DJI Zenmuse L1 system presented misalignment problems, horizontal and/or vertical displacements and noise.

This problem has been reported by several users who already have this system in America and I have been able to confirm it directly.

With Bayesmap, the first thing we did was reconstruct the strips, since DJI Terra generates a single LAS file for the entire block.

Image 3. Checking flight line overlap and at the survey boundary

StripAlign.exe -i C:\ZenmuseL1\*.laz -uav -cut -po C:\ZenmuseL1\*.out -odir C:\ZenmuseL1\cut

StripAlign.exe -i C:\ZenmuseL1\cut\*.laz -uav -align -attrec -po C:\ZenmuseL1\*.out -odir C:\ZenmuseL1\corr

:: We produced a new LAS data set corrected ::

Image 4. Original Z difference and Corrected Z difference

:: Metadata, looking for spurious and noise values ::

lasinfo ^

-i C:\ZenmuseL1\corr\0.new_laz\*.laz ^

-merged ^

-histo z 1 ^

-odir C:\ZenmuseL1\corr\1.lasinfo ^

-o strips_cleaned_merged_info.txt

:: The relevant excerpts of the output of the lasinfo report are shown below ::

[…]

z coordinate histogram with bin size 1.000000

?bin 168 has 29344

?bin 169 has 622895

?bin 170 has 374478

?bin 171 has 315403

?bin 172 has 870666

?bin 173 has 7524467

?bin 174 has 6491927

?bin 175 has 7336407

?bin 176 has 14351979

?bin 177 has 14493398

?bin 178 has 13790638

?bin 179 has 9952524

?bin 180 has 8485887

?bin 181 has 13530489

?bin 182 has 12534225

?bin 183 has 18339524

?bin 184 has 48655482

?bin 185 has 47683441

?bin 186 has 13800749

?bin 187 has 7712764

?bin 188 has 4976386

?bin 189 has 4949985

?bin 190 has 5009978

?bin 191 has 4993566

?bin 192 has 5121673

?bin 193 has 5245282

?bin 194 has 5183363

?bin 195 has 5253325

?bin 196 has 4940315

?bin 197 has 4551300

?bin 198 has 4136114

?bin 199 has 3715658

?bin 200 has 3373698

?bin 201 has 2868903

?bin 202 has 2433275

?bin 203 has 2029894

?bin 204 has 1681965

?bin 205 has 1330334

?bin 206 has 964791

?bin 207 has 588167

?bin 208 has 387342

?bin 209 has 260242

?bin 210 has 178457

?bin 211 has 114902

?bin 212 has 63019

?bin 213 has 26376

?bin 214 has 2152

?bin 215 has 181

?bin 270 has 4

?bin 276 has 3

[…]

:: The few points above 214 meters in elevation can be considered spurious ::

:: With lastile I created tiles of 125 meters with a buffer of 20 meters, while removing the points identified as spurious with the appropriate filters ::

lastile ^

-i C:\ZenmuseL1\corr\0.new_laz\*.laz ^

-drop_z_below 168 -drop_z_above 214 ^

-apply_file_source_ID ^

-tile_size 125 -buffer 20 ^

-odir C:\ZenmuseL1\corr\2.lastiles -olaz

Image 5. This produces 67 buffered tiles.

:: The next three lasthin runs mark a sub set of low candidate points for our lasground filtering. In every 25 cm by 25 cm, every 33 cm by 33 cm and every 50 cm by 50 cm area we reclassify the point closest to the 10th percentile as class 8. In the first call to lasthin we put all other points into class 1 ::

Percentile values (1st, 5th, 10th , 20th, 25th, 30th, 40th, 50th, 60th, 70th, 75th, 80th, 90th, 95th, 99th percentiles)

lasthin ^

-i C:\ZenmuseL1\corr\2.lastiles\*.laz ^

-set_classification 1 ^

-step 0.25 -percentile 10 20 -classify_as 8 ^

-odir C:\ZenmuseL1\corr\3.lasthin1 -olaz ^

-cores 12

lasthin ^

-i C:\ZenmuseL1\corr\3.lasthin1\*.laz ^

-step 0.3333 -percentile 10 20 -classify_as 8 ^

-odir C:\ZenmuseL1\corr\4.lasthin2 -olaz ^

-cores 12

lasthin ^

-i C:\ZenmuseL1\corr\4.lasthin2\*.laz ^

-step 0.5 -percentile 10 20 -classify_as 8 ^

-odir C:\ZenmuseL1\corr\5.lasthin3 -olaz ^

-cores 12

:: Operating only on the points classified as 8 (i.e. ignoring those classified as 1) we then run a ground classification with lasground using the following command line, which creates a “low ground” classification ::

lasground ^

-i C:\ZenmuseL1\corr\5.lasthin3\*.laz ^

-ignore_class 1 ^

-nature -ultra_fine ^

-ground_class 2 ^

-odir C:\ZenmuseL1\corr\6.lasground -olaz ^

-cores 12

:: Using lasheight we then create a “thick ground” by pulling all those points into the ground surface that are between 5 centimeter below and 15 centimeter above the “low ground”. For visualization purposes we temporarily use class 6 to capture this thickened ground ::

lasheight ^

-i C:\ZenmuseL1\corr\6.lasground\*.laz ^

-classify_between -0.05 0.15 6 ^

-odir C:\ZenmuseL1\corr\7.lasheight -olaz ^

-cores 12

:: The “thick ground” is showing orange ::

Image 6. Thick ground

:: We go back to lasthin and reclassify in every 50 cm by 50 cm area the point closest to the 50th percentile as class 8. This is what we call the “median ground" ::

lasthin ^

-i C:\ZenmuseL1\corr\7.lasheight\*.laz ^

-ignore_class 1 ^

-step 0.5 -percentile 50 -classify_as 8 ^

-odir C:\ZenmuseL1\corr\8.lasthin4 -olaz ^

-cores 12

:: These are the points we will use to eventually compute the DTM ::

Image 7. Median ground (Brown color)

:: We complete the fully automated classification available in LAStools by running lasclassify with the following options. See the README file for what these options mean. Note that we move the “thick ground” from the temporary class 6 to the proper class 2. The “median ground” continues to be in class 8 ::

lasclassify ^

-i C:\ZenmuseL1\corr\8. lasthin4\*.laz ^

-change_classification_from_to 6 2 ^

-rugged 0.3 -ground_offset 1.5 ^

-odir C:\ZenmuseL1\corr\9. lasclassify -olaz ^

-cores 12

Image 8. Point cloud classified

:: We should remove the temporary buffers, which is done with lastile – the same tool that created the buffers ::

lastile ^

-i C:\ZenmuseL1\corr\9.lasclassify\*.laz ^

-remove_buffer ^

-odir C:\ZenmuseL1\corr\10.remove_buffer -olaz ^

-cores 12

:: Now we produce the standard product Digital Terrain Model (DTM) and Digital Surface Model (DSM) at a resolution of 0.25 cm. Because the total area is not that big we generate temporary tiles in “raster LAZ” with las2dem and merge them into a single *.asc file with blast2dem ::

:: DSM ::

lasthin ^

-i C:\ZenmuseL1\corr\10.remove_buffer\*.laz ^

-step 0.25 -percentile 95 ^

-odir C:\ZenmuseL1\corr\12.lasthin5 -olaz ^

-cores 12

las2dem ^

-i C:\ZenmuseL1\corr\12.lasthin5\*.laz ^

-step 0.25 -use_tile_bb ^

-odir C:\ZenmuseL1\corr\13.dsm\tiles_dsm_25cm -olaz ^

-cores 12

blast2dem ^

-i C:\ZenmuseL1\corr\13.dsm\*.laz -merged ^

-step 0.25 -hillshade ^

-o C:\ZenmuseL1\corr\13.dsm\dsm.asc

Image 9. Digital Surface Model (DSM)

:: DTM ::

las2dem ^

-i C:\ZenmuseL1\corr\10.remove_buffer\*.laz ^

-keep_class 8 ^

-step 0.25 -use_tile_bb ^

-odir C:\ZenmuseL1\corr\11.dtm\tiles_dtm_25cm -olaz ^

-cores 12

blast2dem ^

-i C:\ZenmuseL1\corr\11.dtm\*.laz -merged ^

-step 0.25 -hillshade ^

-o C:\ZenmuseL1\corr\11.dtm\dtm_25cm.asc

Image 10. Digital Terrain Model (DTM)

I hope this LiDAR workflow using LAStools and Bayesmap is helpful for you, Thanks for read!.

Special thanks! to:

Andre Jalobeanu, PhD from Bayesmap, LLC (CEO of Bayesmap)

Chester P. Walker, PhD from Archaeo-Geophysical Associates, LLC (Owner of the data)

Contacts:

LiDAR Latinoamérica serves as the official distributor for Bayesmap in all Latin American countries, and USA. For commercial inquiries, please contact us at [email protected] and [email protected]