2

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.

3
3

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.

3
  • what do you mean by classified? I don't have Point Clouds but I may be able to generate a point cloud from pixel centers of my DSM. if I do that how can I obtain the DTM?
    – Parsa
    Jul 31 at 13:10
  • Pixel centers from your GeoTiff DSM will not give you what you want. If you want a DTM you'll have to find different data.
    – Pointdump
    Jul 31 at 16:47
  • They exactly give what i want.
    – Parsa
    Jul 31 at 17:03
3

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.

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.