How to replace pixel values in a single band DEM?

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.