9

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?

Glorfindel
  • 21,988
  • 13
  • 81
  • 109
david_p
  • 223
  • 1
  • 6
  • 17

3 Answers3

15

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

Vitaly Olegovitch
  • 3,509
  • 6
  • 33
  • 49
Rutger Kassies
  • 61,630
  • 17
  • 112
  • 97
  • Thanks! But what if I don't want raster to be north oriented? See my question http://gis.stackexchange.com/questions/123532/slanted-raster-along-the-line-wind-erosion-model-and-windbreak-effectivity Is it posible to hande raster creation this way? – david_p Nov 27 '14 at 11:06
  • 1
    Yes that's exactly how you should handle such a case, you need to calculate the y-offset for each x-pixel, and the x-offset for each y-pixel. – Rutger Kassies Nov 27 '14 at 11:30
  • For a 'skew' like in your question at the `GIS SE` you should be able to calculate the y-offset (so GT[4], or the fifth param), with something like `GT[1] * sin(angle)`. GT[2], could be left zero. But check it carefully! – Rutger Kassies Nov 27 '14 at 13:25
  • I will make some tests in near future and then I post what I found. The main problem is that I need 'skew' raster in two different ways. The first is that I need the raster to be 'snapped' od line, even if the line wont be orthogonal. And the second way, I have to 'skew' raster by the wind direction.. I hope I could find the way with your advices. Thanks! – david_p Nov 27 '14 at 14:55
1

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,
])
Glorfindel
  • 21,988
  • 13
  • 81
  • 109
7yl4r
  • 4,788
  • 4
  • 34
  • 46
  • 3
    The link about the data model is broken, it looks like [this would be a good replacement](https://gdal.org/user/raster_data_model.html#affine-geotransform). – olepinto Feb 18 '21 at 17:26
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
shinya
  • 21
  • 2