GDAL SetGeoTransform参数描述

Question

有人能否帮我解决SetGeoTransform的参数问题？我正在使用GDAL创建栅格图层，但找不到SetGeoTransform的第3和第5个参数的说明。应该是单元格x轴和y轴的定义。我试图在这里和这里找到相关信息，但没有任何结果。

我需要找到这两个参数的说明...它们是以度、弧度、米或其他方式表示的值？

- david_p

Answer 1

地理转换（geotransform）是使用仿射变换将地图坐标系转换为像素坐标系和相反方向的过程。第三个和第五个参数（与第二个和第四个参数一起）用于定义旋转，如果您的图像没有“北对齐”。

但是大多数图像都是“北对齐”的，那么第三个和第五个参数都为零。

仿射变换包含六个系数，由GDALDataset :: GetGeoTransform() 返回，它们使用以下关系将像素/行坐标映射到地理参考空间中：

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

请参阅仿射地理变换部分： https://gdal.org/tutorials/geotransforms_tut.html

- Rutger Kassies

---

谢谢！但是如果我不想让栅格朝北呢？请看我的问题http://gis.stackexchange.com/questions/123532/slanted-raster-along-the-line-wind-erosion-model-and-windbreak-effectivity。能否以这种方式处理栅格创建？ - david_p

是的，这正是你应该处理这种情况的方式，你需要为每个x像素计算y偏移量，以及为每个y像素计算x偏移量。 - Rutger Kassies

对于你在 GIS SE 中提出的“倾斜”问题，你应该能够使用类似 GT[1] * sin(angle) 的方法来计算 y 偏移量（因此是 GT[4] 或第五个参数）。GT[2] 可以保持为零。但一定要仔细检查！ - Rutger Kassies

我会在不久的将来进行一些测试，然后发布我的研究结果。主要问题是我需要用两种不同的方式来实现“倾斜”栅格。第一种是，即使线条不正交，我也需要将光栅“捕捉”到线条上。第二种方式是，我必须通过风向“倾斜”栅格。希望能够通过您的建议找到解决方法。谢谢！ - david_p

Answer 2

我做了以下代码。

结果我能够使用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

Answer 3

根据上述 gdal datamodel docs提供的信息，可以通过已知"geo"和"pixel"坐标空间中两个控制点（p1、p2）的已知x和y来计算SatGeoTransform的第三个参数x_skew和第五个参数y_skew。在像素空间中，p1应该在p2的左上方。

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

简而言之，这是地理空间点之间欧几里得距离与像素空间图像高度（或宽度）之比。

参数的单位为“地理长度/像素长度”。

以下是使用存储为控制点（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,
])

- 7yl4r

---

数据模型的链接已经失效，这个链接看起来像 这个链接 可以作为一个很好的替代品。 - olepinto