2

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?

  • In your dataset_out = driver.Create you specify the number of bands that you want, in this case 1 band of type Float32. Change the 1 to the number of bands you want and then use the same in dataset_out.GetRasterBand ( band ).WriteArray – Michael Stimson May 8 '17 at 3:02
  • i'll try! if i get the correct solution i will update! Thanks Michael – Italo Moletto Lobos May 8 '17 at 4:01
3

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)

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.