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")
 # Final argument is optional but will produce much smaller output file
 ds = driver.Create('output.tif', Nx, Ny, 1, gdal.GDT_Byte, ['COMPRESS=LZW'])

 # this assumes the projection is Geographic lat/lon WGS 84
 srs = osr.SpatialReference()
 srs.ImportFromEPSG(4326)
 ds.SetProjection(srs.ExportToWkt())

 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

4 Responses to “GeoTIFF and Python GDAL”

  1. Pete Says:

    I think you may have the Ny and Nx the wrong way round for the driver.Create line. Otherwise a helpful introduction.

  2. rnebot Says:

    Very useful. Just two small corrections: in the tuple returned by “ds.GetGeoTransform()”, “deltaY” should go after “Y” and then “rotation”. The second one is the case of “ImportFromWKT”, which I think should be “ImportFromWkt”.

    • cynici Says:

      Hi, thanks for the correction. I’ve fixed ImportFromWkt() in the text. I have looked up the GDAL reference for GetGeoTransform() again. The order of arguments are indeed counter-intuitive.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Follow

Get every new post delivered to your Inbox.

%d bloggers like this: