[GIS] -9999 (no data value) becomes 0 when writing array to GDAL memory file

gdalmemory-layerpythonpython-2.7

I am trying to reproject and resample (1000m ->> 30m) a raster image (shape 238X406) using gdal.ReprojectImage() function in Python. The input file has -9999 cells as no data value. The result is an array (shape 10834×15055). The data type is float32.

When I write the result to a geotiff file, everything is expected.No data value is set to -9999 and output array has -9999 cells.
However, when I write the result to a gdal memory file to save some time (25 seconds shorter processing time, down from 105 seconds to 80 seconds), this time all -9999s (no data value) become 0 (zero) in the output array of memory file. Both results are exactly same, but GeoTIFF file has -9999, memory file has 0. Even though memory file has no data value of -9999, but the output array is initialized to 0.0 instead of -9999.

I am using the same code to produce the results and the only difference is that I call memory driver when I want to write result to memory file
(driver=gdal.GetDriverByName("MEM")), and GeoTIFF driver when I want to write result to GeoTIFF file (driver=gdal.GetDriverByName("GTiff"))

My gdal version is '1100000' from gdal.VersionInfo(). My operating system us ubuntu 12.04 LTS "precise".

# Create a file
outFileRead=driver.Create(outFilePath,X,Y,1,dataType,options)
print inFileRead.GetRasterBand(1).GetNoDataValue()
print inFileRead.GetRasterBand(1).ReadAsArray()
print inFileRead.GetRasterBand(1).ReadAsArray().shape
# Reproject
gdal.ReprojectImage(inFileRead,outFileRead,inProjection,outProjection,reSamplingType)
print outFileRead.GetRasterBand(1).GetNoDataValue()
print outFileRead.GetRasterBand(1).ReadAsArray()
print outFileRead.GetRasterBand(1).ReadAsArray().shape

the result of prints

-9999.0

[[-9999. -9999. -9999. …, -9999. -9999. -9999.]

[-9999. -9999. -9999. …, -9999. -9999. -9999.]

[-9999. -9999. -9999. …, -9999. -9999. -9999.]

…,

[-9999. -9999. -9999. …, -9999. -9999. -9999.]

[-9999. -9999. -9999. …, -9999. -9999. -9999.]

[-9999. -9999. -9999. …, -9999. -9999. -9999.]]

(406, 238)

-9999.0

[[ 0. 0. 0. …, 0. 0. 0.]

[ 0. 0. 0. …, 0. 0. 0.]

[ 0. 0. 0. …, 0. 0. 0.]

…,

[ 0. 0. 0. …, 0. 0. 0.]

[ 0. 0. 0. …, 0. 0. 0.]

[ 0. 0. 0. …, 0. 0. 0.]]

(15055, 10834)

This shows even if no data value is set to -9999 in memory file, the array is initialized to 0.0 instead of no data value. Array in GeoTIFF is initialized to no data value, -9999, as expected. I think this is a bug in the memory file.

Best Answer

According to the source code, this is a problem that has been fixed in GDAL 2.0. Whether it has or not, you can get around it by pre-filling the new raster with you preferred nodata value:

outFileRead=driver.Create(outFilePath,X,Y,1,dataType,options)

tmp = gdal.AutoCreateWarpedVRT(outFileRead, src_wkt, dst_wkt)

b = outFileRead.GetRasterBand(1)
b.SetNoDataValue(-9999)
a = np.ndarray(shape=(tmp.RasterYSize, tmp.RasterXSize))
a.fill(-9999.0)
b.WriteArray(a)
gdal.ReprojectImage(inFileRead,outFileRead,inProjection,outProjection,reSamplingType)
Related Question