GeoTIFF and Python GDAL

If you program in Python, then the Python binding to the Geospatial Data Abstraction library (GDAL) and Numpy are indispensable tools for reading, writing and manipulating raster image in GeoTIFF format.

This summarizes my brief experience working with them.

A GeoTIFF image file may contain one or more bands. Each band can be viewed as a separate 2D data array but every band in the same file will have the same array dimension (or shape, in Numpy lingo).

All image pixels in a band (data element of array) are of the same data type, e.g. byte, int16, etc., but one band may be of different data type from another in the same file.

By convention, if a GeoTIFF file contains a RGB raster image, it will have 3 bands (4 if it is RGBA, with the addition of alpha channel for pixel transparency) of byte data — band 1 for red, 2 for green and 3 for blue.

GDAL API

To read GeoTIFF metadata (georeference information)

from osgeo import gdal
ds = gdal.Open("sample.tif", gdal.GA_ReadOnly)
(X, deltaX, rotation, Y, rotation, deltaY) = ds.GetGeoTransform()
srs_wkt = ds.GetProjection()
Nx = ds.RasterXSize
Ny = ds.RasterYSize
ds = None

GeoTransform() returns a tuple where (X, Y) are corner coordinates of the image indicating the origin of the array, i.e. data element[0,0]. But, which corner?

If deltaX is positive, X is West. Otherwise, X is East. If deltaY is positive, Y is South. Otherwise, Y is North. In other words, when both deltaX and deltaY is positive, (X, Y) is the lower-left corner of the image.

It is also common to have positive deltaX but negative deltaY which indicates that (X, Y) is the top-left corner of the image.

There are several standard notations for SRS, e.g. WKT, Proj4, etc. GetProjection() returns a WKT string. The meaning and unit-of-measurement for  X, Y, deltaX and deltaY are based on the spatial reference system (SRS) of the image. For example, if the SRS is Stereographic projection, they are in kilometres.

To convert from WKT notation to Proj4

from osgeo import osr
srs = osr.SpatialReference()
srs.ImportFromWkt(srs_wkt)
srs_proj4 = srs.ExportToProj4()

To load a multi-band GeoTIFF image file as Numpy arrays

from osgeo import gdal
arys = []
ds = gdal.Open('sample.tif', gdal.GA_ReadOnly)
for i in xrange(1, ds.RasterCount+1):
arys.append(ds.GetRasterBand(i).ReadAsArray())
# for a single-band GeoTIFF
ary = ds.GetRasterBand(1).ReadAsArray()

Unlike mathematical convention for coordinates (X, Y), Numpy array indexing uses [y-offset][x-offset].

To save a 2-D Numpy array as 1-band GeoTIFF file

Included correction from Pete

from osgeo import gdal
 # data exists in 'ary' with values range 0 - 255
 # Uncomment the next line if ary[0][0] is upper-left corner
 #ary = numpy.flipup(ary)
 Ny, Nx = ary.shape
 driver = gdal.GetDriverByName("GTiff")
 ds = driver.Create('output.tif', Nx, Ny, 1, gdal.GDT_Byte)
 ds.SetGeoTransform( ... ) # define GeoTransform tuple
 ds.GetRasterBand(1).WriteArray(ary)
 ds = None

To slice from a bigger image (array)

subary = ary[4:6, 5:8]

To change coordinates from one SRS to another

Use pyproj
For example, to convert 2 lists of longitudes and latitudes to X, Y coordinates in polar stereographic projection:

import pyproj
# Specify target SRS
pp = pyproj.Proj('+proj=stere +lat_ts=-60.0 +lat_0=-90.0 +lon_0=28.0 +k_0=1.0')

# Origin of 2D-data array: longitudes[0,0] latitudes[0,0]
X, Y = pp(longitudes, latitudes)

# To reverse the operation
longitudes, latitudes = pp(X, Y, inverse=True)

About these ads