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,
])
```