Comment créer un fichier TIFF en utilisant GDAL à partir d'un tableau numpy et en spécifiant la valeur NoData

Les deux fonctions de l'extrait de code ci-dessous, create_raster et numpy_array_to_raster devrait faire l'affaire. En termes de maintien du NoData valeur du tableau dans le raster en sortie, qui est définie sur la ou les bandes d'un raster avec le .SetNoDataValue() méthode qui, dans cet extrait de code, est utilisée dans le numpy_array_to_raster fonction. Pour plus d'informations sur l'utilisation de gdal &numpy pour le traitement raster, je recommanderais vivement le livre de Chris Garrard Geoprocessing with Python et pour une référence rapide, cette page de livre de recettes gdal/ogr est une excellente ressource.

import os
from osgeo import gdal
from osgeo import osr
import numpy

# config
GDAL_DATA_TYPE = gdal.GDT_Int32 
GEOTIFF_DRIVER_NAME = r'GTiff'
NO_DATA = 15
SPATIAL_REFERENCE_SYSTEM_WKID = 4326

def create_raster(output_path,
                  columns,
                  rows,
                  nband = 1,
                  gdal_data_type = GDAL_DATA_TYPE,
                  driver = GEOTIFF_DRIVER_NAME):
    ''' returns gdal data source raster object

    '''
    # create driver
    driver = gdal.GetDriverByName(driver)

    output_raster = driver.Create(output_path,
                                  int(columns),
                                  int(rows),
                                  nband,
                                  eType = gdal_data_type)    
    return output_raster

def numpy_array_to_raster(output_path,
                          numpy_array,
                          upper_left_tuple,
                          cell_resolution,
                          nband = 1,
                          no_data = NO_DATA,
                          gdal_data_type = GDAL_DATA_TYPE,
                          spatial_reference_system_wkid = SPATIAL_REFERENCE_SYSTEM_WKID,
                          driver = GEOTIFF_DRIVER_NAME):
    ''' returns a gdal raster data source

    keyword arguments:

    output_path -- full path to the raster to be written to disk
    numpy_array -- numpy array containing data to write to raster
    upper_left_tuple -- the upper left point of the numpy array (should be a tuple structured as (x, y))
    cell_resolution -- the cell resolution of the output raster
    nband -- the band to write to in the output raster
    no_data -- value in numpy array that should be treated as no data
    gdal_data_type -- gdal data type of raster (see gdal documentation for list of values)
    spatial_reference_system_wkid -- well known id (wkid) of the spatial reference of the data
    driver -- string value of the gdal driver to use

    '''

    print 'UL: (%s, %s)' % (upper_left_tuple[0],
                            upper_left_tuple[1])

    rows, columns = numpy_array.shape
    print 'ROWS: %s\n COLUMNS: %s\n' % (rows,
                                        columns)

    # create output raster
    output_raster = create_raster(output_path,
                                  int(columns),
                                  int(rows),
                                  nband,
                                  gdal_data_type) 

    geotransform = (upper_left_tuple[0],
                    cell_resolution,
                    upper_left_tuple[1] + cell_resolution,
                    -1 *(cell_resolution),
                    0,
                    0)

    spatial_reference = osr.SpatialReference()
    spatial_reference.ImportFromEPSG(spatial_reference_system_wkid)
    output_raster.SetProjection(spatial_reference.ExportToWkt())
    output_raster.SetGeoTransform(geotransform)
    output_band = output_raster.GetRasterBand(1)
    output_band.SetNoDataValue(no_data)
    output_band.WriteArray(numpy_array)          
    output_band.FlushCache()
    output_band.ComputeStatistics(False)

    if os.path.exists(output_path) == False:
        raise Exception('Failed to create raster: %s' % output_path)

    return  output_raster

A lire (de :Comment charger complètement un raster dans un tableau numpy ?) :

import numpy as np
from osgeo import gdal

ds = gdal.Open("mypic.tif")
cols = ds.RasterXSize
rows = ds.RasterYSize
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())

Pour écrire :

# create the output image
driver = ds.GetDriver()
outDs = driver.Create("outimage.tif", cols, rows, 1, gdal.GDT_Float32)
outBand = outDs.GetRasterBand(1)
outBand.SetNoDataValue(15)
outBand.WriteArray(myarray)
outDs.SetGeoTransform(trans)