point density with lasinfo

[Alessandro Montaghi]

The definition of point density given by Prof. Norbert Pfeifer (TU Wien) is:

Point density: Average number of recorded points per square meter within an analysis unit (in the case of multiple echoes per pulse only the points defined by the last echo have to be considered). If the strip overlap is smaller than 50% or in the case of crosswise strips that do not cover the whole analysis unit, the point density has to be estimated with the help of the points of the singular cover. Areas with data holes have to be excluded from the computation and have to be considered in a different way.

I wish to know if somebody know how lasinfo computes and reports the estimate of the point density

lasinfo -i lidar.las -compute_density

1) the algorithm considers all points or only last echoes (following the Norbert Pfeifer's definition) ?

2) lasinfo uses the average number of recorded points per square meter  ( 1 m^2) or  a different grid size (default setting)?

3) is It possibile to set the size in the command lasinfo -i lidar.las -compute_density?

Thanks in advance

Ale

--

Alessandro Montaghi

Just some stuff about climbing:

"Think less, climb more" [Anonymous climber]

"The best route setter is the Nature" [me]

"Climbing is a brilliant form of therapy" [Tim Emmett]

"Promettimi che scalerai" [Silvia Buraggi]

"you are not your job" [a faded poster in J-Tree, California]

"Non vale nemmeno il Gri Gri del Lez" [Alessandro Todeschi]

"La montagna è come un amore: se sei respinto, è meglio tornare indietro e non insistere" [Christian Kuntner]

"Enjoy every single moment" [my Mom]

"Chi più alto sale, più lontano vede; chi più lontano vede, più a lungo sogna" [Walter Bonatti] "Don't eat the yellow snow" [Swedish proverb] "Do what makes you happy" [a yellow post-it in the Copenhagen rail-way station]  "MY HANDS ARE MY GEARS" [Anonymous climber] "When you have time you haven't moneys to climb, when you have moneys to climb you haven't time" [me] “I am losing precious days. I am degenerating into a machine for making money. I am learning nothing in this trivial world of men. I must break away and get out into the mountains to learn the news” [John Muir]

[Roberto Antolín]

Hello Alessandro and all

2013/8/26 Alessandro Montaghi <alessandr...@gmail.com>

Point density: Average number of recorded points per square meter within an analysis unit (in the case of multiple echoes per pulse only the points defined by the last echo have to be considered). If the strip overlap is smaller than 50% or in the case of crosswise strips that do not cover the whole analysis unit, the point density has to be estimated with the help of the points of the singular cover. Areas with data holes have to be excluded from the computation and have to be considered in a different way.

Is this definition published anywhere? Could you please give the reference?

Thank you!

--
Roberto

[Alessandro Montaghi]

Hey Roberto,

This was a Point density definition given by Prof. Norbert Pfeifer during a course that my group and me organized in Sweden last yea for PhD students: NOVA 3D Remote Sensing of Forests, Umeå August 2012

This was the schedule

http://www2.nova-university.org/documents/application/57-111013-436/schedule/57-111013-436_sched.pdf

http://www2.nova-university.org/documents/application/57-111013-436/plan/57-111013-436_plan.pdf

if you are interest about a published definition of "density" you can ask to Prof. Norbert Pfeifer. He is one of the world leader researcher about LiDAR and he is really nice person.

Ale 

--
Download LAStools at
http://lastools.org
http://rapidlasso.com
Be social with LAStools at
http://facebook.com/LAStools
http://twitter.com/LAStools
http://linkedin.com/groups/LAStools-4408378
Manage your settings at
http://groups.google.com/group/lastools/subscribe

--

Alessandro Montaghi

Just some sentences about climbing/mountaineering that i heard during my travels:

"Think less, climb more" [Anonymous climber]

"The best route setter is the Nature" [me]

"Climbing is a brilliant form of therapy" [Tim Emmett]

"Promettimi che scalerai" [Silvia Buraggi]

"you are not your job" [a faded poster in J-Tree, California]

"Non vale nemmeno il Gri Gri del Lez" [Alessandro Todeschi]

"La montagna è come un amore: se sei respinto, è meglio tornare indietro e non insistere" [Christian Kuntner]

"Enjoy every single moment" [my Mom]

"Chi più alto sale, più lontano vede; chi più lontano vede, più a lungo sogna" [Walter Bonatti] "Don't eat the yellow snow" [Swedish proverb] "Do what makes you happy" [a yellow post-it in the Copenhagen rail-way station]  "MY HANDS ARE MY GEARS" [Anonymous climber] "When you have time you haven't moneys to climb, when you have moneys to climb you haven't time" [me]

"The Nut is like the wife while the friend is like the lover" [Emiliano Max]

[Martin Isenburg]

Hello Ale,

lasinfo.exe computes the point density for both, all returns (which I think should really be called return density) and last returns only (which i think should really be called pulse density). how is it done? When you request '-compute_density' in the command line then lasinfo.exe lays an occupany grid over the points that has a cell size of 2 by 2 meters if the projection is specified as meters of 2 by 2 units if nothing is specified or of 6 by 6 (survey) feet if the projection is specified as (survey) feet (either implicitly in the VLR or explicitly in the command line using '-feet' (or '-survey_feet')).

At the end lasinfo.exe counts how many of these 2 by 2 = 4 square meter cells (or 6 by 6 = 36 square (survey) feet) are covered with points and divides the total number of returns as well as the total number of last returns by this area. Because of the partly covered 4 square meter cells or 36 square (survey) feet cells along the boundary of the LAS file the area is usually estimated a bit larger than it really is and therefore the point density will be underestimated. However, that assumes the LiDAR system distributes the points evenly on the ground which is not the case ... especially not for zigzag scanners. For those point distributions the average density computed may be overestimated so that lasinfo may then report about the right result. For more info on the zigzag issue see this discussion here:

https://groups.google.com/forum/#!lastools/QmWrlT1fPEo/zd6g9AJrjsMJ

Regards,

Martin @rapidlasso

http://rapidlasso.com - fast tools to check your LiDARs

[GeorgeM]

Hello,

Is that possible to aggregate multiple las file density values somehow to one single file? Example:
lasfile_1.las ... lasfile_n.las resides in a directory, compute all .las files densities and write it to a csv
(or other file) as filename,density:

lasfile_1, 1.246
lasfile_2, 1.189
...
lasfile_n, 1.157

Thanks.

[Perdue, Michael]

This is relatively easy to do in a script. It is a somewhat painful process in DOS, but in bash it would look like;

# Loop through the files with lasinfo -compute_density
# 2>&1 redirects output from stderr to stdout so it can be piped to grep
# Use grep to isolate the line with the needed info and pipe that line to awk
# Use awk to isolate the 8th white space delimited field (last only density)
# the external loop variable ${i} is assigned to an internal awk variable (var)
# so it can be appended to "list_of_Densities.txt"

for i in *.las; do
lasinfo ${i} -compute_density 2>&1 | \
grep "point density" | \
awk -v var="${i}" '{print var, $8}' >> list_of_Densities.txt
done

To get bash and all of it’s utilities onto windows you could install msys or Cygwin.

Alternatively, I’m sure the same thing could easily be done in perl or python, which you may already have installed on your system. There are a number of users on this forum who are fluent in each and can probably provide an example of how to accomplish it in the respective environment.

Cheers,

Mike

[Martin Isenburg]

Hello,

I like to remind folks that the "average point density" that is computed by lasinfo is a very miss-leading number for the purposes of quality control. Especially when there are multiple flight lines or when oscillating or zigzag scanning systems are used the simple computation 

number of last returns / total area

that is done by lasinfo is a very poor indication of the actual distribution of pulses throughout the surveyed area. I suggest to instead compute either density grids with lasgrid

http://rapidlasso.com/lasgrid

http://lastools.org/download/lasgrid_README.txt

with a small enough pixel size this will enable you to tell whether the average density was achieved - for example - by having approximately half the areas of the survey with a density that is above the tender specification (where flightlines are overlapping) but therefore also half the areas of the survey with a density below the tender specification (where flightlines are *not* overlapping). This can be a serious problem in later processing steps, because you will, for example, get great building classification with lasclassify in the area of flightline overlap but very poor classification in the less dense areas.

If even distribution of pulses on the ground is important to you - especially for single tree detection in forests or clean deliniation of building roofs lines - then you really want to measure the pulse to pulse spacing. Sometimes there spacing is dense along the scan line (perpendicular to the flight direction) but too wide across scan lines (along the flight direction). This gets especially bad for oscillating mirror systems when the "high scan angles" where the mirror slows down to turn around are not "cut off" from the delivered swath of LAS returns as discussed in details here:

http://rapidlasso.com/2014/03/20/density-and-spacing-of-lidar/ 

There is a new way to accurately output scan angle histograms with las2tin that is not yet discussed in the blog post mentioned above.

http://lastools.org/download/las2tin_README.txt

Here an example:

D:\LAStools\bin>las2tin -i ..\data\fusa.laz -keep_last -histo_only min_edge_length 0.1

min edge length histogram histogram with bin size 0.1

  bin [0,0.1) has 1335

  bin [0.1,0.2) has 6693

  bin [0.2,0.3) has 14525

  bin [0.3,0.4) has 238134

  bin [0.4,0.5) has 1715

  bin [0.5,0.6) has 593

  bin [0.6,0.7) has 287

  bin [0.7,0.8) has 65

  bin [0.8,0.9) has 15

  bin [0.9,1) has 1

  bin [1,1.1) has 5

  average min edge length histogram 0.330636 for 263368 element(s)

D:\LAStools\bin>las2tin -i ..\data\fusa.laz -keep_last -histo_only max_edge_length 0.1

max edge length histogram histogram with bin size 0.1

  bin [0.1,0.2) has 2

  bin [0.2,0.3) has 27

  bin [0.3,0.4) has 624

  bin [0.4,0.5) has 2277

  bin [0.5,0.6) has 16630

  bin [0.6,0.7) has 36825

  bin [0.7,0.8) has 68327

  bin [0.8,0.9) has 84470

  bin [0.9,1) has 36070

  bin [1,1.1) has 6475

  bin [1.1,1.2) has 2740

  bin [1.2,1.3) has 2148

  bin [1.3,1.4) has 2385

  bin [1.4,1.5) has 902

  bin [1.5,1.6) has 495

  bin [1.6,1.7) has 378

  bin [1.7,1.8) has 1187

  bin [1.8,1.9) has 283

  bin [1.9,2) has 98

  bin [2,2.1) has 90

  bin [2.1,2.2) has 54

  bin [2.2,2.3) has 22

  bin [2.3,2.4) has 34

  bin [2.4,2.5) has 28

  bin [2.5,2.6) has 26

  bin [2.6,2.7) has 22

  bin [2.7,2.8) has 26

  bin [2.8,2.9) has 25

  bin [2.9,3) has 20

  bin [3,3.1) has 14

  bin [3.1,3.2) has 14

  bin [3.2,3.3) has 15

  bin [3.3,3.4) has 13

  bin [3.4,3.5) has 15

  bin [3.5,3.6) has 8

  bin [3.6,3.7) has 12

  bin [3.7,3.8) has 11

  bin [3.8,3.9) has 17

  bin [3.9,4) has 9

  bin [4,4.1) has 13

  bin [4.1,4.2) has 11

  bin [4.2,4.3) has 11

  bin [4.3,4.4) has 11

  bin [4.4,4.5) has 6

  bin [4.5,4.6) has 8

  bin [4.6,4.7) has 9

  bin [4.7,4.8) has 8

  bin [4.8,4.9) has 6

  bin [4.9,5) has 6

  bin [5,5.1) has 9

  bin [5.1,5.2) has 6

  bin [5.2,5.3) has 4

  bin [5.3,5.4) has 8

  bin [5.4,5.5) has 4

  bin [5.5,5.6) has 7

  bin [5.6,5.7) has 4

  bin [5.7,5.8) has 3

  bin [5.8,5.9) has 5

  bin [5.9,6) has 7

  bin [6,6.1) has 4

  bin [6.1,6.2) has 6

  bin [6.2,6.3) has 4

  bin [6.3,6.4) has 8

  bin [6.4,6.5) has 6

  bin [6.5,6.6) has 5

  bin [6.6,6.7) has 4

  bin [6.7,6.8) has 11

  bin [6.8,6.9) has 4

  bin [6.9,7) has 3

  bin [7,7.1) has 7

  bin [7.1,7.2) has 3

  bin [7.2,7.3) has 6

  bin [7.3,7.4) has 5

  bin [7.4,7.5) has 6

  bin [7.5,7.6) has 5

  bin [7.6,7.7) has 3

  bin [7.7,7.8) has 4

  bin [7.8,7.9) has 5

  bin [7.9,8) has 6

  bin [8,8.1) has 6

  bin [8.1,8.2) has 4

  bin [8.2,8.3) has 6

  bin [8.3,8.4) has 2

  bin [8.5,8.6) has 1

  bin [8.6,8.7) has 3

  bin [8.7,8.8) has 6

  bin [8.8,8.9) has 2

  bin [8.9,9) has 3

  bin [9,9.1) has 1

  bin [9.1,9.2) has 6

  bin [9.2,9.3) has 3

  bin [9.3,9.4) has 5

  bin [9.4,9.5) has 2

  bin [9.5,9.6) has 3

  bin [9.6,9.7) has 4

  bin [9.7,9.8) has 3

  bin [9.8,9.9) has 2

  bin [9.9,10) has 4

  bin [10,10.1) has 3

  bin [10.1,10.2) has 4

  bin [10.2,10.3) has 1

  bin [10.3,10.4) has 3

  bin [10.4,10.5) has 4

  bin [10.5,10.6) has 4

  bin [10.6,10.7) has 1

  bin [10.7,10.8) has 2

  bin [10.8,10.9) has 1

  bin [10.9,11) has 4

  bin [11,11.1) has 6

  bin [11.1,11.2) has 4

  bin [11.3,11.4) has 2

  bin [11.4,11.5) has 1

  bin [11.5,11.6) has 2

  bin [11.6,11.7) has 3

  bin [11.7,11.8) has 3

  bin [11.8,11.9) has 3

  bin [11.9,12) has 1

  bin [12,12.1) has 4

  bin [12.2,12.3) has 2

  bin [12.4,12.5) has 3

  bin [12.5,12.6) has 3

  bin [12.6,12.7) has 1

  bin [12.7,12.8) has 6

  bin [12.8,12.9) has 2

  bin [12.9,13) has 1

  bin [13,13.1) has 3

  bin [13.1,13.2) has 1

  bin [13.2,13.3) has 4

  bin [13.3,13.4) has 1

  bin [13.4,13.5) has 4

  bin [13.5,13.6) has 2

  bin [13.6,13.7) has 1

  bin [13.7,13.8) has 3

  bin [13.8,13.9) has 4

  bin [13.9,14) has 1

  bin [14,14.1) has 3

  bin [14.1,14.2) has 1

  bin [14.2,14.3) has 1

  bin [14.3,14.4) has 1

  bin [14.5,14.6) has 2

  bin [14.6,14.7) has 2

  bin [14.7,14.8) has 1

  bin [14.8,14.9) has 2

  bin [14.9,15) has 1

  bin [15,15.1) has 2

  bin [15.4,15.5) has 3

  bin [15.7,15.8) has 1

  bin [15.8,15.9) has 1

  bin [16,16.1) has 1

  bin [16.1,16.2) has 2

  bin [16.3,16.4) has 1

  bin [16.5,16.6) has 2

  bin [16.6,16.7) has 1

  bin [16.9,17) has 1

  bin [17.2,17.3) has 1

  bin [17.3,17.4) has 1

  bin [17.4,17.5) has 1

  bin [17.7,17.8) has 2

  bin [17.8,17.9) has 1

  bin [17.9,18) has 1

  bin [18,18.1) has 1

  bin [18.1,18.2) has 1

  bin [18.2,18.3) has 1

  bin [18.3,18.4) has 2

  bin [18.4,18.5) has 2

  bin [18.5,18.6) has 2

  bin [18.7,18.8) has 1

  bin [18.8,18.9) has 3

  bin [18.9,19) has 1

  bin [19.1,19.2) has 1

  bin [19.4,19.5) has 1

  bin [19.5,19.6) has 2

  bin [19.6,19.7) has 1

  bin [19.8,19.9) has 2

  bin [19.9,20) has 1

  bin [20.2,20.3) has 1

  bin [20.6,20.7) has 1

  bin [21,21.1) has 1

  bin [21.1,21.2) has 3

  bin [21.3,21.4) has 1

  bin [21.7,21.8) has 1

  bin [21.8,21.9) has 2

  bin [21.9,22) has 1

  bin [22.2,22.3) has 3

  bin [22.4,22.5) has 2

  bin [22.5,22.6) has 2

  bin [22.9,23) has 2

  bin [23,23.1) has 1

  bin [23.4,23.5) has 1

  bin [23.9,24) has 1

  bin [24.4,24.5) has 2

  bin [24.5,24.6) has 1

  bin [24.9,25) has 1

  bin [25.2,25.3) has 1

  bin [25.5,25.6) has 1

  bin [25.9,26) has 1

  bin [26.8,26.9) has 1

  bin [27.2,27.3) has 2

  bin [27.6,27.7) has 1

  bin [28,28.1) has 1

  bin [28.3,28.4) has 1

  bin [28.6,28.7) has 1

  bin [29.6,29.7) has 2

  bin [30.7,30.8) has 1

  bin [32.9,33) has 1

  bin [33.4,33.5) has 1

  bin [33.7,33.8) has 1

  bin [34.1,34.2) has 1

  bin [34.4,34.5) has 1

  bin [34.8,34.9) has 1

  bin [35.1,35.2) has 1

  bin [35.5,35.6) has 1

  bin [35.7,35.8) has 1

  bin [35.9,36) has 1

  bin [36.2,36.3) has 1

  bin [36.4,36.5) has 2

  bin [36.5,36.6) has 1

  bin [36.9,37) has 1

  bin [37,37.1) has 2

  bin [37.2,37.3) has 1

  bin [37.6,37.7) has 1

  bin [37.9,38) has 1

  bin [38.3,38.4) has 1

  bin [38.6,38.7) has 2

  bin [38.8,38.9) has 1

  bin [39.2,39.3) has 2

  bin [39.9,40) has 2

  bin [41.6,41.7) has 2

  bin [43.4,43.5) has 2

  bin [44.9,45) has 1

  bin [46,46.1) has 2

  bin [49.4,49.5) has 2

  average max edge length histogram 0.838941 for 263368 element(s)

Regards,

Martin @rapidlasso

[Jeffrey Osborne]

Does the occupancy grid that gets created observe any filters that may be applied? For example if I drop a water class when calculating the density is the area of a water-body (assuming all returns within it are classified as water) included or excluded in the total area, or is only the number of water returns subtracted from the number of returns when the density calculation is performed such that the area of the water-body will cause the density to be underestimated as well.

[Martin Isenburg]

Hello,

All of the filters listed when you type, for example, "lasview -h" are applied *before* any LAStools is actually given the point cloud so that all *filtered out* points are for all practical purposes not there. An exception is the '-filtered_transform' option that gives the filter to the transform engine and does not filter out any points.

The area will be subtracted but only in increments of 2 by 2 meter cells (or 6 by 6 feet). A single non-water return in a cell means the cells entire area is counted. If the pulse density calculation is of great concern to use then create a '-point_density' raster of the last returns with lasgrid and evaluate the result manually.

http://rapidlasso.com/lasgrid

http://lastools.org/download/lasgrid_README.txt

Here an example (with buildings taking the role of water). First we create a last return density raster for each 2.5 by 2.5 meter cell. 

lasgrid -i ..\data\fusa.laz ^

            -drop_class 6 ^

            -keep_last ^

            -point_density -step 2.5 ^

            -odix _d_250cm -oasc

Below we create a histogram of the values stored in each raster cell using lasinfo by loading the ASC as if it was LAS and ask for a histogram of the z coordinate as this is where the raster value gets stored when an ASC (or BIL or DTM) raster is read with on-the-fly conversion to LAS by LASlib. Note that only 8945 cells are occupied of this 100 by 100 = 10000 cell raster. The water (aka buildings) are leaving holes in the raster.

You can also visualize the point density with lasview and you'll see that the occasional non-water points (aka non-building points) along the river banks (aka roof edges) are counted and lower the average. You can apply filters '-drop_z_below 0.5' as well as '-drop_z_above 6' to kick out the outliers from the average computation. Or just report the entire histogram.

lasview -i ..\data\fusa_d_250cm.asc

Also see this blog post:

http://rapidlasso.com/2014/03/20/density-and-spacing-of-lidar/

lasinfo -i ..\data\fusa_d_250cm.asc ^

           -nh -nv -nr -nmm ^

           -histo z 0.1

lasinfo (170322) report for ..\data\fusa_d_250cm.asc

WARNING: for return 1 real number of points by return is 8945 but header entry was not set.

z coordinate histogram with bin size 0.1

  bin [0.1,0.2) has 82

  bin [0.3,0.4) has 55

  bin [0.4,0.5) has 34

  bin [0.6,0.7) has 47

  bin [0.8,0.9) has 35

  bin [0.9,1) has 49

  bin [1.1,1.2) has 62

  bin [1.2,1.3) has 53

  bin [1.4,1.5) has 24

  bin [1.6,1.7) has 23

  bin [1.7,1.8) has 48

  bin [1.9,2) has 50

  bin [2,2.1) has 58

  bin [2.2,2.3) has 75

  bin [2.4,2.5) has 95

  bin [2.5,2.6) has 81

  bin [2.7,2.8) has 96

  bin [2.8,2.9) has 108

  bin [3,3.1) has 130

  bin [3.2,3.3) has 187

  bin [3.3,3.4) has 705

  bin [3.5,3.6) has 673

  bin [3.6,3.7) has 469

  bin [3.8,3.9) has 352

  bin [4,4.1) has 294

  bin [4.1,4.2) has 335

  bin [4.3,4.4) has 382

  bin [4.4,4.5) has 1365

  bin [4.6,4.7) has 1230

  bin [4.8,4.9) has 779

  bin [4.9,5) has 402

  bin [5.1,5.2) has 161

  bin [5.2,5.3) has 80

  bin [5.4,5.5) has 77

  bin [5.6,5.7) has 79

  bin [5.7,5.8) has 67

  bin [5.9,6) has 28

  bin [6,6.1) has 32

  bin [6.2,6.3) has 13

  bin [6.4,6.5) has 3

  bin [6.5,6.6) has 7

  bin [6.7,6.8) has 5

  bin [6.8,6.9) has 3

  bin [7,7.1) has 2

  bin [7.2,7.3) has 2

  bin [7.3,7.4) has 2

  bin [7.5,7.6) has 2

  bin [7.6,7.7) has 1

  bin [8,8.1) has 1

  bin [8.3,8.4) has 1

  bin [8.9,9) has 1

  average z coordinate 3.97626 for 8945 element(s)

--

[Jeffrey Osborne]

Thanks Martin, that is the behavior I was expecting but I wanted to be certain. I've seen that blog post and have created density rasters to avoid relying on the lasinfo density calculation alone. I hadn't thought to do a histogram of the rasters in that way, that is very helpful as I wanted to do a % of data below a particular density and that will get me there. One other question I was wondering about is the difference between the 16bit and 32bit counters for the pulse density. When would they be used rather than just the 8bit counter?

[Martin Isenburg]

Hello,

 One other question I was wondering about is the difference between the 16bit and 32bit counters for the pulse density. When would they be used rather than just the 8bit counter?

If you expect your cells to contain more then 255 last returns (e.g. you use large cells or high point densities) then you should use the 16 bit version as that can count up to 65535 last returns per cell. You'll rarely need the 32 bit ones. For pulse densities of 20 a cell size of 5 meters would result in an average of 20 * 5 * 5 = 500 last returns per cell so an 8 bit counter would overflow so '-pulse_density_16bit' is needed. But for the example shown with a average pulse density of 5 and a cell size of 2.5 meters we mostly count on average around 5 * 2.5 * 2.5 = 31.25 last returns per cell so that 8 bit should be fine. The memory requirements of lasgrid essentially double when you use wider (more bit) counters.

Regards,

Martin