5

I have two rasters: Raster A and B. Raster A has has only values of 0. The white area in Raster B are No Data values. I want to merge them to one raster (Raster C). So basically I want all pixels in Raster B with a No Data value to have a value of 0.

I tried to do this with the raster calculator in QGIS by adding Raster A to Raster B, but it only returns again Raster B.

I am looking for solutions with QGIS, GDAL, or Grass GIS.

enter image description here

  • 3
    You would benefit by reviewing the GRASS map calculator expressions that are available, paying special attention to how they support null data. They provide many solutions, such as combining isnull with the conditional operator ?:. – whuber Feb 5 '14 at 18:55
4

One solution may be this batch script:

set nodatavalue=-99999
gdalwarp -srcnodata None rasterB_with_NoData.tif rasterB_without_NoData.tif
python gdal_calc.py -A rasterB_without_NoData.tif --outfile=rasterB_updated_with_NoData.tif --calc="(A<>%nodatavalue%)*A"
gdalwarp -srcnodata None rasterB_updated_with_NoData.tif rasterB_updated_without_NoData.tif

The trick is -srcnodata None offered by gdalwarp, which is very useful in order to ignore intrinsic nodata settings on the source dataset (from the doc). In detail, we need to apply it two times: firstly because we need to substitute the NoData value with 0 and so we have to include it in the calculation, secondly because gdal_calc.py introduces NoData value again.

1

I ended up writing a GDAL python script, that replaces the NoDataValue with a new value (in my case 0):

import gdal, ogr, osr, os
import numpy as np

def raster2array(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.ReadAsArray()

def getNoDataValue(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.GetNoDataValue()

def array2raster(rasterfn,newRasterfn,array):
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3] 
    pixelWidth = geotransform[1] 
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


rasterfn = 'Slope.tif'
newValue = 0
newRasterfn = 'SlopeNew.tif'

# Convert Raster to array
rasterArray = raster2array(rasterfn)

# Get no data value of array
noDataValue = getNoDataValue(rasterfn)

# Updata no data value in array with new value
rasterArray[rasterArray == noDataValue] = newValue

# Write updated array to new raster
array2raster(rasterfn,newRasterfn,rasterArray)

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.