How to check if grid cells are empty?

555 views
Skip to first unread message

Ellon Mendes

unread,
Oct 24, 2017, 12:19:53 PM10/24/17
to LAStools - efficient tools for LiDAR processing
    Hello Martin and all,

I'm wondering if there is a way to count the number of cells in a grid over a LAS dataset that contains zero pulses. To give some context, I want to use small step sizes in lasgrid and lascanopy, and I want to find a compromise between the step size and the possible number of empty cells inside the dataset.

From what I tested, cells with no points are not present in the output for the XYZ format, and other formats as ASC or BIL add them as nodata, but since my bounding box is not completely filled with points (see attached image) I have many regions with no data, so I can not assess the number of empty cells inside the data.

For example, consider we use the following with the lake.laz file inside LASTools/data (which have a region in the middle without points):

lasgrid -i lake.laz -last_only -step 1 -counter -oxyz

python
Python 2.7.6 (default, Oct 26 2016, 20:30:19)
[GCC 4.8.4] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> lake_count = np.loadtxt('lake.xyz',delimiter=',',usecols=(2,))
>>> lake_count.all()
True

It means that there are no cells with zero points in the XYZ, even if we know that there is a hole in the dataset (the lake).

I also tried to get the output in CSV format to see if I can get cells with zero points, but if I use -ocsv I get an error:

lasgrid -i lake.laz -last_only -step 1 -counter -ocsv
ERROR: zero column pointer not supported by SRwriter_csv
ERROR: cannot open SRwriter_csv
ERROR: could not open srwriter

Is there any way to check if there are empty cells with LASTools?


Best,

Ellon
data_bounding_box.png

Tobias K Kohoutek

unread,
Oct 25, 2017, 7:45:52 AM10/25/17
to LAStools - efficient tools for LiDAR processing
Hi Ellon,

you could use:

lasgrid -i "C:\LAStools\data\lake.laz" -point_density -false -set_min_max 1 4 -odir "C:\LAStools\data\test" -otif

than load your *.tif into into your favorite GIS software, which should give you the option to "Create Area Features from Equal Values ..." and then you select the color "BLACK" and a new layer is created including only those pixel (1m raster) without information.



Cheers,
Tobias

Perdue, Michael

unread,
Oct 25, 2017, 10:05:29 PM10/25/17
to last...@googlegroups.com

Hello,

 

You can use a grid format that supports “data voids” by the use of a nodata value (bil, tif, etc). After that, you can use a separate call to gdal_translate to remove the no data value assigned in the header by lasgrid. Alternately you can just count the number of nodata cells (in python/GIS/whatever) or manually edit the nodata value in the header if you use bil format.

 

> lasgrid -counter -last_only -i lake.laz -o lake_grid.bil

 

> gdal_translate -a_nodata none -of EHdr lake_grid.bil lake_grid_nodata.bil

 

> lasinfo -i lake_grid_nodata.bil -histo z 1

…

z coordinate histogram with bin size 1

  bin 0 has 27878

  bin 1 has 13250

  bin 2 has 14512

  bin 3 has 8536

  bin 4 has 3145

  bin 5 has 910

  bin 6 has 323

  bin 7 has 170

  bin 8 has 89

  bin 9 has 68

  bin 10 has 51

bin 11 has 41

  bin 12 has 30

  bin 13 has 32

  bin 14 has 15

  bin 15 has 11

  bin 16 has 8

  bin 17 has 17

  bin 18 has 5

  bin 19 has 6

bin 20 has 7

  bin 21 has 7

  bin 22 has 1

  bin 23 has 1

  bin 24 has 2

  bin 25 has 3

  bin 26 has 2

  bin 27 has 3

  bin 28 has 3

  bin 29 has 2

  bin 30 has 1

  bin 31 has 3

  bin 33 has 1

  bin 34 has 2

  bin 35 has 2

  bin 38 has 1

  bin 39 has 2

bin 41 has 2

  bin 44 has 1

  bin 50 has 1

  average z coordinate 1.35244 for 69144 element(s)

 

Cheers,

 

Mike

Carlos Alberto Silva

unread,
Oct 26, 2017, 2:04:02 AM10/26/17
to last...@googlegroups.com

Hi everyone,

Is there any workflow of how thinning lidar data by reducing pulse density (not point density) systematically based on GPS time in Lastools?

Thanks a lot!


Best Regards,

Carlos A. Silva
Forest Engineer 
MSc Forest Resource 
Ph.D. Student at University of Idaho
Current Research: LiDAR Remote Sensing of Forest Resources
USDA Forest Service
Rocky Mountain Research Station (RMRS)
Moscow, Idaho, USA
Office Phone Number: (208) 882-3558

Web-LiDAR forest inventory applications:

rLiDAR: An R package for reading, processing and visualizing LiDAR (Light Detection and Ranging) data.

rForest: An R package for Forest Inventory Analysis


 

Tobias K Kohoutek

unread,
Oct 26, 2017, 7:35:20 AM10/26/17
to LAStools - efficient tools for LiDAR processing
Dear Carlos,

Yes there is!

Sorry if this sounds a bit harsh, but:
First of all, you're out of topic here in this post.
Second, when you write in your signature, that you're a Ph.D. Student, people normally expact, that you at least use the search function in a forum/group, which would have helped you to find at least those 2 posts:


We're all here to help each other, but research, and this is what a Ph.D. is about, includes "search" which isn't equal to "I'll ask and see who will tell me what to do."

Once again I'm sorry is this sounds harsh in a group forum, but check the links and the search option and you'll find algorithms. When you've specific question towards those algorithms, that aren't clear or give you errors, of course in this way you should ask.

Cheers,
Tobias

Ellon Mendes

unread,
Oct 30, 2017, 11:39:03 AM10/30/17
to last...@googlegroups.com
    Hello Tobias and Michael,

Thanks for your answers. I managed to do a solution extending Michael's idea with a gdalwap call to crop the bil file to the extent of my target area. Here's it is an example of the work-flow, which I repeat for several step values in a loop:
  1. grid into a bil file:

    lasgrid -i all_14_Riegl.las
            -last_only
            -counter_16bit
            -step 1.95
            -odir find_best_step_tmp
            -obil
            -odix _1.95
            -nodata 0
            -nbits 16

  2. change no data to none:


  1. gdal_translate -a_nodata none
                   -of EHdr
  1.                find_best_step_tmp/all_14_Riegl_1.95_nodata.bil
                   find_best_step_tmp/all_14_Riegl_1.95_cropped.bil

  2. crop to field extents using gdalwarp and a shapefile. Cells outside the shapefile will be set to 65535, which will be the value used for nodata:

    gdalwarp --config SHAPE_RESTORE_SHX true
             -cutline corn_field.shp
             -crop_to_cutline
             -dstnodata 65535
             -of EHdr
             find_best_step_tmp/all_14_Riegl_1.95_nodata.bil
             find_best_step_tmp/all_14_Riegl_1.95_cropped.bil

  3. create histograms:

    lasinfo -i find_best_step_tmp/all_14_Riegl_1.95_cropped.bil
            -histo z 1
            -odix _count_histo
            -otxt

Thanks for the help!

Best,

Ellon

Martin Isenburg

unread,
Oct 31, 2017, 12:07:48 PM10/31/17
to LAStools - efficient command line tools for LIDAR processing
Hello Ellon,

Here another approach. For tiles that are not fully filled this will require a few steps. In the first step you compute - as usual - a raster of densities. 

:: create a suitable test file by keeping just one flight line

E:\LAStools\bin>las2las -i ..\data\france.laz 
                                       -keep_point_source 3 ^
                                      -o france_fl3.laz

:: have a look at the result (see attached pic 1)

E:\LAStools\bin>lasview -i france_fl3.laz

:: compute a 1m density grid of last returns in BIL format

E:\LAStools\bin>lasgrid -i france_fl3.laz ^
                                      -keep_last ^
                                      -step 1 -point_density ^
                                      -odix _density -obil

:: have a look at the result (see attached pic 2)

E:\LAStools\bin>lasview france_fl3_density.bil

:: trick: create a raster of zeros for the area covered by the LiDAR

E:\LAStools\bin>las2dem -i france_fl3.laz ^
                                         -keep_last ^
                                         -set_intensity 0 ^
                                         -step 1 -intensity ^
                                         -odix _zero -obil

:: have a look at the result (see attached pic 3)

E:\LAStools\bin>lasview -i france_fl3_zero.bil

:: merge the two rasters keeping the higher value per grid cell

E:\LAStools\bin>lasgrid -i france_fl3_density.bil ^
                                      -i france_fl3_zero.bil ^
                                      -merged ^
                                      -step 1 -highest ^
                                      -o france_fl3_merged.bil

:: have a look at the result (see attached pic 4)

E:\LAStools\bin>lasview -i france_fl3_merged.bil

:: print a histogram of point densities

E:\LAStools\bin>lasinfo -i france_fl3_merged.bil ^
                                      -nh -nv -nmm -nr ^
                                      -histo z 1
[...]
z coordinate histogram with bin size 1
  bin 0 has 1006
  bin 1 has 1200
  bin 2 has 2232
  bin 3 has 924
  bin 4 has 751
  bin 5 has 200
  bin 6 has 121
  bin 7 has 49
  bin 8 has 36
  bin 9 has 23
  bin 10 has 10
  bin 11 has 9
  bin 12 has 11
  bin 13 has 2
  bin 14 has 1
  bin 16 has 3
  bin 17 has 1
  bin 20 has 1
  bin 21 has 1
  bin 24 has 1
[...]

:: but this overestimates the zero grid cells along the oversimplified border
:: trick: create a *tighter* raster of zeros for the area covered by the LiDAR

E:\LAStools\bin>las2dem -i france_fl3.laz ^
                                         -keep_last ^
                                         -kill 10 ^
                                         -set_intensity 0 ^
                                         -step 1 -intensity ^
                                         -odix _zero_tight -obil

:: have a look at the result (see attached pic 5)

E:\LAStools\bin>lasview -i france_fl3_zero_tight.bil

:: merge the two rasters keeping the higher value per grid cell

E:\LAStools\bin>lasgrid -i france_fl3_density.bil ^
                                      -i france_fl3_zero_tight.bil ^
                                      -merged ^
                                      -step 1 -highest ^
                                      -o france_fl3_merged_tight.bil

:: have a look at the result (see attached pic 6)

E:\LAStools\bin>lasinfo -i france_fl3_merged_tight.bil ^
                                      -nh -nv -nmm -nr ^
                                      -histo z 1
[...]
z coordinate histogram with bin size 1
  bin 0 has 438
  bin 1 has 1200
  bin 2 has 2232
  bin 3 has 924
  bin 4 has 751
  bin 5 has 200
  bin 6 has 121
  bin 7 has 49
  bin 8 has 36
  bin 9 has 23
  bin 10 has 10
  bin 11 has 9
  bin 12 has 11
  bin 13 has 2
  bin 14 has 1
  bin 16 has 3
  bin 17 has 1
  bin 20 has 1
  bin 21 has 1
  bin 24 has 1
[...]

Regards,

Martin @rapidlasso

lasgrid_las2dem_lasinfo_finding_empty_grid_cells_pic1.png
lasgrid_las2dem_lasinfo_finding_empty_grid_cells_pic2.png
lasgrid_las2dem_lasinfo_finding_empty_grid_cells_pic3.png
lasgrid_las2dem_lasinfo_finding_empty_grid_cells_pic4.png
lasgrid_las2dem_lasinfo_finding_empty_grid_cells_pic5.png
lasgrid_las2dem_lasinfo_finding_empty_grid_cells_pic6.png

Jeffrey Osborne

unread,
Feb 9, 2018, 4:32:06 PM2/9/18
to LAStools - efficient tools for LiDAR processing
Has anyone had trouble running las2dem on a large number of files? If I follow Martin's example workflow on a directory of many thousands of tiles las2dem always stops working after it processes 2043 files. It is not due to a bad tile as I've tested this on different datasets and this always happens no matter what file the 2044th happens to be. My workaround so far is to separate the tiles in directories of <2044 files and run the command multiple times but I'd like to find out why las2dem is behaving this way. I've also tried using a list of files or *.las as the input and it fails for each.   

las2dem.exe -lof filelist.txt -keep_last -kill 10 -set_intensity 0 -step 1 -intensity -odix _zero_tight -obil -odir D:\LAS2DEM_Test\bil  
las2dem_error.png

Martin Isenburg

unread,
Feb 11, 2018, 6:04:24 PM2/11/18
to LAStools - efficient command line tools for LIDAR processing
Hello Jeffrey,

maybe there is a bug in my BIL writer that forgets to close all the files that a BIL is composed of and I believe the maximal number of files that can be open for a process on a Windows Operating System is around that number where it stops for you. If instead you generate TIF or IMG or DTM or ASC or FLT or XYZ rasters ... does the issue with the stopping after 2043 files go away?

Also ... simply adding a second code (or a few more) with '-cores 3' (assuming your PC/laptop has multiple CPUs) will not only circumvent the bug from occurring but also speed up your DEM creation.

Regards,

Martin

Jeffrey Osborne

unread,
Feb 13, 2018, 6:31:44 AM2/13/18
to LAStools - efficient tools for LiDAR processing
Hi Martin, thanks for the suggestions. I tested the different outputs and they all stop after 2043 files. I'll look at using the cores option however it would be nice if a single process could run on all the files for cases where the number of input files would exceed the number of cores I have available.

Jeffrey Osborne

unread,
Mar 9, 2018, 1:56:13 PM3/9/18
to LAStools - efficient tools for LiDAR processing
I've also noticed a similar problem with lasgrid. If running a single core on a large number of files it stops after processing 507 files.
Reply all
Reply to author
Forward
0 new messages