I have a .las file and want to show the outline of the big rocks located in this area.

So far I have done the following:

  1. Obtain a ground model from the .las file including the rocks.
  2. Calculate the slopes in the ground model.
  3. Only show slopes greater than 30 degrees.

Step 2 and 3. are fine. My problem is to obtain a reasonable ground model in step 1. (Step 3 unfortunately removes some of the smaller rocks, but this is not very important).

The .las file is downloaded from hoydedata.no. It is available here. In the linked map there are two big rocks showing up in the DSM Hillshade layer. One rock next to the road, and one next to a small lake. The rock next to the lake does not show up in the DTM Slope degrees layer.

Further information:
The points in the .las files are classified as 1, 2, 7. It is somewhat arbitrary if big rocks are assigned class 1 or 2.

A rock is big if it is greater than 2 meters in at least 2 directions (and usually in all 3 directions). However, a method which can outline smaller rocks as well would be great.

The rocks are typically located in forests, swampland or mountains. There are usually no buildings around. There might be cliffs around, but these I want to show as well. Hence, an alternative could be to identify and remove the vegetation.

Are there any tools to remove vegetation points from .las files?

I have tried various tools for step 1:
pdal works when I increase the slope parameter slightly, but I worry that there are some trees showing up:

pdal ground -i /data/in.las -o /data/out.las --slope 0.7 --denoise

lidR has an impressive tree segmentation algorithm which color the rocks correctly, see 2. However, I was unable to extract the rocks or remove the vegetation. It looks like the rocks are classified as trees.

lastools works reasonably well with the following parametres (suggested by the data providers):

lasground.exe -i *.laz -ultra_fine -ignore_class 7 -step 5  -spike 1 -offset 2.5  -odir ground -olaz -cores 5
las2dem.exe -i ground\*.laz -slope -keep_class 2 -step 0.25 -opng -odir dtmslope

Unfortunately I do not have access to a licensed version of lastools, and I only have limited access to a Windows environment for running lastools.

Preferably I would like to use open source tools.

I am completely new to GIS and LiDAR.

How do I tweak the parameters or should I use different tooling for step 1?

For step 2 I currently use lidR and the following code:

las <- readLAS("out.las", select="xyzc")
las <- lasfilter(las, Classification == 2)
gnd <- grid_terrain(las, res = 1, tin())
sl <- terrain(gnd, opt="slope")
png("tmp.png", width = 1920, height=1080)

This produces the picture below from the output of pdal. In the picture you can see the outline of two rocks at the bottom center next to the road and the small lake.

Slopes from ground model.

  • Your question is "how to segment rocks?". It is a too broad question. I have never heard about that algorithm. You should better ask a single and focused question such as "how to extract my segmented rocks with lidR" or "how to get a highly details DTM that includes rocks with pdal?" for examples. – JRR Oct 15 '19 at 16:12
  • @JRR "how to get a highly details DTM that includes rocks with pdal" (and not trees) would solve my problem. Is there a way to achieve this with pdal? – YellowBird Oct 15 '19 at 20:24
  • Have you looked at intensity values of points associated with rocks? Any significance there? – Aaron Oct 15 '19 at 21:01
  • My initial thought is to train a feature extraction model based on derived LiDAR products. How much training data can you generate? – Aaron Oct 15 '19 at 21:08
  • @Aaron Most of the rock points have intensity < 1200, and > 800, while most of the points which make up trees have intensity < 600. I will see if I can use this to remove the trees. – YellowBird Oct 16 '19 at 9:22

I tried naively to use nearest neighbors to separate the ground and the trees. I iteratively set a point to be a tree, ground, or other type of point. A point is a tree or ground point if a percentage of the k=10 nearest neighbors are tree or ground points. The percentage decreases with each iteration. The initial guess of trees are points with intensity < 600 (I also tried points where the convex hull of the nearest neighbors has volume > 1e6, but this produces worse results).

import numpy as np
import pylas
from sklearn.neighbors import NearestNeighbors
import sys

def summary(i, istree):
    print("{0}: trees: {1}, ground: {2}, other: {3}"
                  (istree == 1).sum() / len(istree),
                  (istree == -1).sum() / len(istree),
                  (istree == 0).sum() / len(istree)))

def main(fname1, fname2):
    with pylas.open(fname1) as f:
        outf = f.read()
    pts = outf.points
    data = np.asarray([pts[i] for i in ("X",
    coords = data[:, :3]

    Nnbs = 10
    nbrs = NearestNeighbors(n_neighbors=Nnbs,
    distances, indices = nbrs.kneighbors(coords)

    tree_pts = np.logical_or(pts["intensity"] < 600, False)
    ground_pts = True ^ tree_pts

    istree = tree_pts * 1 + -1 * ground_pts

    for i in range(3):
        istree = solve(istree, indices,
        summary(i, istree)

    niter = 10
    for i in range(niter):
        istree = solve(istree, indices,
                       bound=0.2 * (niter - 1 - i) / (niter - 1))
        summary(i, istree)

    classification = np.zeros(len(istree), dtype=np.uint8)
    classification += np.uint8(3) * (istree == 1)
    classification += np.uint8(2) * (istree == -1)
    classification += np.uint8(1) * (istree == 0)

    outf.points["raw_classification"] = classification

def solve(istree, indices, bound):
    n_istree = np.zeros(shape=(istree.shape[0],), dtype=np.int8)

    for pt, inds in enumerate(indices):
        s = istree[inds].sum() / len(inds)
        if s > bound:
            n_istree[pt] = 1
        elif s < -bound:
            n_istree[pt] = -1
            n_istree[pt] = 0

    return n_istree

if __name__ == "__main__":
    main(sys.argv[1], sys.argv[2])

Unfortunately this code produces worse results than pdal.

The .las file for testing is available here.

Combining the above code with the individual tree segmentation of lidR gives good results.


las = readLAS("test.las")
las = lasground(las, csf())
algo = pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
chm  = grid_canopy(las, 0.5, algo)
algo = watershed(chm, th = 4)
las  = lastrees(las, algo)
writeLAS(las, "tmp2.las")

For each point classified to be part of a tree by lidR I vote if this is a tree or ground by taking a weighted average of the result of the nearest neighbor python code. This is done in the following python code (tmp.las is the output of the nearest neighbor code).

import numpy as np
import pylas
from collections import defaultdict

with pylas.open("tmp.las") as f:
    x = f.read()

with pylas.open("tmp2.las") as f:
    y = f.read()

xx = np.array([x[v] for v in ("X", "Y", "Z", "intensity")]).transpose().copy()
yy = np.array([y[v] for v in ("X", "Y", "Z", "intensity")]).transpose().copy()

valid_points = []
ii = 0
for i in range(len(x.points)):
    if np.count_nonzero(xx[i] - yy[ii]) == 0:
        ii += 1
        if ii > len(yy): break

assert len(valid_points) == len(yy)

x.points = np.array(valid_points, dtype=x.points.dtype)

treeID = y["treeID"]
classification = x["classification"]
ground = y["classification"] == 2

istree = (-0.25 * (classification == 2) +
          1 * (classification == 3) +
          0 * (classification == 1))

d = defaultdict(int)
for i in range(len(x.points)):
    d[treeID[i]] += istree[i]

for i in range(len(y.points)):
    if d[treeID[i]] > 0 and not ground[i]:
        classification[i] = 3
        classification[i] = 2

y["classification"] = classification

Still there are some trees which are classified as ground. I suspect the reason is that the treeID of lidR is an 8 bit integer and there seems to be more than 256 trees.

| improve this answer | |
  • 8 bits integer does not exist in R – JRR Oct 18 '19 at 12:40
  • @JRR in rdrr.io/github/Jean-Romain/lidR/man/lastrees.html it is stated that 1 byte is used to store treeID. – YellowBird Oct 18 '19 at 14:49
  • 1
    It is said that it is added as an extra byte attribute. See cran.r-project.org/web/packages/lidR/vignettes/… section Extra attributes and extra bytes in a LAS object. – JRR Oct 18 '19 at 14:58
  • Thanks! I thought treeID was stored as an integer (looking at the output with pylas confused me), but apparently treeID is a double. So my guess for why some trees are classified as ground is wrong. The reason must then be that too large a proportion of the points classified as a particular tree are actually ground points. – YellowBird Oct 18 '19 at 15:56
  • If you want details about treeID and so on please open a new question. But treeID should be an integer (on 32 bits) – JRR Oct 18 '19 at 16:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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