2

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.

1

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
  • I am passing the DataType as "GDT_Float32" the error that I get is File "C:\Users\Gurminder\Documents\Quadium\SPI\Code\testing_stack_arr.py", line 39, in CreateGeoTiff3D DataSet = driver.Create(Name, xsize, ysize, Array.shape[0], DataType) File "C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\osgeo\gdal.py", line 388, in Create return _gdal.Driver_Create(self, *args, **kwargs) TypeError: in method 'Driver_Create', argument 6 of type 'GDALDataType' – Gurminder Bharani Oct 13 '16 at 5:10
  • @GurminderBharani From what I can tell there shouldn't be any problems. – Logan Byers Oct 13 '16 at 12:21

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.