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.
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 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) # 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
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)