Convert DSM to DTM and nDSM

Question

I have some DSM GeoTIFF files and I want to convert them to DTM and then obtain nDSM by subtracting DTM from DSM. Is there any way to do it?

I want to apply this to several different DSMs so it's best if I find an automatic code or bash tool.

Answer 1

So I found the answer myself.

This is the pipeline:

1. Convert DSM to LAS.
2. Classify LAS ground points.
3. Create DTM from classified LAS file.

There is LAStools in which there is a command txt2las. So if I have a txt file which contains XYZ values in each row of the file, then I can use this command to obtain LAS.

So first I have to drag XYZ values from my Geotiff image to a text file. To do this I used answers from this question. I extracted x and y values from the method there and z values from the DSM itself.

def convert_dsm_to_txt(out_folder_path, input_dsm_path):
# Read raster
with rasterio.open(input_dsm_path) as r:
    T0 = r.transform  # upper-left pixel corner affine transform
    dsm = r.read()  # pixel values

# All rows and columns
cols, rows = np.meshgrid(np.arange(dsm.shape[2]), np.arange(dsm.shape[1]))

# Get affine transform for pixel centres
T1 = T0 * Affine.translation(0.5, 0.5)
# Function to convert pixel row/column index (from 0) to easting/northing at centre
rc2en = lambda r, c: (c, r) * T1

# All eastings and northings (there is probably a faster way to do this)
eastings, northings = np.vectorize(rc2en, otypes=[float, float])(rows, cols) 

# only get non NaN values. -32767 in my file represent NaN cells so I don't want them
dsm = dsm[0, ...]
x = eastings[dsm != -32767]
y = northings[dsm != -32767]
z = dsm[dsm != -32767]

d = np.stack((x, y, z)).T
del x, y, z, dsm
np.savetxt(os.path.join(out_folder_path, os.path.basename(dsm_path) + '.txt'), d)

Now that I have a text file containing xyz values, I can run the command I specified earlier.

txt2las -parse xyz -i input.txt -o output.las

here we have a LAS file. To convert it to DTM I used lidR library in R. I used the function classify_ground to classify the ground points and then I used grid_terrain to obtain DTM.

las <- readLAS(las_path)
las <- classify_ground(las, csf())
dsm <- raster(x = dsm_path)
dtm <- grid_terrain(las, algorithm = tin(), res = dsm)
writeRaster(dtm, dtm_out, overwrite=TRUE)

At last subtracted DTM from DSM to get nDSM. here is the result.

The resulted nDSM is exactly what I wanted.

I hope this answer helps the next person with same question.

Answer 2

GeoTiffs are a derivitive product. If you have the original LiDAR Point Clouds and the points are classified, then you can filter the ground points using software like LasTools or CloudCompare. But if all you have are the GeoTiffs, I don't think there's any way.