8

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?

Improve this question
19
  • 2
    It can't be a binding problem, you didn't use them. You execute the gdalwarp command. You should check the error status of the gdalwarp command, see: stackoverflow.com/questions/3791465/…
    – Zoltan
    Jan 19, 2016 at 18:09
  • 2
    To expand on what Zoltan is saying, you are essentially calling the same program (gdalwarp). os.system actually starts up a shell and executes the command given. In addition to the return value, you should probably verify the variables that you are trying to pass in (ie. any quoting issues? etc.)
    – Evil Genius
    Jan 19, 2016 at 18:23
  • You should also consider moving from os.system to subprocess which is the newer module and contains a number of (security) improvements.
    – Kersten
    Jan 20, 2016 at 9:38
  • @Zoltan what error message? I didn't get an error message because the command technically worked, just the output was not right (values of 0).plz see my update. thanks!
    – user20408
    Jan 20, 2016 at 18:40
  • 2
    One issue is that I think you're going to need to specify -ot Float32 or something similar for the quality mask, because you can't represent fractions with an integer raster.
    – Nat Wilson
    Jan 27, 2016 at 3:06

2 Answers 2

Reset to default
3

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?

Improve this answer
2

You don't say how you start your python/gdalwarp script. I discovered that a cronjob will not always have the same environment as my command line environment. I had to start creating runtime environments for these kinds-of scripts. With that being said, if you start you script from, say, an icon on your desktop, then it may not have the same runtime environment as your command line environment. "PYTHONPATH" may be one of the environment variables that you have to set. In addition, you may have to set variables for gdalwarp. Finally, your data files may not be in the correct location. You may have to set an absolute path such as /xxx/xxxx/NDVI_raster or use tilda ~/NDVI_raster. Like PYTHONPATH, you may also have to setup up PATH and other environment variables. Make sure to "export or source" the sames settings at the beginning of your script.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.