2

I'm trying to use GDAL/python API to open a raster, do some math, then write the two resulting arrays as rastes. One is to be an 8-bit raster and the other 16-bit signed. However, only the 8-bit raster is coming out correctly. The 16-bit raster comes back as a raster of 0's. Here is the relevant code:

import numpy as np
import gdal
gdal.UseExceptions() # enable exceptions to report errors
drvtif = gdal.GetDriverByName("GTiff")

# open raster, read as array, do some math, return two arrays: arr1 and arr2

# arr1 and arr2 are both float32, though not necessary

print np.unique(arr1)
# prints [-15000, -10000, ..., 10000]
print np.unique(arr2)
# prints [1, 2, ..., 11]

# now write arr1 and arr2 as .tif:

outNDVI = '/path/to/outputs/output1.tif':
ndvidrv = drvtif.Create(outNDVI, ncols, nrows, 1, 3) # 3 = 16 bit signed
ndvidrv.SetGeoTransform(hc_gt)
ndvidrv.SetProjection(proj2)
ndvidrv.GetRasterBand(1).SetNoDataValue(-15000)
ndvidrv.GetRasterBand(1).WriteArray(arr1)
del ndvidrv

outQA = '/path/to/outputs/output2.tif'
qadrv = drvtif.Create(outQA, ncols, nrows, 1, 1) # 1= 8 bit
qadrv.SetGeoTransform(hc_gt)
qadrv.SetProjection(proj2)
## qadrv.GetRasterBand(1).SetNoDataValue(-15000)
qadrv.GetRasterBand(1).WriteArray(arr2)
del qadrv

outQA turns out as expected, with the correct range of values and correct projection, but outNDVI is just a raster of 0's (in the correct projection). The values should range from -10000 to 10000 (with NoData= -15000).

Any thoughts?

Deleting ndvidrv seemed to do the trick and I've updated above to reflect that.

4
  • 2
    Are you running this code from an IDE or debugger? They often keep the python process alive after the script finishes so the datasets do not get closed properly. You may be seeing this if the dataset is not getting closed and written to disk - trac.osgeo.org/gdal/wiki/…
    – user2856
    Mar 11, 2016 at 4:45
  • How are you viewing the results? A NoData value -15000 is bogus for a Byte raster (qadrv)
    – Mike T
    Mar 11, 2016 at 5:38
  • @MikeT hi Mike, I updated the post with some information about the arrays. You're right, I didn't mean to set the NoData value to -15000 in outQA (fixed), but that raster turns out to be correct anyway. It's the first case (outNDVI) that results in a raster of 0s.
    – user20408
    Mar 11, 2016 at 6:00
  • @Luke hi Luke, I'm calling my Python script from the command line of a Linux server. I don't normally have this problem though I could be doing something different that i just can't see.
    – user20408
    Mar 11, 2016 at 6:08

1 Answer 1

2

I think you might need to flush the band out and then close ("delete") the datasource try adding:

ndvidrv.GetRasterBand(1).FlushCache() ## flush to finish writing

del ndvidrv
3
  • 2
    @user20408 del ndvidrv is sufficient. The wiki notes that not all drivers support flushing and that it doesn't guarantee that the data are written to disk, so the preferred method is to deallocate (i.e. del ndvidrv or ndvidrv = None)
    – user2856
    Mar 11, 2016 at 8:59
  • 1
    Yes, it is enough del ndvidrv or ndvidrv = None after ndvidrv.GetRasterBand(1).WriteArray(arr1).
    – xunilk
    Mar 11, 2016 at 10:19
  • Thanks Josh, Luke, @xunilk... del ndvidrv worked and I've updated the question. Thank you all
    – user20408
    Mar 11, 2016 at 16:06

Your Answer

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

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