GDAL

Overview

The Geospatial Data Abstraction Library (GDAL) is a open-source software library for handling geospatial imagery data.

RPC Metadata

The command-line tool gdalinfo will print the RPC metadata if present.

Basic Reading/Writing Of GeoTiff Files

Updating the same file is the simplest example, as you don’t have to worry about creating a driver:

1
2
3
4
5
6

import gdal
ds = gdal.Open('image.tif', gdal.GA_Update)
img = ds.GetRasterBand(1).ReadAsArray()
# Change img as you see fit (it's a Numpy array)
ds.GetRasterBand(1).WriteArray(img)
ds = None

Always Make Sure To dataset=None

Until a dataset object is set to None in Python, you cannot guarantee that all of the data has been written to the file system.

1
2
3
4
5
6
7

ds = gdal.Open('image.tif')

## Do stuff here...

## Make sure to set the dataset to None, this will make sure all data
# is written out of buffers into file
ds = None

Converting x, y Pixel Coordinates To Lat/Lons Using RPCs

RPC metadata can be embedded within an non-orthorectified image file. GDAL can use this RPC metadata to convert (or transform) co-ordinates in x, y pixels space into latitude/longitudes (and also to perform the reverse, to calculate x/y pixel co-ordinates from lat/lons).

You can do this in GDAL with a Transformer and the TransformPoints() function. However, these are quite poorly documented online (the benefit of TransformPoints() over TransformPoint is that you can do many points at once, making things more efficient).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

ds = gdal.Open('image_with_rpc_metadata.tif')

tr = gdal.Transformer(ds, None, ['METHOD=RPC'])

## TransformPoints requires the points to be in the same array, in the form:
# [ [x0, y0], [x1, y1], ... ]
xy_points = np.array([[0, 0], [5, 10], [10, 20]])
transformed_points, shape = tr.TransformPoints(0, xy_points)

# tranformed_points will be a list of tuples, in the form:
# [ (lon0, lat0, height0), (lon1, lat1, height1), ... ]
# where (lon0, lat0, height0) maps to (x0, y0)

The geotransform

1. Top left X
2. W-E pixel resolution (pixel size y*sin(res))
3. Rotation
4. Top left Y
5. Rotation
6. N-S pixel resolution

The upper left corner of the upper left pixel will be placed at position gt[0], gt[3].

If your projection was to WGS84 (lat, lon) then:

gt[0]: Longitude of top left corner
gt[1]: cos(alpha)*W-E resolution
gt[2]: -sin(alpha_*W-E resolution)
gt[3]: Latitude of top left corner
gt[4]: sin(alpha)*N-S resolution
gt[5]: cos(alpha)*N-S resolution