I have a single band DEM as a GeoTiff. I would like to replace values above or below a certain threshold with 0. Can I do that with GDAL?


I'm not sure if you can do this with any of the gdal cli tools, but I wrote something in python which accomplishes it:

from osgeo import gdal
from osgeo.gdalconst import GDT_Float32

import sys
import numpy as np

def fix_dem_nodata(raster_input, raster_output, nodata=0, threshold=-900):

        in_data, out_data = None, None

        # open input raster
        in_data = gdal.Open(raster_input)
        if in_data is None:
            print 'Unable to open %s' % raster_input
            return None

        # read in data from first band of input raster
        band1 = in_data.GetRasterBand(1)
        rows = in_data.RasterYSize
        cols = in_data.RasterXSize
        vals = band1.ReadAsArray(0, 0, cols, rows)

        # create single-band float32 output raster
        driver = in_data.GetDriver()
        out_data = driver.Create(raster_output, cols, rows, 1, GDT_Float32)
        if out_data is None:
            print 'Could not create output file %s' % raster_output
            return None

        # set values below nodata threshold to nodata
        dem_data = np.array(vals)
        dem_data[dem_data < threshold] = nodata

        # write the data to output file
        out_band = out_data.GetRasterBand(1)
        out_band.WriteArray(dem_data, 0, 0)

        # flush data to disk, set the NoData value and calculate stats

        # georeference the image and set the projection

        return raster_output

    except Exception as e:
        print "Failure to set nodata values on %s: %s" % raster_input, repr(e)
        return None

        del in_data
        del out_data

if __name__ == "__main__":
    infile = sys.argv[1]
    outfile = sys.argv[2]
    print fix_dem_nodata(infile, outfile)

I set the threshold to -900 as it seemed a reasonable limit for realistic dems, but you can tune for your dataset. N.b. you may want to use a value other than 0 for your nodata, as 0 is a realistic elevation.

| improve this answer | |
  • Brilliant saved me a couple of hours. – HeikkiVesanto Feb 13 '18 at 20:22

You can do something similar using gdal_calc.py, e.g.:

gdal_calc.py -A dtm.tif --calc='((A>=100)*(A<=200))*A+((A<100)*0)+((A>200)*0)' --outfile=dtm_reclass.tif --NoDataValue=-32767

This calc expression would:

  1. Assign a value of 0 to all pixel values less than 100 ((A<100)*0)
  2. Assign a value of 0 to all pixel values greater than 200 ((A>200)*0)
  3. Keep existing pixel values greater than or equal to 100 and less than or equal to 200 ((A>=100)*(A<=200))*A
| improve this answer | |

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.