[GIS] In gdal how to save 3 dimensional array as stacked raster image

gdalpython

I have constructed a 3-dimensional numpy array by stacking multiple 2-dimentional rasters.

>>> print arr_list_stack.shape
(30L, 370L, 400L)

As the arr_list_stack.shape shows that this 3-dimensional numpy array has 30 bands, by using gdal I want to save this array as raster image which has all the 30 bands.

I have a method which can save 2-dimentional arrays 

def CreateGeoTiff(Name, Array, driver, NDV, 
                  xsize, ysize, GeoT, Projection, DataType):
    Array[np.isnan(Array)] = NDV
    DataSet = driver.Create(NewFileName, xsize, ysize, 1, DataType)
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    DataSet.GetRasterBand(1).WriteArray( Array )
    DataSet.GetRasterBand(1).SetNoDataValue(NDV)
    return NewFileName

but CreateGeoTiff() fails to save 3-dimentional arrays and returns the error 

  File "C:\Users\Documents\Quadium\SPI\Code\testing_stack_arr.py", line 39, in CreateGeoTiff
    DataSet.GetRasterBand(1).WriteArray( Array )
  File "C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\osgeo\gdal.py", line 1235, in WriteArray
    callback_data = callback_data )
  File "C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\osgeo\gdal_array.py", line 359, in BandWriteArray
    raise ValueError("expected array of dim 2")

What changes should I make in CreateGeoTiff() method that will enable me to save 3-dimentional array.

Best Answer

Look at the documentation for driver.Create. Your current code is trying to write a 3D array to a single band. That is what raise ValueError("expected array of dim 2") is telling you. You need to create a tiff with bands equal to the number of "stacked images". Then write a single image to each band. You cannot write all the bands at once.

Here is the same function you wrote, updated to reflect this change. There are also modifications because NewFileName isn't the parameter name....

def CreateGeoTiff(Name, Array, driver, NDV, 
                  xsize, ysize, GeoT, Projection, DataType):
    Array[np.isnan(Array)] = NDV
    DataSet = driver.Create(Name, xsize, ysize, Array.shape[0], DataType)
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    for i, image in enumerate(Array, 1):
        DataSet.GetRasterBand(i).WriteArray( image )
        DataSet.GetRasterBand(i).SetNoDataValue(NDV)
    DataSet.FlushCache()
    return Name

You can make this less verbose because your array will tell you the xsize and ysize.

def CreateGeoTiff(Name, Array, driver, NDV, 
                  GeoT, Projection, DataType):
    Array[np.isnan(Array)] = NDV
    DataSet = driver.Create(Name, Array.shape[2], Array.shape[1], Array.shape[0], DataType)
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    for i, image in enumerate(Array, 1):
        DataSet.GetRasterBand(i).WriteArray( image )
        DataSet.GetRasterBand(i).SetNoDataValue(NDV)
    DataSet.FlushCache()
    return Name
Related Question