Python 中 Lat/Lon Meshgrids 的 GDAL 仿射系数

Question

我在创建新的 GeoTIFF 文件时遇到了一些仿射变换系数问题。我正在做的是对一个科学数据集进行 ETL，该数据集产生一个 2D Ndarray 以及一组包含 Lat 和 Lon 的网格网格 Ndarray。网格网格和数据集数组都具有相同的 645 x 980 尺寸。据我了解，当通过SetGeoTransform() 方法从 Python GDAL 创建时，GeoTIFF 需要一个仿射系数列表。该列表的格式为[xllcorner, xrotation, x_cellsize, yllcorner, yrotation, y_cellsize]。我的处理方法与此处概述的类似：http://adventuresindevelopment.blogspot.com/2008/12/python-gdal-adding-geotiff-meta-data.html

此时我遇到了问题。我分别使用min() 方法为 lat 和 lon 的两个 meshgrid 数组计算 xllcorner 和 yllcorner，我通过应用公式 [max-min]/dimension size 手动计算 x 和 y 单元格大小，x 维度大小为 x 轴大小对于 lons 网格网格，y 尺寸大小是 lats 网格网格的 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 的左下角。当将数组作为打印矩阵与映射栅格进行比较时，它只会颠倒过来。