Python中从经纬度网格获取GDAL仿射系数

Question

我在创建新的GeoTIFF文件时遇到了一些关于仿射变换系数的问题。我正在对一个科学数据集进行ETL处理，结果是生成一个二维的Ndarray，以及一组包含经纬度的网格Ndarrays。这些网格和数据集的数组大小都是645 x 980。从我理解的情况来看，GeoTIFF在通过Python的GDAL库创建时，需要使用SetGeoTransform()方法提供一组仿射系数。这组系数的格式是[xllcorner, xrotation, x_cellsize, yllcorner, yrotation, y_cellsize]。我的做法和这里描述的类似：http://adventuresindevelopment.blogspot.com/2008/12/python-gdal-adding-geotiff-meta-data.html

在这个过程中，我遇到了问题。我使用min()方法分别计算了两个网格数组（经纬度）的xllcorner和yllcorner，并且通过公式[max-min]/dimension size手动计算了x和y的单元格大小，其中x的维度大小是经度网格的x轴大小，y的维度大小是纬度网格的y轴大小。当我应用这个方法并尝试通过GetRasterBand().WriteArray()写出数组时，出现了这个错误信息：

Traceback (most recent call last):
    ...
    raise ValueError("array larger than output file, or offset off edge")
ValueError: array larger than output file, or offset off edge

因此，我认为我的仿射系数可能计算错了，但根据数据来看，这让我很困惑。我甚至确保在尝试创建仿射系数之前，空间参考系统设置为WGS:84。那么，我的问题是如何正确地使用经纬度网格和共享相同维度的数据数组来创建仿射系数？我觉得我的单元格大小计算不能仅仅是经纬度的差值，但我不太确定。

Answer 1

这个错误通常是因为你期待的数组形状和实际的不匹配。比如，你可以用下面的代码查看一下你期待的形状是什么：

band = src.GetRasterBand(1)
arr = band.ReadAsArray()
print(arr.shape)  # (656L, 515L)

这段代码需要写成的numpy数组的形状是这样的：

assert other_array.shape == arr.shape
band.WriteArray(other_array)

如果你想引发同样的ValueError错误，可以把数组的形状改得在某个维度上更长，比如：

band.WriteArray(other_array.T)

---

至于仿射变换，这个通常不会引发错误，因为它通常只是作为数据存储。地理信息系统（GIS）的栅格图像通常会把世界坐标放在左上角，并使用一个负的dy值来向下计算行数。不过，使用左下角和正的dy值在大多数软件中也是可以的。这样做的话，比较打印出来的矩阵和映射的栅格图像时，可能会出现上下颠倒的情况。