r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
596 stars 131 forks source link

Clarification on the dsmtin algorithm #583

Closed floresteiro closed 2 years ago

floresteiro commented 2 years ago

Hi Jean-Romain,

I tried to usedsmtin to generate a DSM similar to the one I can get when I use only first returns in lastools' las2dem (

this tool reads LIDAR points from the LAS/LAZ format, triangulates lidar points temporarily into a TIN, and then rasters the TIN onto a DEM

), because it seemed to work like (similarly) las2dem works. However, the output of dsmtin is not as detailed as that of las2dem's and looks more like the outputs of algorithms that use the highest point found in each pixel, such as p2r or the lasgrid tool (with the -highest setting) on lastools.

I am particularly interested in a DSM that does not remove/ignore lower first returns and, therefore, provides a more detailed crown (rougher canopy). In this regard, If possible, I would like to know if I can get a DSM using dsmtinthat reflects the first return penetration in the canopy, such as the one I got with las2dem.

In the figure below, the left side shows the output of dsmtin and the right side shows the output of las2dem. image The mapswipe tool, on QGIS, was used to divided the same area in two. Therefore, the picture shows adjacent areas.

Jean-Romain commented 2 years ago

According to the doc of dsmtin

This function is made to be used in rasterize_canopy. It implements an algorithm for digital surface model computation using a Delaunay triangulation of first returns with a linear interpolation within each triangle.

There is no mention of first returns in las2dem. Did you called las2dem with a first return filter? Or did you include all returns?

floresteiro commented 2 years ago

Yes, I called it with a first return filter. So, I was expecting that they would produce similar outputs.

Jean-Romain commented 2 years ago

Can you run the same on the example file included with lidR Megaplot.laz. Show me the lidR output + the R source code + the lastools output + the lastools command. Remember that lastools is close source. We can't know exactly what it does. Maybe the rasterization step is different. But I'll try to provide an explanation and an option in lidR to reproduce more closely lastools output.

floresteiro commented 2 years ago

lidR R source code:

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)
dsm = rasterize_canopy(las, res=0.5, algorithm = dsmtin())
writeRaster(dsm,"D:/Megaplot_50cm_dsmtin.tif")

lastoolscommand:

las2dem -i D:\Megaplot.laz -o d:\Megaplot_dsm_50cm_las2dem.tif -otif^ 
        -first_only ^
        -step 0.5 ^
        -elevation 

In the picture, the left side shows the output of dsmtin and the right side shows the output oflas2dem. image The mapswipe tool, on QGIS, was used to divided the same area in two. Therefore, the picture shows adjacent areas.

Jean-Romain commented 2 years ago

I'm not able to answer. One step of dsmtin consists in keeping only the highest point per pixel before the triangulation to reduce computation time. This is not documented and I will add and option to include all points and I will update the doc. Yet if I produce a CHM skipping this step I still do not have a las2dem-like CHM.

dsm1 = rasterize_canopy(las, res=0.5, algorithm = dsmtin(highest = TRUE))
dsm2 = rasterize_canopy(las, res=0.5, algorithm = dsmtin(highest = FALSE))

Rplot

I can hardly see were those las2dem pits are coming from. Look at the 1-m resolution point-to-raster CHM. There is no extra pits

dsm3 = rasterize_canopy(las, res=1, algorithm = p2r())

Rplot