9

If I read in a raster image as an array, then make some changes to the values in the array, how do I then save the array as a raster with the same projection information as the original array?

In particular I am conducting processing on some ISIS3 cubes from Mars. These are not projected in any of the nice SetWellKnownGeogCS options. Perhaps this makes my problem somewhat unusual, but I thought it worth documenting my solution all the same.

2 Answers 2

17

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.

1
  • 1
    Looks like PyGDAL has moved beyond using strings for things like datatype, and using None for No Data Values. Needed to tweak some things here. Aug 17, 2015 at 17:34
1

In complement to the @EddyTheB ‘s answer, mine is if the line NewArray = Array + 20 # If only it were that easy was not that easy.

Lets say, for example, if the process we want to perform do not accept the NaN values, such as an np.fft (which does NOT allow NaN values). Then, the converted array (TIFF to numpy) must to have not only no datas different to NaNs, but also a number that is not present in the image itself, simply because if you set NoData values with zeros, which would allow the FFT to work fine, but you will end up with NoData (zeros) value also inside the image, which during convertion of numpy to TIFF, might result something like:

Setting NoData as zero, when zero is also inside the images

Thus, the alternative would be a number that you could ensure that is not present in the image at all, in addition to not being NaN (e.g. 9999). So, with these considerations, the solution would be:

    NO_DATA = 9999

    def tif2array(self, input_file):
        """
        Read GeoTiff and convert to numpy.ndarray

        :param input_file: absolute path to input GeoTiff file
        :return : (np.array) image for each bands

        Source:
            - https://gist.github.com/jkatagi/a1207eee32463efd06fb57676dcf86c8
        """
        logging.info(">>>> Converting geographic format to array...")

        dataset = gdal.Open(input_file, gdal.GA_ReadOnly)
        datatype = dataset.GetRasterBand(1).DataType
        ndv = dataset.GetRasterBand(1).GetNoDataValue()

        image = np.zeros((dataset.RasterYSize, dataset.RasterXSize, dataset.RasterCount), dtype=float)

        for b in range(dataset.RasterCount):
            band = dataset.GetRasterBand(b + 1)
            image[:, :, b] = band.ReadAsArray()

        image[image == ndv] = NO_DATA

        return image, datatype, dataset

    def array2raster(self, new_rasterf_fn, dataset, array, dtype):
        """
        Save GTiff file from numpy.array

        :param new_rasterf_fn: the output image filename
        :param dataset: the original tif file, with spatial metadata
        :param array: image in numpy.array
        :param dtype:

        Source:
            - https://gist.github.com/jkatagi/a1207eee32463efd06fb57676dcf86c8
            - https://gdal.org/development/rfc/rfc58_removing_dataset_nodata_value.html
        """
        logging.info(">>>> Converting array to geographic format...")

        cols = array.shape[1]
        rows = array.shape[0]
        origin_x, pixel_width, b, origin_y, d, pixel_height = dataset.GetGeoTransform()

        driver = gdal.GetDriverByName('GTiff')

        if array.ndim == 2:
            band_num = 1
        else:
            band_num = array.shape[2]

        output_dir = os.path.dirname(new_rasterf_fn)
        if not os.path.isdir(output_dir):
            try:
                os.mkdir(output_dir)
            except OSError:
                logging.info(">>>>>> Creation of the directory {} failed".format(output_dir))
            else:
                logging.info(">>>>>> Successfully created the directory {}".format(output_dir))

        out_raster = driver.Create(new_rasterf_fn, cols, rows, band_num, dtype)
        out_raster.SetGeoTransform((origin_x, pixel_width, 0, origin_y, 0, pixel_height))
        prj = dataset.GetProjection()
        out_raster_srs = osr.SpatialReference(wkt=prj)
        out_raster.SetProjection(out_raster_srs.ExportToWkt())

        for b in range(band_num):
            outband = out_raster.GetRasterBand(b + 1)
            outband.SetNoDataValue(NO_DATA)

            if band_num == 1:
                outband.WriteArray(array)
            else:
                outband.WriteArray(array[:, :, b])
            outband.FlushCache()

        del out_raster

which will let you convert your TIFF to np.array, do your processing, then, convert it again to a geographic form, displaying no data "transparent" only outside the image, where the values have been set to 9999. Only pay attention in images with higher radiometric resolutions (16bits, 32bits, so on).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.