Pit-free CHM using lastools

[Irf]

Hi Martin,

Thanks for responding. I am trying to create the pit-free CHM. Before we create the series of CHMs using different height thresholds, we need to normalise the lidar points and remove points above a certain height. My data is one laz file with about 134m points. I got the 'abnormal program termination' message, when I attempted to thin the data with the following:

lasthin -i lidar_normalised_below50m.laz -subcircle 0.1 -step 0.16667 -highest -o temp.laz

I followed some of your posts on this, from which I got the impression that it had something to do with the memory allocation. Coarser steps (e.g -step 1) worked, but I am interested in finer resolution. I also found that las2las has the -thin_with_grid argument. My question first question is that:

Q1: In case of a big file, can we do the following instead of using lasthin?

#Thinning

las2las -i lidar_normalised_below50m.laz -thin_with_grid 0.16667 -o thinned.laz

lasthin -i thinned.laz -subcircle 0.1 -o thinned_splat.laz

but we can't add the argument -highest here.

I also wanted to know that:

Q2: How is thin_with_grid (in las2las) different from -step x (in lasthin)?

My third question is about the pit-free CHM. The journal article explicitly mentions the use of 'first returns' in each height threshold; however "the pit-free method of Khosravipour et al." at the end of this blog and the "pit_free_script.bat" in "example_batch_scripts" (pasted below in light blue colour) do not use the '-keep_first' argument when using las2dem (please note here that my data is classified: report.txt attached). I am curious to know that:

Q3: Does las2dem / blast2dem create the raster surfaces (e.g. chm_02.bil. chm_05.bil...chm_65.bil)  using the 'first returns' in their respective height threshold, even if -keep_first is not explicitly mentioned in the arguments?

I tested both, using -keep_first and the blog method, and the results are different. -keep_first created a coarser CHM with many pits as compared to not using the argument -keep_first.

To create a CHM, I am using the script below. 

mkdir tmp_chm

lasheight -i lidar_normalised.laz -drop_z_above 50 -odix _below50 -olaz

## no thinning because the point density was already low (report.txt attached)

## instead of -drop_z_below, I used -keep_class, because the data was already classified

blast2dem -i lidar_normalised_below50.laz -keep_class 2 -step 1 -o tmp_chm/chm_ground.bil

blast2dem -i lidar_normalised_below50.laz -step 1 -kill 1.0 -o tmp_chm/chm_00.bil

blast2dem -i lidar_normalised_below50.laz -keep_class 4 5 -step 1 -kill 1.0 -o tmp_chm/chm_02.bil

blast2dem -i lidar_normalised_below50.laz -keep_class 5 -step 1 -kill 1.0 -o tmp_chm/chm_05.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 10 -step 1 -kill 1.0 -o tmp_chm/chm_10.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 15 -step 1 -kill 1.0 -o tmp_chm/chm_15.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 20 -step 1 -kill 1.0 -o tmp_chm/chm_20.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 25 -step 1 -kill 1.0 -o tmp_chm/chm_25.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 30 -step 1 -kill 1.0 -o tmp_chm/chm_30.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 35 -step 1 -kill 1.0 -o tmp_chm/chm_35.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 40 -step 1 -kill 1.0 -o tmp_chm/chm_40.bil

blast2dem -i lidar_normalised_below50.laz -drop_z_below 45 -step 1 -kill 1.0 -o tmp_chm/chm_45.bil

lasgrid -i tmp_chm/chm*.bil -merged -step 1 -highest -o tmp_chm/chm_pit_free_TEST.bil

lasgrid -i tmp_chm/chm*.bil -merged -step 1 -highest -false -set_min_max 0 50 -o tmp_chm/chm_pit_free_TEST.png

## rmdir tmp_chm

Thank you for bearing the length of my email :)

Kind regards,

Irfan

"pit_free_script.bat" in "example_batch_scripts"

rmdir temp_dir /s /q

mkdir temp_dir

las2dem -i normalized.laz -drop_z_above 0.1 -step %STEP% -o temp_dir/chm_ground.bil

lasthin -i normalized.laz -subcircle 0.025 -step %HALF_STEP% -highest -o temp.laz

las2dem -i temp.laz -step %STEP% -kill %KILL% -o temp_dir/chm_00.bil

las2dem -i temp.laz -drop_z_below  2 -step %STEP% -kill %KILL% -o temp_dir/chm_02.bil

las2dem -i temp.laz -drop_z_below  5 -step %STEP% -kill %KILL% -o temp_dir/chm_05.bil

las2dem -i temp.laz -drop_z_below 10 -step %STEP% -kill %KILL% -o temp_dir/chm_10.bil

las2dem -i temp.laz -drop_z_below 15 -step %STEP% -kill %KILL% -o temp_dir/chm_15.bil

las2dem -i temp.laz -drop_z_below 20 -step %STEP% -kill %KILL% -o temp_dir/chm_20.bil

las2dem -i temp.laz -drop_z_below 25 -step %STEP% -kill %KILL% -o temp_dir/chm_25.bil

las2dem -i temp.laz -drop_z_below 30 -step %STEP% -kill %KILL% -o temp_dir/chm_30.bil

las2dem -i temp.laz -drop_z_below 35 -step %STEP% -kill %KILL% -o temp_dir/chm_35.bil

las2dem -i temp.laz -drop_z_below 40 -step %STEP% -kill %KILL% -o temp_dir/chm_40.bil

las2dem -i temp.laz -drop_z_below 45 -step %STEP% -kill %KILL% -o temp_dir/chm_45.bil

las2dem -i temp.laz -drop_z_below 50 -step %STEP% -kill %KILL% -o temp_dir/chm_50.bil

las2dem -i temp.laz -drop_z_below 55 -step %STEP% -kill %KILL% -o temp_dir/chm_55.bil

las2dem -i temp.laz -drop_z_below 60 -step %STEP% -kill %KILL% -o temp_dir/chm_60.bil

las2dem -i temp.laz -drop_z_below 65 -step %STEP% -kill %KILL% -o temp_dir/chm_65.bil

lasgrid -i temp_dir/chm*.bil -merged -step %STEP% -highest -false -set_min_max 0 65 -o chm_pit_free_d05.png

rmdir temp_dir /s /q

[Terje Mathisen]

Irf wrote:
> Hi Martin,
>
> Thanks for responding. I am trying to create the pit-free CHM. Before
> we create the series of CHMs using different height thresholds, we
> need to normalise the lidar points and remove points above a certain
> height. My data is one laz file with about 134m points. I got the
> 'abnormal program termination' message, when I attempted to thin the
> data with the following:
>
> lasthin -i lidar_normalised_below50m.laz -subcircle 0.1 -step 0.16667
> -highest -o temp.laz

I'm pretty confident that if you change to a lastile-based approach you
will get better and faster results. :-)

Terje

>
> I followed some of your posts on this, from which I got the impression
> that it had something to do with the memory allocation. Coarser steps
> (e.g -step 1) worked, but I am interested in finer resolution. I also
> found that las2las has the -thin_with_grid argument. My question first
> question is that:
>
> Q1: In case of a big file, can we do the following instead of using
> lasthin?
>
> #Thinning

> las2las -i lidar_normalised_below50m.laz -thin_with_grid 0.16667 -o

> thinned.laz
> lasthin -i thinned.laz -subcircle 0.1 -o thinned_splat.laz
>
> but we can't add the argument -highest here.
>
> I also wanted to know that:
>
> Q2: How is thin_with_grid (in las2las) different from -step x (in
> lasthin)?
>

> My third question is about the pit-free CHM. The journal article
> <https://www.researchgate.net/publication/273663100_Generating_Pit-free_Canopy_Height_Models_from_Airborne_Lidar>Â explicitly
> mentions the use of 'first returns' in each height threshold; however
> "the pit-free method of Khosravipour et al." at the end of this blog
> <http://rapidlasso.com/2014/11/04/rasterizing-perfect-canopy-height-models-from-lidar/>Â and

> the "pit_free_script.bat" in "example_batch_scripts" (pasted below in
> light blue colour) do not use the '-keep_first' argument when using
> las2dem (please note here that my data is classified: report.txt
> attached). I am curious to know that:
>
> Q3: Does las2dem / blast2dem create the raster surfaces (e.g.

> chm_02.bil. chm_05.bil...chm_65.bil) Â using the 'first returns' in

> their respective height threshold, even if -keep_first is not
> explicitly mentioned in the arguments?
>

> I tested both, using -keep_first and the blog
> <http://rapidlasso.com/2014/11/04/rasterizing-perfect-canopy-height-models-from-lidar/>Â method,

> and the results are different. -keep_first created a coarser CHM with
> many pits as compared to not using the argument -keep_first.
>

> To create a CHM, I am using the script below.Â

> *"pit_free_script.bat" in "example_batch_scripts"
> *

> rmdir temp_dir /s /q
> mkdir temp_dir
> las2dem -i normalized.laz -drop_z_above 0.1 -step %STEP% -o
> temp_dir/chm_ground.bil
> lasthin -i normalized.laz -subcircle 0.025 -step %HALF_STEP% -highest
> -o temp.laz
> las2dem -i temp.laz -step %STEP% -kill %KILL% -o temp_dir/chm_00.bil

> las2dem -i temp.laz -drop_z_below  2 -step %STEP% -kill %KILL% -o
> temp_dir/chm_02.bil
> las2dem -i temp.laz -drop_z_below  5 -step %STEP% -kill %KILL% -o

> --
> 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


--
- <Terje.M...@tmsw.no>
"almost all programming can be viewed as an exercise in caching"

[Irf]

Hi Terje,

Thanks for the tip. I agree that the lastile approach will work on portion of data and we can merge them later. I would be glad if you shared your thoughts on my other queries too.

Q2: How is thin_with_grid (in las2las) different from -step x (in lasthin)?

Q3: Does las2dem / blast2dem create the raster surfaces (e.g. chm_02.bil. chm_05.bil...chm_65.bil)  using the 'first returns' in their respective height threshold, even if -keep_first is not explicitly mentioned in the arguments?

(PS: Please refer to the details in my OP)

Kind regards,

Irf

[Martin Isenburg]

Hello,

 

Q2: How is thin_with_grid (in las2las) different from -step x (in lasthin)?

The *filter* in LASlib that is available for all LAStools called '-thin_with_grid 2.0' keeps the *first* point (not eliminated by an earlier filter such as '-keep_class 2') that occupies each 2.0 by 2.0 grid cell. The '-step 2.0' of lasthin also sets the grid cell size for a similar grid but then keeps not the first but the '-lowest', the '-highest', the most '-central', or a '-random' point depending on the selected option.

 

Q3: Does las2dem / blast2dem create the raster surfaces (e.g. chm_02.bil. chm_05.bil...chm_65.bil)  using the 'first returns' in their respective height threshold, even if -keep_first is not explicitly mentioned in the arguments?

las2dem and blas2dem triangulate *all* points that are used as input into a TIN surface and then sample the elevation of that TIN surface at the center of each of the grid cells. if you want to use only the first or only the last returns or only the ground plus keypoints then you need to use the approriate filter such as '-keep_first', '-keep_last', or '-keep_class 2 8'. For coarser rasters of a step size of, for example, '-step 2.0' is may be useful to combine this with an on-the-fly thinning '-thin_with_grid 0,5' to lower the number of input points to a reasonable praction and speed up the overall computation as fewer points will have to be triangulated.

Regards,

Martin @rapidlasso

[Irf]

Thanks again Martin. Having read your explanation about the las2dem  and its filters, I wonder why the -first_only filter wasn't used in the 'pit-free' algorithm. Also, instead of the -drop_z_below, shouldn't it be -drop_below? Because as explained in the ReadMe files, -drop_z_below removes the 'elevations' below a certain threshold. Whereas, -drop_below removes the 'computed heights' (after the data is normalised) below a certain threshold. Since the data was already normalised (-replace_z), do you agree that instead of -drop_z_below it should be -drop_below?

Many thanks,

Irf

[Martin Isenburg]

Hello,

The algorithm described on the blog post [1] is a slight variation of the original pit-free CHM algorithm of (Khosravipour et al 2014). Instead of using the first returns it uses the highest returns that are found by lasthin with the '-highest' option on a finer grid after turning all returns into small disks with '-subcircle 0.05'. This "don't use the first returns" idea is somewhat more fleshed out in the SilviLaser 2015 presentation of (Khosravipour et al 2015) on generating pit-free DSM with the new spike-free algorithm that is currently under journal review. Watch for an upcoming blog post on the issue.

About your second question: There are no elevations in a normalized LiDAR file. The z coordinate stores the height and not the elevation after normalization. And those z coordinates will be dropped at the corresponding set of height cutoffs. All is correct. The '-drop_below' option only exists in the lasheight tool.

[1] http://rapidlasso.com/2014/11/04/rasterizing-perfect-canopy-height-models-from-lidar/

Khosravipour, A., Skidmore, A.K., Isenburg, M., Wang, T.J., Hussin, Y.A., (2014). Generating pit-free Canopy Height Models from Airborne LiDAR. PE&RS = Photogrammetric Engineering and Remote Sensing 80, 863-872.

Khosravipour, A., Skidmore, A.K., Isenburg, M. and Wang, T.J. (2015). Development of an algorithm to generate pit-free Digital Surface Models from LiDAR. Proceedings of SilviLaser 2015, (available from http://www.isprs-geospatialweek2015.org/files/Proceedings_Silvilaser2015.pdf)

Regards,

Martin @rapidlasso

[Anahita Khosravipour]

Hi everybody,

I've already received couple of emails about the script file of the "pit-free CHM algorithm" which is described in our paper (Khosravipour et al. 2014). I would like to share it in the group then everybody can directly use it .

please find the attached *.txt file and save as it as *.bat file.  

Best regards

Anahita

[Satej Panditrao]

Hello,

Which programming language is this??

--

SATEJ PANDITRAO

Project Assistant

Advisory Services and Satellite Oceanography Group (ASG),

Indian National Center for Ocean Information Services (INCOIS),

Ministry of Earth Sciences, Govt. of India,

"Ocean Valley", Pragathi Nagar (BO),

 Nizampet (SO), Hyderabad- 500090, Telangana, India.

Email: satej.p...@gmail.com   (Personal)

          satejpa...@rocketmail.com   (Personal)

          pandi...@incois.gov.in   (Official)

Phone No: 040- 23886038 (Office)

Mobile no: +918978490266

Fax: +91-40-23895001

 

[Raja Ram Aryal]

Hello, 

If i would like to compare the pit-free CHM obtained as prescribed by Anahita text file with CHM with pit (Normal Pit). CHM with pit may be calculated like this or not?

blast2dem -i %INPUT_CANOPY_HEIGHT_POINTS%.%INPUT_FORMAT% ^

          -keep_first ^

          -step %STEP% ^

          -odir %TEMPORARY_CHM_DIR% -odix _00 -obil

Then i would further compared the DBH, Volume and tree height obtained from LiDAR metrics from both CHM using field inventory data. Is it possible or not?, Please suggest me. Thank you.

Best Regards

Raja Ram Aryal

On Wed, Jan 20, 2016 at 3:00 PM, Raja Ram Aryal <rajara...@gmail.com> wrote:

Hello, 

From those CHM how i can obtain LiDAR metrics to calute further tree structure tree height, DBH and volume. Thank you.

[Martin Isenburg]

Hello Raja,

Yes, a first-return CHM (with pits) at the highest resolution that the LiDAR supports be computed by simply interpolating all first returns into a triangulated surface (TIN) and then sampling this TIN at the centers of all raster cells. This is exactly the computation that happens when you call

las2dem -i height_normalized_lidar.laz ^

                 -keep_first ^

                 -step 0.5 ^

                 -o chm_with_pits.bil

blast2dem -i huge_height_normalized_lidar.laz ^

                 -keep_first ^

                 -step 0.5 ^

                 -o huge_chm_with_pits.bil

assuming the average pulse spacing is around 0.5 meters. Often you may want to add a -kill 1.5 so you do not triangulate across large voids in the data and create "fake" canopies by connecting the top of the vegetation to the other side of a river or a lake with huge slanted triangles. 

However, if you want to compute a CHM whose resolution is much *coarser* than that of the LiDAR pulse spacing (e.g. a 2 meter CHM for an 0.5 meter average pulse spacing) then often you will find that simple use of binning with lasgrid with lasgrid is sufficient to create a pit-free CHM as discussed here:

http://rapidlasso.com/2014/11/04/rasterizing-perfect-canopy-height-models-from-lidar/

Regards,

Martin @rapidlasso

[Irf]

Hi Anahita,

Thanks for sharing the script. I believe this is going to be useful for many. However, it would be worthwhile to mention that the data used in the "executing commands" are ground classified, height normalised and not splatted.


Regards,

Irf