Georeference and rotate an image using GDAL and Python

Question

I have a drone image with no spatial reference (ungeoreferenced). I know from some calculations the coordinates of the upper left corner, and X and Y resolutions. I was able to georeference the image using Python by setting the geotransform. See example below;

ulx = 5
uly = 10
xres = 0.1
yres = -0.1
xskew = 0
yskew = 0
ulx, xres, xskew, uly, yskew, yres = gtf
raster.SetGeoTransform(gtf)

The above code placed the image in its approximate location. Then I noticed that the image is supposed to be rotated (about 120 degrees). I tried to set the x and y skew to 120

ulx = 5
uly = 10
xres = 0.1
yres = -0.1
xskew = 120
yskew = 120
ulx, xres, xskew, uly, yskew, yres = gtf
raster.SetGeoTransform(gtf)

This takes the image away from its location and the pixels are extremely larger. What went wrong and how to fix it?

• PyMapr: Have a look at affine geotransform. gdal.org/user/raster_data_model.html#affine-geotransform

Answer 1

I just faced a similar problem and found a solution for my use case.

You can interpret the GeoTransform as the following:

• X (or easting) coordinate of the image
• meters (or other unit) going east per pixel in x axis of image
• meters (or other unit) going north per pixel in x axis of image
• Y (or northing) coordinate of the image
• meters (or other unit) going east per pixel in y axis of image
• meters (or other unit) going north per pixel in y axis of image

So, if your image was aligned to the SRS and each pixel represents an area of 0.1 meters to east and 0.2 meters to north and is located at coordinate (100, 200), then the Geotransform has the following shape:

• 100 - the image is 100 meters east of the origin of the SRS
• 0.1 - each pixel goes 0.1 meters to the east on x axis
• 0.0 - each pixel doesn't go north on the x axis
• 200 - the image is 300 meters north of the origin of the SRS
• 0.0 - each pixel doesn't go east on the y axis
• 0.2 - each pixel goes 0.2 meters to the borth on y axis

Now, if the image is rotated, in your case 120 degrees, the following calculation can be applied.

from math import cos, sin, radians
rotation = 120.0
ulx = 100
uly = 200
xres = cos(radians(rotation)) * 0.1
yres = cos(radians(rotation - 90)) * 0.2
xskew = sin(radians(rotation)) * 0.1
yskew = sin(radians(rotation - 90)) * 0.2

Please be aware, that this only works on Cartesian SRS. Consider to warp the image into something like UTM first.

At least, it worked in my (rather simple) use case.

• seb007: I think that there is a mistake in the order of the transform components. For example, the third component of the transform GT(2) is: meters (or other unit) going EAST per pixel in Y axis of image. See gdal.org/en/latest/tutorials/geotransforms_tut.html