Python and R – Convert DSM to DTM and nDSM

dempythonr

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.

Best Answer

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. results

The resulted nDSM is exactly what I wanted.

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

Related Question