地理转换(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我做了以下代码。
结果我能够使用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
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,
])
GIS SE
中提出的“倾斜”问题,你应该能够使用类似GT[1] * sin(angle)
的方法来计算 y 偏移量(因此是 GT[4] 或第五个参数)。GT[2] 可以保持为零。但一定要仔细检查! - Rutger Kassies