Replacing nodata values by a constant using GDAL CLI tools

Question

I have a raster DEM file containing 1 band: f='DEM.tif'
It has some NODATA values which I want to replace by a true, actual constant value (i.e. no more nodata stuff), let's say, zero.

Here are all the things I've tested:

1. gdal_translate; not working:

gdal_translate -of GTiff -a_nodata 0 ${f} ${f%.*}_cleaned.tif

well... it actually replaces the nodata value of 3.399999^38 to a nodata value of 0, but it doesn't replace it by a true 0 (which would be an actual zero value, no a nodata value).

2. I also tried the accepted answer here using gdalbuildvrt: Redefining nodata value into zero in QGIS? but with no success. Nodata values still show as ... well, nodata values.

3. Also, there is the gdal_fillnodata tool, which seemed promising: https://gdal.org/programs/gdal_fillnodata.html but it actually interpolates values from the closest borders, which doesn't make sense in my case.

4. Finally, gdalwarp was my last candidate, but it is not working better:

gdalwarp -srcnodata 3.39999995214436425e+38 -dstnodata 0 ${f} ${f%.*}_cleaned.tif

because, here again, it replaces the nodata value by another nodata value, which is not what I would like to achieve.

I've done that successfully in no time with QGIS using the "Fill NoData cells" on my local machine as a test, but I must work on my full dataset on a headless server, without QGIS installed.

So, do you know how to fill nodata values by a constant using any of the GDAL CLI (bash) based tools? I hope I simply missed something obvious...

Edits

5. Removed as it was a malformed version of @user2856's working answer.

6. I forgot to mention this when I wrote the post, but I also tried @Mike T's first part of his answer:

gdal_edit.py -unsetnodata ${f}

(Disclaimer: this is a non reversible operation as it's actually working 'inplace', so backup your raster before running this command!)

Here, on the contrary, gdalinfo doesn't report NoData anymore. But do not rejoice too quickly, as a loading into QGIS shows:

The second part of Mike T's answer:

gdal_calc.py -A ${f}.tif --outfile=result.tif --calc="numpy.where(A>=3e34, 0, A)"

is actually giving a result.tif raster with:

gdalinfo result.tif | grep No
  NoData Value=3.40282346600000016e+38

That's funny? Probably some weird numerical issue under the hood.
But the most funny part is that if I first run gdal_edit.py -unsetnodata and then the latter gdal_calc.py command, the resulting raster has so nodata again (heck, where do these actually come from?)!

---

General info on my setup just in case:

# gdalinfo --version
GDAL 3.3.0dev-33cf0e31a992be112b3091f012368d15605ed51d, released 2021/03/17

from osgeo/gdal:ubuntu-small-latest Docker image available at https://hub.docker.com/r/osgeo/gdal

Answer 1

One (slightly hacky) way of doing it is to use an intermediate VRT with the gdalbuildvrt command and specifying a new NoData value with the srcnodata and vrtnodata arguments. This does actually change the underlying values.

First, get your actual nodata value from gdalinfo:

$ gdalinfo input.tif | grep No
NoData Value=3.39999995214436425e+38

Assign it to a bash variable: na=3.39999995214436425e+38

Or if you have jq installed (in ubuntu/debian based Linux distros sudo apt install jq):

na=$(gdalinfo input.tif -json | jq '.bands[0].noDataValue')

And use it in the following examples:

gdalbuildvrt -srcnodata ${na} -vrtnodata 0 output.vrt input.tif

However, it also sets those values to NoData in the output, so you need to specify a different value to be considered as NoData when converting to your final GeoTIFF.

gdal_translate -a_nodata ${na} output.vrt output.tif  # original 3.39999995214436425e+38 values from input.tif are now 0

Chaining it together with a pipe so you don't need to clean up any intermediate VRT files:

gdalbuildvrt -srcnodata ${na} -vrtnodata 0 /vsistdout/ input.tif | gdal_translate -a_nodata -9999 /vsistdin/ output.tif

Answer 2

Here are some untested ideas:

1. Use gdal_edit.py to remove the NoData metadata from the file. This is an in-place edit, so make sure you have a back-up, or use a copy of the original file. E.g.:

   gdal_edit.py -unsetnodata file.tif

2. Use gdal_calc.py to modify the values, by replacing one value with another. Numpy's where function does this. For example to replace a previous nodata 3.399999e+38 with 1.1:

   gdal_calc.py -A file.tif --outfile=result.tif --calc="numpy.where(A>=3e38, 1.1, A)"

The condition A>=3e38 will avoid floating point issues while trying to pick exact values.

Answer 3

I would recommend the following which works for me:

1. replace your_original_nodata_value (e.g. 0) with a_new_nodata_value (e.g., -9999):

gdal_calc.py -A input.tif --outfile=output.tif --calc='(A==0)*(-9999)+(A!=0)*A'

2. set the new_nodata_value as nodata:

gdal_translate -a_nodata -9999 output.tif output_2.tif

Then output_2.tif should be your desired output.

Answer 4

In the meantime, I build a small Python function depending on rasterio that may be helpful for people having access to Python. But I consider this only as a quick'n'dirty alternative:

import os, glob, rasterio

folder = '/path/to/your/rasters/'
files_list = glob.glob(os.path.join(folder,'*.tif'))

def replace_nodata_by(input_file, replacement_value=0):
    epsg = 31983
    outfile = input_file.split('.')[0]+'_fixed.tif'
    with rasterio.open(input_file, mode='r', driver='GTiff',
                       count=None, crs='EPSG:'+str(epsg), transform=None,
                       dtype='float32', nodata=None) as source_dataset:
        data_array = source_dataset.read()
        data_array[data_array > 1e5] = replacement_value # kick out nodata values
        # This is a bonus step in my case to fill under ocean level's 
        # values with zeros as well:
        data_array[data_array < 0] = replacement_value
        output_metadata = source_dataset.meta
   
    output_metadata.update({"driver": "GTiff", "nodata": -9999})
    with rasterio.open(outfile, "w", **output_metadata) as destination_dataset:
        destination_dataset.write(data_array.astype(rasterio.float32))
    
    return True
        
for file in files_list:
    print(f"Processing {file}...")
    replace_nodata_by(file, replacement_value=0)

gdalinfo will tell there still is a nodata value of -9999 now, but if you carefully look at the code, I've set up this new -9999 value just to be able to differentiate it from the original nodata value in case of any kind of failure/bug. But normally, there will be no such values in the output raster files.

Answer 5

An alternative is to use gdal_calc and numpy.nan_to_num:

gdal_calc.py -A input.tif --outfile=output.tif \
    --calc="numpy.nan_to_num(A)"

Variations on this should let you specify the value to replace NaNs:

gdal_calc.py -A input.tif --outfile=output.tif \
    --calc="numpy.nan_to_num(A, nan=123)"

Answer 6

I was looking for an answer to this question, and went back to the documentation of gdal_calc.py here: https://gdal.org/programs/gdal_calc.html

If you are using version 3.3 or newer, you can combine the following to get the desired results:

• use numpy.where to replace the old NoDataValue with new value
• option --hideNoData to ignore the input band NoDataValue
• option --NoDataValue='none' will result in output with no NoDataValue

Your command would look like this:

gdal_calc.py -A DEM.tif --calc="numpy.where(A>=3e38, 0, A)" \
  --outfile DEM_clean.tif --hideNoData --NoDataValue='none'