Home Articles

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

Asked
Viewed 962 times
1

在创建新的GeoTIFF文件时,我遇到了一些Affine转换系数问题 . 我正在做的是在一个科学数据集上的ETL,它产生一个2D Ndarray以及一组包含Lat和Lon的meshgrid Ndarray . meshgrids和数据集数组都具有645 x 980的相同尺寸 . 据我所知,GeoTIFF需要一个Affine系数列表,当通过 SetGeoTransform() 方法从Python GDAL创建时 . 该列表的形式为 [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维度大小为lons meshgrid的x轴大小并且y维度大小是lats meshgrid的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 . 所以我的问题是如何使用lat / lon meshgrids和共享公共维度的数据数组正确创建Affine系数?我认为我的细胞大小计算不能简单地说是纬度/经度差异;但我不确定 .

1 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的左下角 . 将阵列作为打印矩阵与映射栅格进行比较时,它只会颠倒过来 .

Related