Working with rasters using GDAL and Python

01.04.2010 09:10 ·  GIS  ·  gdal, python, howto

GDAL is a free library for working with raster and vector data. The command line utilities included in the library are widely used to perform a variety of tasks. Thanks to the developed API, you can work with GDAL functions from many programming languages. This article is dedicated to using the GDAL API in Python and is based on the GDAL Raster API Tutorial.

Preparation

I assume that all software is installed using the OSGeo4W installer. To use GDAL from Python, the appropriate language bindings are required. In the OSGeo4W the required package is called gdal-python (for GDAL 1.5.x) or gdal16-python (for GDAL 1.6.x). I recommend using the 1.6.x version. We also need the NumPy package to perform various operations on the raster, which is also available via OSGeo4W.

The GDAL Python API consists of five core modules and five optional modules (exist for compatibility with older versions):

They can be imported using the following code:

# main modules
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo import gdalconst

# compatibility modules. They will be removed in GDAL 2.0
import gdal
import ogr
import osr
import gdalnumeric
import gdalconst

When using GDAL version 1.5 and above, it is recommended to import “main” modules. If the code needs to be compatible with older GDAL versions, modules can be imported using the following approach:

try:
    from osgeo import gdal
except ImportError:
    import gdal

Opening a file

The gdal.Open() function is used to open a raster. It takes the full path to the raster as the first argument and an optional constant describing the opening mode as the second argument. If the constant is omitted, a read-only mode is assumed. On success, the function returns a GDALDataset object, otherwise it returns None.

import sys
from osgeo import gdal
gdalData = gdal.Open("/home/alex/test/input.tif", GA_ReadOnly)
# or without specifying open mode
# gdalData = gdal.Open("/home/alex/test/input.tif")

# checking that everything is fine
if gdalData is None:
    print "ERROR: can't open raster"
    sys.exit(1)

Obtaining raster information

The GDALDataset object allows you to get information about the number of raster bands, extract raster data and get metadata. All the functions used have self-explanatory names, so I hope there is no need to describe them in detail. Instead, I will give a small example:

print "Driver short name", gdalData.GetDriver().ShortName
print "Driver long name", gdalData.GetDriver().LongName
print "Raster size", gdalData.RasterXSize, "x", gdalData.RasterYSize
print "Number of bands", gdalData.RasterCount
print "Projection", gdalData.GetProjection()
print "Geo transform", gdalData.GetGeoTransform()

The last function used in the example above, GetGeoTransform(), returns a list of 6 numbers describing the georeferencing transformation of the raster:

Extracting raster/band data

Any raster can be represented as an array. A single-band raster is represented as a two-dimensional array of size X×Y; a multi-band raster is represented as a multidimensional array. GDAL allows you to access either the whole grid or each band separately. The result of such an operation is a NumPy array of the corresponding dimension:

# get the whole raster
raster = gdalData.ReadAsArray()
# get specific band
gdalBand = gdalData.GetRasterBand(1)
band_1 = gdalBand.ReadAsArray()

You can also process the image on a line-by-line basis or on a block-by-block basis:

line = gdalBand.ReadAsArray(xoffeset, yoffset, xsize, ysize)

where:

That is, if both offsets are equal to zero, and xsize and ysize are equal to the height and width of the raster, we will get the first row of the array in the line variable. And by changing the value of yoffset we can iterate through all the rows in the array.

When extracting data, remember that the data type in the array is the same as the raster data type, and in some cases this can lead to errors. For example, if you subtract two unsigned integer arrays, you may well get an incorrect result due to an integer overflow. Therefore, you must be careful when extracting data, and convert to a different data type if necessary. For example, the following line will extract the raster band and force the datatype to Float32 (32-bit float).

band_1 = gdalBand.ReadAsArray().astype(numpy.float32)

Conversion to any other data type can be done in the same way. A full list of supported data types can be found in the NumPy documentation.

Raster processing

To work with arrays (and as we already know, rasters are represented as arrays), we will use NumPy — a Python module that provides large multidimensional arrays and matrices, as well as an extensive library of mathematical functions for manipulating these arrays.

Almost all functions in NumPy operate element-wise, which makes array processing much easier. For example, to get the sum of two arrays with the same dimensions, it is enough to write a single line of code:

result = numpy.add(array1, array2)

Subtraction, multiplication, division, and other operations are performed in the same way. Unary operations are also performed element-wise, e.g. after executing the following line:

r = sin(array1)

the array r will contain the sine values of each element of the array array1.

Here is a more complex example — processing a raster on a pixel-by-pixel basis:

gdalData = gdal.Open("/home/alex/test/input.tif")
# raster size
xsize = gdalData.RasterXSize
ysize = gdalData.RasterYSize
# get raster data as array
raster = gdalData.ReadAsArray()
# iterate over pixels
for col in range(xsize):
    for row in range(ysize):
        # if the pixel value is 5, change it to 10, otherwise do nothing
        if raster[row, col] == 5:
            raster[row, col] = 10

The same result can be achieved with NumPy, and it will be much faster:

mask = numpy.equal(raster, 5)
numpy.putmask(raster, mask, 10)

Writing files

There are two ways to create a new raster: using the CreateCopy() or Create() functions. When using CreateCopy() you should specify the “reference” raster whose metadata will be used to create a new file. When using the Create() function, it is necessary to provide all the metadata manually and only then create a new file.

Note that not all drivers support the Create() method, so you should first check if the driver supports the required method:

format = "GTiff"
driver = gdal.GetDriverByName(format)
metadata = driver.GetMetadata()
if metadata.has_key(gdal.DCAP_CREATE) and metadata[gdal.DCAP_CREATE] == "YES":
    pass
else:
    print "Driver %s does not support Create() method." % format
    sys.exit(1)
# similarly check for CreateCopy support
if metadata.has_key(gdal.DCAP_CREATECOPY) and metadata[gdal.DCAP_CREATECOPY] == "YES":
    pass
else:
    print "Driver %s does not support CreateCopy() method." % format
    sys.exit(1)

Some drivers may be read-only and may not support either of these methods.

When using CreateCopy(), all necessary information is taken from the “reference” image, but it is possible to pass additional format specific parameters.

# "reference" raster
inputData = gdal.Open("/home/alex/test/input.tif")
# create a new raster with the same parameters as inputData has
outputData = driver.CreateCopy("/home/alex/test/output.tif", inputData, 0)
# the same as above, but with passing additional creation options
# outputData = driver.CreateCopy("/home/alex/test/output.tif",inputData, 0, ["TILED=YES", "COMPRESS=PACKBITS"])

# close all datasets and free memory
inputData = None
outputData = None

When using Create() it is necessary to manually set the number of channels, the image size and the desired data type (byte, long integer, float…).

# raster size
cols = 512
rows = 512
# band count
bands = 1
# raster data type
dt = gdal.GDT_Byte
# create new raster
outData = driver.Create("/home/alex/test/out.tif", cols, rows, bands, dt)

When the raster is created, we can add projection information:

outData.SetProjection(proj)
outData.SetGeoTransform(transform)

and write the data to the band:

outData.GetRasterBand(1).WriteArray(raster)
# closing dataset when done
outData = None

If there are several bands, it is more convenient to write data using a loop:

bands = 3
for i in range(bands):
    outData.GetRasterBand(i + 1).WriteArray(raster[i])
outData = None

Example. Simple raster calculator

Using the information from the post, let’s try to develop a simple Python script to add two raster bands. The script takes two parameters: the source image and the path to save the result.

# -*- coding: utf-8 -*-

#!/usr/bin/env python

import sys
import numpy
import osgeo.gdal as gdal

def isMultiband(path):
    """Checks whether raster is a multiband raster"""
    gdalData = gdal.Open(path)
    if gdalData.RasterCount < 2:
        print "ERROR: raster must contain at least 2 bands"
        return False
    return True


def bandAsArray(path, bandNum):
    """Extracts the band data from the raster"""
    gdalData = gdal.Open(path)
    gdalBand = gdalData.GetRasterBand(bandNum)
    array = gdalBand.ReadAsArray().astype(numpy.float32)
    gdalBand = None
    gdalData = None
    return array


def saveRaster(outPath, referencePath, array):
    """Saves array as GeoTiff raster to the outPath.
    Projection information is taken from the raster located at referencePath
    """
    gdalData = gdal.Open(referencePath)
    projection = gdalData.GetProjection()
    transform = gdalData.GetGeoTransform()
    xsize = gdalData.RasterXSize
    ysize = gdalData.RasterYSize
    gdalData = None

    format = "GTiff"
    driver = gdal.GetDriverByName(format)
    metadata = driver.GetMetadata()
    if metadata.has_key(gdal.DCAP_CREATE) and metadata[gdal.DCAP_CREATE] == "YES":
        outRaster = driver.Create(outPath, xsize, ysize, 1, gdal.GDT_Float32)
        outRaster.SetProjection(projection)
        outRaster.SetGeoTransform(transform)
        outRaster.GetRasterBand(1).WriteArray(array)
        outRaster = None
    else:
        print "Driver %s does not support Create() method." % format
        return False


if __name__ == "__main__":
    args = sys.argv[1:]
    inPath = args[0]
    outPath = args[1]
    if isMultiband(inPath):
        band1 = bandAsArray(inPath, 1)
        band2 = bandAsArray(inPath, 2)
        res = numpy.add(band1, band2)

        if not saveRaster(outPath, inPath, res):
            print "ERROR: saving failed"
            sys.exit(1)

        print "Success"
    else:
        print "Input raster has only 1 band"
        sys.exit(1)

For those who want to learn more about raster processing with GDAL and Python, I suggest the following links as a starting point:

⮜ Prev
Next ⮞