Python GDAL: écrire un nouveau raster en utilisant la projection de l'ancien

Si je lis dans une image raster en tant que tableau, puis modifie les valeurs du tableau, comment puis-je enregistrer le tableau en tant que raster avec les mêmes informations de projection que le tableau d'origine?

En particulier, je procède au traitement de certains cubes ISIS3 de Mars. Ceux-ci ne sont projetés dans aucune des belles options SetWellKnownGeogCS. Cela rend peut-être mon problème quelque peu inhabituel, mais j'ai pensé qu'il valait quand même la peine de documenter ma solution.

C'est la routine que j'ai développée pour convertir des cubes ISIS3 en GTiffs. J'espère qu'une approche similaire devrait fonctionner entre tous les types de pilotes (bien que je pense que la méthode driver.Create () pourrait limiter le choix du fichier de sortie).

import numpy as np
import gdal
from gdalconst import *
from osgeo import osr

# Function to read the original file's projection:
def GetGeoInfo(FileName):
    SourceDS = gdal.Open(FileName, GA_ReadOnly)
    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)
    return NDV, xsize, ysize, GeoT, Projection, DataType

# Function to write a new file.
def CreateGeoTiff(Name, Array, driver, NDV, 
                  xsize, ysize, GeoT, Projection, DataType):
    if DataType == 'Float32':
        DataType = gdal.GDT_Float32
    NewFileName = Name+'.tif'
    # Set nans to the original No Data Value
    Array[np.isnan(Array)] = NDV
    # Set up the dataset
    DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType )
            # the '1' is for band 1.
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    # Write the array
    DataSet.GetRasterBand(1).WriteArray( Array )
    DataSet.GetRasterBand(1).SetNoDataValue(NDV)
    return NewFileName

# Open the original file
FileName = 'I29955002trim.cub'    # This is the ISIS3 cube file
                                  # It's an infra-red photograph
                                  # taken by the 2001 Mars Odyssey orbiter.
DataSet = gdal.Open(FileName, GA_ReadOnly)
# Get the first (and only) band.
Band = DataSet.GetRasterBand(1)
# Open as an array.
Array = Band.ReadAsArray()
# Get the No Data Value
NDV = Band.GetNoDataValue()
# Convert No Data Points to nans
Array[Array == NDV] = np.nan

# Now I do some processing on Array, it's pretty complex 
# but for this example I'll just add 20 to each pixel.
NewArray = Array + 20  # If only it were that easy

# Now I'm ready to save the new file, in the meantime I have 
# closed the original, so I reopen it to get the projection
# information...
NDV, xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)

# Set up the GTiff driver
driver = gdal.GetDriverByName('GTiff')

# Now turn the array into a GTiff.
NewFileName = CreateGeoTiff('I29955002trim', NewArray, driver, NDV, 
                            xsize, ysize, GeoT, Projection, DataType)

Et c'est tout. Je peux ouvrir les deux images dans QGIS. Et gdalinfo sur l'un ou l'autre fichier montre que j'ai les mêmes informations de projection et de géoréférencement.