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