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.
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:
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
dataset object is set to
None in Python, you cannot guarantee that all of the data has been written to the file system.
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
TransformPoint is that you can do many points at once, making things more efficient).
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)
- Top left X
- W-E pixel resolution (pixel size y*sin(res))
- Top left Y
- N-S pixel resolution
The upper left corner of the upper left pixel will be placed at position gt, gt.
If your projection was to WGS84 (lat, lon) then:
gt: Longitude of top left corner gt: cos(alpha)*W-E resolution gt: -sin(alpha_*W-E resolution) gt: Latitude of top left corner gt: sin(alpha)*N-S resolution gt: cos(alpha)*N-S resolution