[GIS] gdalwarp producing empty rasters when called from python script, but not via command line

gdalgdalwarppython

I'm trying to aggregate some raster data (all .tif) up to a coarser resolution (from .05 degree to .25 degree) using gdalwarp in python, but the command is not working. Instead of getting an output with a wide range of values, all values of the output are 0. The resolution and pixel depth/type are correct, but the values are not.

Here is the documentation to the gdalwarp command: http://www.gdal.org/gdalwarp.html

I have two input files which I wish to aggregate up to .25 degree resolution, producing several outputs:

  • 'NDVI_raster': The first input is a 16-bit signed raster representing NDVI, with values ranging from -10,000 to 10,000 and nodata values of -15,000.

  • 'nodata_mask': The second is a NoData mask, 32-bit float, where 1 = values of "good" data, and 0 = NoData.

With 'NDVI_raster' as input, I want to produce 7 different outputs, each representing a different statistic. I do this by calling gdalwarp 7 times, each time setting the resampling method (-r) to one of the following: average, mode, max, min, median, q1, q2. I'll call the outputs NDVI_ave.tif, NDVI_mode.tif, etc. Right now, I am using GDAL 1.10.1, which only allows average and mode, so I'm testing with just these two stats now.

Using 'nodata_mask' as input, I want to ultimately produce a QAL (quality assurance layer). To do this, I use gdalwarp, with the resampling mode set to 'average' to aggregate up to .25 degree. This results in each pixel representing the ratio of good pixels/total pixels from the input. Let's call the output QAL.

Here is what's in my code (using mode as example for the first input):

os.system('gdalwarp -tr .25 .25 -r mode -srcnodata -15000 %s %s' % (NDVI_raster, NDVI_mode))

And for the QA layer:

os.system('gdalwarp -tr .25 .25 -r average -srcnodata -15000 %s %s' % (nodata_mask, QAL))

The results are rasters with the correct resolution, projection, and pixel depth, but the pixel values are all 0.

Anyone familiar with python/gdal know what's going on?

When calling the gdalwarp command from the command line (linux), I get the desired result. When using os.system to call gdalwarp from python, I get empty rasters. So maybe something is wrong with my gdal/python bindings?

Instead of calling the command via os.system, I used subprocess. The tool via this method also appears to run smoothly, but the result is the same: a raster full of 0's.

I tried putting the gdalwarp call in a bash shell script and calling that shell script from python, but the result is a bunch of -1s instead of 0s. Weirdly enough, I had tested it before and am pretty sure it worked, but the test got deleted off my server and now I can't recreate it for some reason.

Putting the gdalwarp command in a bash shell script and then calling that shell script from the command line gives me the desired result. Calling the same shell script from python, though, doesn't. It's looking like something is wrong with python, but what and how do I fix it?

Best Answer

I had this problem too. I eventually figured out in my case it was because I had just created the on-disk image using Python gdal bindings, but I hadn't closed the gdal.Dataset in-memory object, so the write-to-disk had only partially completed. Bizarrely enough, the only way I could find to close a gdal.Dataset in Python is: del variable_name_of_dataset - so ugly!

In C there is a GDALClose() method that at the moment is not implemented by the Python GDAL API, but Rasterio does: https://github.com/sgillies/rasterio/blob/876b9a1e2bf04e349b485e05ebc4a8674ace3cf0/rasterio/_io.pyx#L1463

Also see: Why close a dataset in GDAL Python?

Related Question