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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
user20408
  • 1,535
  • 1
  • 22
  • 41
  • 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: http://stackoverflow.com/questions/3791465/python-os-system-for-command-line-call-linux-not-returning-what-it-should – Zoltan Jan 19 '16 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 '16 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 '16 at 09: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 '16 at 18:40
  • @EvilGenius hi, thanks. i printed the statement to the terminal and used the exact same command on the command line, where it worked – user20408 Jan 20 '16 at 18:40
  • @Kersten thank you. At your suggestion, I switched to using subprocess but am still getting the same problem – user20408 Jan 20 '16 at 18:41
  • I tried this, and it worked for me using both subprocess.call and the obsolete os.system. I'm using GDAL 1.11.2. Can you provide some test data demonstrating the problem? – Nat Wilson Jan 27 '16 at 03:05
  • 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 '16 at 03:06
  • Is file-locking an issue here? I have an unsolved Python GDAL question which could possibly share some characteristics with this one: http://gis.stackexchange.com/questions/175952/pyqgis-save-raster-as-rendered-image-then-use-gdal-tools-on-it – Tom Chadwin Jan 27 '16 at 12:11
  • @NatWilson hi Nat, thanks. I can send you the raster I'm using as input and the entire python script, so you can test it in python and on the command line, but how do I do that? – user20408 Jan 27 '16 at 17:18
  • @user20408 I think it would be best to make the data available as part of your question, rather than sending it to me personally. Ideally, a script that generates a tiny raster that exhibits the problem, but if that's not feasible, perhaps a link to it (e.g. on Dropbox) would be acceptable? – Nat Wilson Jan 27 '16 at 17:49
  • @NatWilson ok will work on that. in the meantime, i did try specifying the output type -ot Float32 and it still didn't work. Also, please see my 4th update. Thanks much. I'd love to give you these 50 reputation points as I really need to figure this out. – user20408 Jan 27 '16 at 19:15
  • I just tested this with a simple 3 band Byte color image and noticed that gdalwarp V1.11.3 on windows gives me strange results. Not only is it different between python os.system and the prompt but also the prompt gives slightly incorrect answers. I also tested with gdalwarp v2.0.1 on linux and that seems to work allright. By any chance you can upgrade to 2.0.1? The strange behaviour on 1.11 might be worth a new question, with some simple test data. – tilt Jan 31 '16 at 10:05
  • Did you try to wrap your -srcnodata -15000 value argument into quotes? Because it may look like a command option if not escaped on the command line. Like -srcnodata "-15000"? You didn't show how you used your subprocess command, but the arguments should always be strings. – Frank Feb 01 '16 at 07:12
  • @tilt Hi Tilt, I should have mentioned this but I actually have been working on GDAL 2.0.1 – user20408 Feb 01 '16 at 16:13
  • Then indeed we need some sample data and more information on the platform (from your bash usage I understand it must be some *ix flavour) and python version. – tilt Feb 01 '16 at 20:47
  • try using subprocess with shell = True https://docs.python.org/2/library/subprocess.html – RutgerH Feb 02 '16 at 21:35
  • With GDAL trunk version you have a new option to use GDAL utilities as a library https://trac.osgeo.org/gdal/wiki/rfc59.1_utilities_as_a_library. – user30184 Feb 15 '16 at 14:21
  • To complete @EvilGenius, I usually find that errors with gdalwarp come from not having the GDAL_DATA environmental variable set – Jose Apr 10 '17 at 11:13

2 Answers2

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?

Hazzles
  • 131
  • 2
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.

Community
  • 1
Greg
  • 1,035
  • 6
  • 9