[GIS] Save multiband raster with gdal

exportgdalgeotiff-tiffmulti-bandpython

I would like to know how to export a GeoTiff image with n bands in Python.

I have a function but it only works with one band raster.

def save_raster ( output_name, dataset, raster_data, driver ,NaN_Value,arr_projection=None):
    if arr_projection is None:
        arr_projection=[]
    else:
        if str(type(arr_projection)) == "<class 'osgeo.gdal.Dataset'>":
            srs=arr_projection.GetProjectionRef ()
        else:
            srs = arr_projection
    """
    A function to save a 1-band raster using GDAL to the file indicated
    by ``output_name``. It requires a GDAL-accesible dataset to collect 
    the projection and geotransform.
    """
    # Open the reference dataset
    g = ( dataset )
    # Get the Geotransform vector

    if raster_data is False:
        raster_data=g.ReadAsArray()
    if type(raster_data) == tuple:
        raster_data = np.array(raster_data[0])
    if str(type(g)) == "<class 'osgeo.gdal.Dataset'>":
        geo_transform = g.GetGeoTransform ()
        x_size = g.RasterXSize # Raster xsize
        y_size = g.RasterYSize # Raster ysize
        srs = g.GetProjectionRef () # Projection
    elif str(type(g)) == "<class 'affine.Affine'>":
        geo_transform = (g[2],g[0],g[1],g[5],g[3],g[4])
        RastArr = raster_data
        x_size = int(RastArr.shape[1])
        y_size = int(RastArr.shape[0])
        #PROCESS RASTERIO NUMPY
    else:
        geo_transform = (g[1][2],g[1][0],g[1][1],g[1][5],g[1][3],g[1][4])
        RastArr = np.array(g[0])
        x_size = int(RastArr.shape[2])
        y_size = int(RastArr.shape[1])
    if raster_data.ndim > 2:
        raster_data=raster_data[0]
    NaN_rast = NaN_Value
    # raster_data[raster_data == NaN_rast] = 'NaN'
    raster_data[raster_data == NaN_rast] = np.NaN
    # Need a driver object. By default, we use GeoTIFF
    driver = gdal.GetDriverByName ( driver )
    dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
            gdal.GDT_Float32 )
    dataset_out.SetGeoTransform ( geo_transform )
    dataset_out.SetProjection ( srs )
    dataset_out.GetRasterBand ( 1 ).WriteArray ( \
            raster_data.astype(np.float32) )

Any advice?

Best Answer

i've found a way using a gdal object supported by numpy arrays..

import gdal
import numpy as np
from affine import Affine
def save_multiband (output_name, dataset, raster_data, driver,
                    NaN_Value, nband, arr_projection=None):
    if arr_projection is None:
        arr_projection = []
    else:
        if str(type(arr_projection)) == "<class 'osgeo.gdal.Dataset'>":
            srs = arr_projection.GetProjectionRef()
        elif str(type(arr_projection)) == "<type 'int'>":
            srs = osr.SpatialReference()
            srs.ImportFromEPSG(arr_projection)
            srs = srs.ExportToPrettyWkt()
        else:
            srs = arr_projection
    NaN_rast = NaN_Value
    # Open the reference dataset
    g = dataset
    # Get the Geotransform vector
    if raster_data is False:
        raster_data = g.ReadAsArray()
        rast_arr = np.copy(raster_data)
    #auxiliar
    band = 0
    if type(raster_data) == tuple:
        rast_arr = np.array(raster_data[band])
    if str(type(g)) == "<class 'osgeo.gdal.Dataset'>":
        geo_transform = g.GetGeoTransform()
        x_size = g.RasterXSize  # Raster xsize
        y_size = g.RasterYSize  # Raster ysize
        srs = g.GetProjectionRef()  # Projection
    elif str(type(g)) == "<class 'affine.Affine'>":
        geo_transform = (g[2], g[0], g[1], g[5], g[3], g[4])
        rast_arr = raster_data[band,:,:]
        x_size = int(rast_arr.shape[1])
        y_size = int(rast_arr.shape[0])
    driver = gdal.GetDriverByName(driver)
    dataset_out = driver.Create(output_name, x_size, y_size, nband,
                                gdal.GDT_Float32)
    #end auxiliar
    for band in range(nband):
        if type(raster_data) == tuple:
            rast_arr = np.array(raster_data[band])
        if str(type(g)) == "<class 'osgeo.gdal.Dataset'>":
            geo_transform = g.GetGeoTransform()
            x_size = g.RasterXSize  # Raster xsize
            y_size = g.RasterYSize  # Raster ysize
            srs = g.GetProjectionRef()  # Projection
        elif str(type(g)) == "<class 'affine.Affine'>":
            geo_transform = (g[2], g[0], g[1], g[5], g[3], g[4])
            rast_arr = raster_data[band,:,:]
            x_size = int(rast_arr.shape[1])
            y_size = int(rast_arr.shape[0])
        #PROCESS RASTERIO NUMPY
        else:
            geo_transform = (g[1][2], g[1][0], g[1][1], g[1][5], g[1][3], g[1][4])
            rast_arr = np.array(g[0])
            x_size = int(rast_arr.shape[2])
            y_size = int(rast_arr.shape[1])
        rast_arr[rast_arr == NaN_rast] = np.NaN
        dataset_out.SetGeoTransform(geo_transform)
        dataset_out.SetProjection(srs)
        dataset_out.GetRasterBand(band + 1).WriteArray(
            rast_arr.astype(np.float32))

using a example data:

#assuming the data had the same extent and spatial resolution
B1 = gdal.Open(r'C\B1.tif')
B2 = gdal.Open(r'C\B2.tif')
B3 = gdal.Open(r'C\B3.tif')
B4 = gdal.Open(r'C\B4.tif')
srs = B1.GetProjectionRef()
b1, b2, b3, b4 = (B1.ReadAsArray(), B2.ReadAsArray(),
                  B3.ReadAsArray(), B4.ReadAsArray())

def GetTransform(Transform):
    transform = Transform.GetGeoTransform()
    Aff = Affine(transform[1], transform[2], transform[0], transform[4],
                 transform[5], transform[3])
    return(Aff)

B = np.stack([b1, b2, b3, b4])
save_multiband(path+'\\SR.tif', GetTransform(B1), B, "GTiff", -9999, B.ndim + 1, srs)
Related Question