3

I think I understood the steps to write an array to a raster file using python. however it doesn't work in my case and I don't know why. I have an array that I want to write to a Geotiff with Swiss Oblique Mercator Projection. I have a DEM as Geotiff in the exact same location with same projections and Geotransform, so I thought I just read all parameters from that DEM and use them to write my new raster. Here's the code that I'm, using:

driver.Register()
fn = '/media/LACIE_SHARE/davos2015.tif'
inDs = gdal.Open(fn)
if inDs is None:
    print 'Could not open ' + fn 
    sys.exit(1)

rows=inDs.RasterYSize
cols=inDs.RasterXSize

driver=inDs.GetDriver()
outDs=driver.Create('/home/grassdata/test.tif',cols,rows,1,gdal.GDT_Float32)

if outDs is None:
    print "crap"
    sys.exit(1)
outBand=outDs.GetRasterBand(1)
outData=coh0 #coh0 is the array I want to convert to raster. it's an array with float32 entries. 

outBand.WriteArray(outData) # this gives 0
outBand.FlushCache()
outBand.SetNoDataValue(-99) #this gives 0
outDs.SetGeoTransform(inDs.GetGeoTransform()) #this gives 0
outDs.SetProjection(inDs.GetProjection()) #this gives 0
del outData

the lines where I commented accordingly raise 0 as result and I don't understand why.

  • 1) Is test.tif created at all? 2) Try gdal.GetDriverByName('GTiff') instead of getting the driver from the input dataset. 3) You can try gdal.UseExceptions() at the beginning of your script to catch exceptions on the fly. 4) HAve you had a look at GDALs Python Gotchas – Kersten Sep 2 '15 at 14:46
  • 1) Yes it's created 2)tried but gives same result. The driver I read from input is the same as the one given by gdal.GetDriverByName('GTiff'). 3) No exceptions raised, but result of writeArray is still 0 :-( 4) nothing useful concerning this problem in the Gotchas.. – uetli Sep 3 '15 at 6:07
  • The only other thing I would suspect in the code is that FlushCache() is usually called on the dataset and not the band. Are you able to provide a small subset of the DEM as well as the corresponsing coherence(?) coh0 array? – Kersten Sep 3 '15 at 11:19
  • @kersten Then how would I have to do it to Flush the whole dataset?.. I'm confused :-) I provided you with the DEM and coh0 here: [link] dropbox.com/sh/iyfmo1ax6ww8j96/AADk3jQ2R2PtQ93rjcXgMjeXa?dl=0 hope that works. you can read the coherence using the following line of code: coh0=np.fromfile('/yourpath/coh0.npy',dtype='<f32') coh0=coh0[20:].reshape(3200,3700) thanks for trying to help – uetli Sep 3 '15 at 13:47
2

Your issue is rooted in a well known GDAL Python gotcha - a dataset needs to be closed for it to be written to disk.

In Python this happens when the object goes out of scope and is garbage collected or when you manually dereference it. This is usually done by setting it to None or deleting it.

In your particular code the error is in the last line: del outData only closes your coherence numpy array. You want to close outDs so it is written to disk. FlushCache() is not necesary in your case since that only influences the behaviour of the blockcache which GDAL handles internally.

Your last lines should look like this:

outBand.WriteArray(coh0)
outBand.SetNoDataValue(-99)
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())

outBand = None
outDs = None

Additional notes:

  • If you want to create a copy of an existing dataset and just fill it with new data you could also use the CreateCopy() method to save some lines of code.
  • GDAL is very mighty but at the same time very unpythonic. rasterio is an alternative way to use the GDAL bindings in a much more pythonic way.
  • You are right, such a stupid mistake! Thank you so much, it works now! you saved my weekend – uetli Sep 4 '15 at 15:16
  • I've made that mistake at least a dozen times myself and still took some time to catch it. The devil really is in the details with GDAL ;) – Kersten Sep 4 '15 at 15:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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