More

Write python ndarray to raster

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.


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 toNoneor deleting it.

In your particular code the error is in the last line:del outDataonly closes your coherence numpy array. You want to closeoutDsso 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 theCreateCopy()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.

Watch the video: Batch Processing of Raster Images using Iterate Raster ModelBuilder - ArcGIS