Python GDAL: Write new raster using projection from old
This is the routine I developed to convert ISIS3 cubes to GTiffs. I expect a similar approach should work between any types of drivers (though I think the driver.Create() method might limit the choice of output file).
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)
And that's it. I can open both images in QGIS. And gdalinfo on either file shows that I have the same projection information and georeferencing.