Extracting (big) rocks from LiDAR data

Question

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:

require(lidR)
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)
plot(sl)
dev.off()

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.

.

Answer 1

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}"
          .format(i,
                  (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",
                                        "Y",
                                        "Z")]).transpose().copy()
    coords = data[:, :3]

    Nnbs = 10
    nbrs = NearestNeighbors(n_neighbors=Nnbs,
                            algorithm='ball_tree').fit(coords)
    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,
                       bound=0.2)
        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
    outf.write(fname2)


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
        else:
            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.

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

require(lidR)

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:
        valid_points.append(x.points[i])
        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
    else:
        classification[i] = 2

y["classification"] = classification
y.write("tmp3.las")

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.