Ctrl+K

The Geotransformation

The Geotransformation#

This 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 the notes of the class, we have introduced how to use gdal to access raster data. In addition, we looked 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, 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 second lab. Let’s start by opening the raster file in the same way as in the other materials.

import os
from osgeo import gdal

path = "D:/OneDrive - Conservation Biogeography Lab/_TEACHING/__Classes-Modules_HUB/M8_Geoprocessing-with-python/Week_02/WS_23-24/Assignment02/"
os.chdir(path)

ds = gdal.Open(path + "tileID_410_y2000.tif")

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

pr = ds.GetProjection()

gt = ds.GetGeoTransform()

InvGeoTransform() and ApplyGeoTransform()#

These two functions are key for our 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. The only, yet important, assumption/pre-condition is that multiple raster files and/ 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/line 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().

# 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 we need to do is convert these floating point numbers into integers, as columns and rows in arrays are defined as such. We use the map() 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. 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!

Contents