Classify DEM pixels

Question

I have a differential surface model dataset. That was made: DSM subtracted from DTM, so for example trees have negative height. It is a raster (TIFF) dataset, Type=Float32.

I want to create a shapefile that will have three different values:

• 0: error (all the pixels in the original dataset with values > 0)
• 1: OK (all the pixels in the original dataset with values < 0)
• 2: nodata (all the pixels in the original dataset with values = -10000000000)

I assume first thing to create a raster with pixel values 0, 1, 2 only and than vectorize it. Probably the vectorization will not be a problem, using gdal_polygonize.py but creating the raster is quite a challenge for me, although I was able to achieved some minor success using gdal_calc.

I'd like to do this using GDAL, do you have an idea?

Answer 1

If you can use R, it's raster package supports intuitive raster algebra that will make your task very simple. Try this:

library(raster)

rred <- brick(system.file("external/rlogo.grd", package="raster"))[['red']] ## MWE, use your data
values(rred) <- values(rred) - mean(values(rred))

rbad <- rred > 0
rgood <- rred <= 0
rmiss <- rred == -100 ## please use your true value of -10e9, this is only for MWE

writeRaster(rgood, filename='good.tif')
writeRaster(rbad, filename='bad.tif')
## etc

now polygonize them with GDAL, if you wish. Alternatively install and load gdalUtils and do it all inside R. The raster package can do polygonization too, but it's insufferably slow.

Answer 2

You can use the gdal python bindings to read the array and then numpy to change the pixel values. Then you can use gdal again to overwrite your original raster.

import gdal

fn = 'path_to_your_raster.tif'
ds = gdal.Open(fn, 1)   # open TIFF file in writing mode
arr = ds.ReadAsArray()  # read the data as a numpy array

# specify your conditions
con1 = (arr > 0)
con2 = (arr < 0)
con3 = (arr == -10000000000)

# change values according to the conditions
arr[con1] = 0
arr[con2] = 1
arr[con3] = 2

# overwrite the TIFF data
band = ds.GetRasterBand(1)
band.WriteArray(arr)

# save file and close
ds.FlushCache()
del ds, band