Description of parameters of GDAL SetGeoTransform

Question

Can anyone help me with parameters for SetGeoTransform? I'm creating raster layers with GDAL, but I can't find description of 3rd and 5th parameter for SetGeoTransform. It should be definition of x and y axis for cells. I try to find something about it here and here, but nothing.

I need to find description of these two parameters... It's a value in degrees, radians, meters? Or something else?

Answer 1

I did do like below code. As a result I was able to do same with SetGeoTransform.

# new file
dst = gdal.GetDriverByName('GTiff').Create(OUT_PATH, xsize, ysize, band_num, dtype)

# old file
ds = gdal.Open(fpath)
wkt = ds.GetProjection()
gcps = ds.GetGCPs()

dst.SetGCPs(gcps, wkt)
...
dst.FlushCache()
dst = Nonet

Answer 2

The geotransform is used to convert from map to pixel coordinates and back using an affine transformation. The 3rd and 5th parameter are used (together with the 2nd and 4th) to define the rotation if your image doesn't have 'north up'.

But most images are north up, and then both the 3rd and 5th parameter are zero.

The affine transform consists of six coefficients returned by GDALDataset::GetGeoTransform() which map pixel/line coordinates into georeferenced space using the following relationship:

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

See the section on affine geotransform at: https://gdal.org/tutorials/geotransforms_tut.html

Answer 3

Given information from the aforementioned gdal datamodel docs, the 3rd & 5th parameters of SatGeoTransform (x_skew and y_skew respectively) can be calculated from two control points (p1, p2) with known x and y in both "geo" and "pixel" coordinate spaces. p1 should be above-left of p2 in pixelspace.

x_skew = sqrt((p1.geox-p2.geox)**2 + (p1.geoy-p2.geoy)**2) / (p1.pixely - p2.pixely)`
y_skew = sqrt((p1.geox-p2.geox)**2 + (p1.geoy-p2.geoy)**2) / (p1.pixelx - p2.pixelx)`

In short this is the ratio of Euclidean distance between the points in geospace to the height (or width) of the image in pixelspace.

The units of the parameters are "geo"length/"pixel"length.

Here is a demonstration using the corners of the image stored as control points (gcps):

import gdal
from math import sqrt

ds = gdal.Open(fpath)
gcps = ds.GetGCPs()

assert gcps[0].Id == 'UpperLeft'
p1 = gcps[0]
assert gcps[2].Id == 'LowerRight'
p2 = gcps[2]

y_skew = (
    sqrt((p1.GCPX-p2.GCPX)**2 + (p1.GCPY-p2.GCPY)**2) /
    (p1.GCPPixel - p2.GCPPixel)
)
x_skew = (
    sqrt((p1.GCPX-p2.GCPX)**2 + (p1.GCPY-p2.GCPY)**2) /
    (p1.GCPLine - p2.GCPLine)
)

x_res = (p2.GCPX - p1.GCPX) / ds.RasterXSize
y_res = (p2.GCPY - p1.GCPY) / ds.RasterYSize

ds.SetGeoTransform([
    p1.GCPX,
    x_res,
    x_skew,
    p1.GCPY,
    y_skew,
    y_res,
])