SPACE

GDAL

Article by:

Overview

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

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

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


Authors

Geoffrey Hunter

Dude making stuff.