# How to replace pixel values in a single band DEM?

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?

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`

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):
try:

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
out_band.FlushCache()
out_band.SetNoDataValue(nodata)

# georeference the image and set the projection
out_data.SetGeoTransform(in_data.GetGeoTransform())
out_data.SetProjection(in_data.GetProjection())

return raster_output

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

finally:
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.

• Brilliant saved me a couple of hours. Feb 13, 2018 at 20:22