1

I am trying to create a canopy height map from LIDAR point cloud data (LAS format) but I get some jagged edges/holes in my results (see photo). Is there a parameter I can tweak to fill these in? I tried adding the subcircle parameter to pitfree but still have the holes/edges. Here is the code I am using:

        library(lidR)
        library(terra)
        library(sf)
        print("Calculating CHM")
        las = readLAS(input)
        nonnoise <- filter_poi(las, Classification != 18) # remove high noise
        nlas <- normalize_height(nonnoise, knnidw())
        chm <- rasterize_canopy(nlas, res = 1, algorithm = pitfree(subcircle = 0.15), pkg = "terra")
        writeRaster(chm, sprintf('./out/%s_chm.tif', filename), overwrite=TRUE)

chm

chm2

Edit: I tweaked my code to the following to fix this issue:

chm <- rasterize_canopy(nlas, res = 1, algorithm = pitfree(max_edge = c(3, 3), subcircle = 1.4))
fill.na <- function(x, i=5) { if (is.na(x)[i]) { return(mean(x, na.rm = TRUE)) } else { return(x[i]) }}
w <- matrix(9,9,9)
filled <- terra::focal(chm, w, fun = fill.na)
smoothed <- terra::focal(filled, w = matrix(1,3,3), fun = mean, na.rm = TRUE)
4
  • What do you mean? What is the problem? I'm not even sur of what I'm seeing. I see some urban image under green vegetation that looks nice.
    – JRR
    Nov 29, 2023 at 2:33
  • @JRR Thank you for the reply. I edited my question to include a second photo of my CHM in grayscale with some more examples of the holes I get in my resulting CHM. I am wondering about methods I can try to remove the holes, or to avoid creating them in the first place
    – Max Ogden
    Nov 29, 2023 at 5:30
  • Ok I see, you are running a triangulation-based algorithm for CHM you should have a continuous CHM without any gaps. Yet you have gaps in... gaps. This is not expected indeed. Please share data so we can reproduce.
    – JRR
    Nov 29, 2023 at 10:01
  • @JRR Thank you. Here is a link to the raw .las file drive.google.com/file/d/1rYDgvrmrqqxy9VcYn31mMHahPw_tBmvR/… and accompanying documentation lariac-lacounty.hub.arcgis.com/pages/lariac4-documents-data
    – Max Ogden
    Nov 29, 2023 at 17:54

1 Answer 1

1

I ran the following. Everything look good except the geometric deformation due to steep slopes. Please provide a reproducible example.

library(lidR)
library(terra)
library(sf)

ctg = readLAScatalog("L4_6577_1894.las",  select = "xyzrc", filter = "-drop_class 18")
las = clip_rectangle(ctg, 6579000,1897000, 6579000+1000,1897000+1000)
nlas <- normalize_height(las, knnidw())
chm <- rasterize_canopy(nlas, res = 1, algorithm = pitfree(subcircle = 0.15))
plot(chm, col = height.colors(25))
4
  • Thank you, I ran this code in RStudio, saved the raster to tif writeRaster(chm, 'chmtest.tif'), and then opened in QGIS, zoomed in and am seeing this for example in the bottom left corner: imgur.com/a/18sXgtn, does that match what you are seeing? By the way, I am using sf_1.0-14 terra_1.7-55 lidR_4.0.4 with R version 4.3.2
    – Max Ogden
    Nov 29, 2023 at 22:49
  • The pit free algorithm is design to avoid pits. You have pits, which means that you must tune some parameters I guess. To be able to tune the parameters, you must understand the algorithm. To understand the algorithm, you must read the publication that is associated with it. I myself need to read it again and make some trial and error. By the way, I also wrote an illustrated demonstration of the algorithm here r-lidar.github.io/lidRbook/chm.html#pitfree
    – JRR
    Nov 29, 2023 at 23:36
  • thanks I'll play around some more and report results if I get any success!
    – Max Ogden
    Nov 30, 2023 at 4:26
  • I played around with the pitfree settings and decided these look best for my data: pitfree(max_edge = c(3, 3), subcircle = 1.4) and combined with the fill and smoothing steps from chapter 6.4 of your book I am now getting good looking output: imgur.com/a/sBjRYcw I updated my question with my working code.
    – Max Ogden
    Nov 30, 2023 at 6:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.