Sunday 27 December 2015

write python ndarray to raster


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.



Answer




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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...