2

I am a first time user trying to help my community inventory trees over a threshold height using public Lidar data. My goal is to extract hulls and/or approximate GPS coordinates of each tree exceeding the target height.

I followed along with the individual tree segmentation examples (1,2) but in my results trees are generally segmented into 4-8 parts based on the colorization in the image (included).

I first thought this had to do with the watershed parameters, but I didn't hit on the right combo in the workflow. Later I tried the li2012 algorithm and my results were no better.

Can someone suggest either the key parameters I should be changing, or a different workflow I should be following?

print(Sys.time())
library(lidR)

las = readLAS("/Users/brianvanvoorst/Desktop/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.las")
las = classify_ground(las, csf())
print ("Classify done")
las = normalize_height(las, tin())
print ("Normalize done")
algo = pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
print ("Pitfree done")
chm  = grid_canopy(las, 0.5, algo)

plot(chm, col = height.colors(50))
# smoothing post-process (e.g. two pass, 3x3 median convolution)
ker = matrix(1,3,3)
chm = focal(chm, w = ker, fun = median)
chm = focal(chm, w = ker, fun = median)

plot(chm, col = height.colors(50)) # check the image
algo = watershed(chm, th = 4)
las  = segment_trees(las, algo)

# remove points that are not assigned to a tree
trees = filter_poi(las, !is.na(treeID))

plot(trees, color = "treeID", colorPalette = pastel.colors(100))
print(Sys.time())

Inset of segmentation over CHM, "easy" trees are fractioned into parts


Okay here is my latest source and result

require(lidR)
require(rlas)
require(rgdal)
require(tictoc)

las <- readLAS("/Users/brianvanvoorst/Desktop/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.las", filter="-keep_class 1L")

#dtm <- grid_terrain(las, algorithm = knnidw(k = 8, p = 2))
# Error: No ground points found. Impossible to compute a DTM.
#las_normalized <- normalize_height(las, dtm)

# Create a filter to remove points above 95th percentile of height
filter_noise = function(las, sensitivity)
{
  p95 <- grid_metrics(las, ~quantile(Z, probs = 0.95), 10)
  las <- merge_spatial(las, p95, "p95")
  las <- filter_poi(las, Z < p95*sensitivity)
  las$p95 <- NULL
  return(las)
}

las_denoised <- filter_noise(las, sensitivity = 1.2)

chm <- grid_canopy(las_denoised, 0.5, pitfree(c(0,2,5,10,15), c(3,1.5), subcircle = 0.2))

plot_dtm3d(chm)

ker <- matrix(1,5,5)
chm_s <- focal(chm, w = ker, fun = median)

algo <- watershed(chm_s, th = 4)
las_watershed  <- segment_trees(las_denoised, algo)

# remove points that are not assigned to a tree
trees <- filter_poi(las_watershed, !is.na(treeID))

# View the results
plot(trees, color = "treeID", colorPalette = pastel.colors(100))

enter image description here

5
  • 1
    Can you edit your post to include your code? Is this the tutorial you went through?: github.com/Jean-Romain/lidR/wiki/…
    – Aaron
    Oct 4, 2020 at 3:49
  • Sorry that took so long. I have edited to include the code. I didn't go through the example you linked to, I had gone through another example - (part 1). Oct 7, 2020 at 15:59
  • Please try the tutorial that I linked to and see if that solves your issue.
    – Aaron
    Oct 7, 2020 at 16:12
  • Ok. Source included. I think the results are somewhat improved but still not distinct trees where I would have hoped/expected them. I included a shot of the same area so you could see for comparison. I am wondering if a different kernel would have helped, but I am not sure if I should go larger or smaller. Also not sure about parameters to grid_canopy should be different. Also of note: 1) filtering points was interesting and I think helpful, but unsure. I ended up with different classes. 2) grid_terrain didn't work due to no ground points - is this my problem? Oct 9, 2020 at 1:04
  • You need to classify the ground points if they are not already classified. You can use lasground in lidR for that.
    – Aaron
    Oct 9, 2020 at 1:48

1 Answer 1

1

So I eventually found the data online:

  1. You are working in an urban context. lidR's algorithms are designed to work in a forest context. A tree is a tree, but a point is also a point. You will unavoidably segment building as trees because there is no way to make the distinction between a tree and a building. The point cloud must be classified upstream if you want to get a chance to filter out the buildings.

  2. Your point cloud is in feet. You are providing algorithm parameters in meters. No chance to get a good output. The csf performs unreasonably slow because it believes it is processing a 6.25 km2 tile while it is actually a 0.7 km2 tile and your CHM has a 15 cm resolution for a point cloud with a density of 3 points/m2 (roughly). All the parameters are irrelevant.

After using feet-based parameters + parameters more carefully chosen for this dataset the CHM looks nicer and you can start doing something with it. However building and other human made structure such as wire conductors will still be segmented as trees.

las = readLAS("/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.laz", filter = "-drop_y_above 523580 -drop_x_below 19381500")
las = normalize_height(las, tin())
thresholds = round(c(0,5,10,15,20,25)/0.3048,0) # The highest point is ~80 feet ~= 25 m
algo = pitfree(thresholds = thresholds, max_edge = c(0, 2), subcircle = 0.2/0.3048)
chm = grid_canopy(las, 2, algo)
plot(chm, col = height.colors(50))

Before

enter image description here

After

enter image description here

2
  • Thank you for this help! It is really appreciated. Since this is a one time community project I am happy to live with the buildings. This has been a very interesting learning experience. I will take what you have given me and see what progress I can make towards my goals. MUCH APPRECIATED. Oct 14, 2020 at 1:56
  • Please mark the question as anwsered
    – JRR
    Oct 14, 2020 at 6:43

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.