The Affine Transformation

This week’s lab requires you to think about overlap areas, how to identify and define them, and how to make use of the affine transformation to avoid loading unnecessary data into an array. In addition, you will have to think about reprojecting images between different coordinate systems. In the class, we introduced how to use gdal to access raster data. We also look at several functions that allow you to access raster properties, such as the dimensions of the raster, the projection or the affine transformation. The latter can be identified using the function .GetGeoTransform() (Have a look at the location where we get the .GetGeoTransform() and inspect again the outcomes).

Now, in this chapter here we apply two new functions: (1) gdal.InvGeoTransform() and (2) gdal.ApplyGeoTransform(). This small tutorial explains what these functions do, and how you can make use of them in the context of the lab. Let’s start by opening the raster file in the same way as in the other materials.

import os
from osgeo import gdal
gdal.DontUseExceptions()

path = "D:/OneDrive - Conservation Biogeography Lab/_TEACHING/__Classes-Modules_HUB/M8_Geoprocessing-with-python/02_RasterData/WS_23-24/Assignment02/"
ds = gdal.Open(path + "tileID_410_y2000.tif")
ds

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x0000023FD14D3F60> >

cols = ds.RasterXSize
rows = ds.RasterYSize
nbands = ds.RasterCount
cols, rows, nbands

(1529, 1380, 1)

pr = ds.GetProjection()
pr

'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'

gt = ds.GetGeoTransform()
gt

(-63.86859973338657,
 0.00026949458523585647,
 0.0,
 -24.446662310500248,
 0.0,
 -0.00026949458523585647)

InvGeoTransform() and ApplyGeoTransform()

These two functions are key for the lab. When I discovered this functionality, I immediately became a fan. The biggest advantage is that their application represents very simple arithmetic calculations, which is always good when working with many datasets, as it reduces processing time immensely. The only, yet important, assumption/pre-condition is that multiple raster files and/or individual locations are in the same coordinate system. Some of the code below is inspired by the book from Garrard. We have provided the raster chapters 9 and 10 in moodle in case you want/need to dig deeper.

Let’s have a look first at gdal.InvGeoTransform(). The only argument that this function needs is the affine transformation of a raster image, in our case the variable gt:

invGT = gdal.InvGeoTransform(gt)
invGT

(236994.0, 3710.6496931091187, 0.0, -90713.0, 0.0, -3710.6496931091187)

As one would expect, the function returns a tupel of six elements. However, the values actually represent here the coefficients needed to transform geographic coordinates back into pixel/array coordinates. In other words: if you have a geo-coordinate (e.g., in lat / lon) and you want to convert it into array coordinates (which represent “row / col*), then you can use this rule and apply it to this set of coordinates. You may guess: the function we need to do so is gdal.ApplyGeoTransform(). The example below shows this for a random location represented by longitude and latitude, but you can of course imagine that this works as well with the image coordinates stored in the affine transformation gt:

# Define some example lat/lon
lat = -25.0
lon = -63.0
# Convert these geo-coordinates into array coordinates
x, y = gdal.ApplyGeoTransform(invGT, lon, lat)
x, y

(3223.0693341255246, 2053.242327727974)

Now, these are the array coordinates of the latitude/longitude values I defined. And since it is a simple multiplication of two floating point numbers, it went incredibly fast. We can interpret these coordinates also as offsets from the upper left coordinate of the original raster in x and y direction. Remember the term offset when considering the function to load a raster band into a numpy array. The last thing that is left to do is convert these floating point numbers into integers, as columns and rows in arrays are defined as such. We use the map()-function to do so. map() applies a function (first argument) to a list or a tupel of values (second argument). In our case, we want to convert our values (x, y) into int:

off_x, off_y = map(int, (x,y))
off_x, off_y

(3223, 2053)

This is all it takes. We have calculated for a random coordinate inside our image boundaries the cooresponding pixel/coordinate, expressed in row/column. Feel free to read a bit more in the Garrard book if you want to dig deeper. Yet, for the current lab you are - from a technical perspective - ready to go. Good luck and have fun!