GDAL ReadRaster 获取像素块地图范围

Question

GDAL 关于像素行/列与地图坐标之间的地理变换的用户指南说使用以下参数数组：

GT[0] x-coordinate of the upper-left corner of the upper-left pixel.
GT[1] w-e pixel resolution / pixel width.
GT[2] row rotation (typically zero).
GT[3] y-coordinate of the upper-left corner of the upper-left pixel.
GT[4] column rotation (typically zero).
GT[5] n-s pixel resolution / pixel height (negative value for a north-up image)

因此，如果我有一个从左上角坐标向右 10 列、向下 20 行的像素块，即 GT[0]、GT[3]，如何计算 10 x 20 块的范围像素？

我认为像素块的范围将由4个角的地图（x，y）组成，例如：

the upper left vertex (x, y) would be:
  (x_map = GT[0], 
   y_map = GT[3])

the upper right vertex (x, y) would be:
  (x_map = GT[0] + 10 * GT[1] + 20 * GT[2], 
   y_map = GT[3])

the lower right vertex (x, y) would be:
  (x_map = GT[0] + 10 * GT[1] + 20 * GT[2],
   y_map = GT[3] + 10 * GT[4] + 20 * GT[5])

the lower left vertex (x, y) would be:
  (x_map = GT[0],
   y_map = GT[3] + 10 * GT[4] + 20 * GT[5])

这是一个真实的例子： 这看起来还好吗？有没有 GDAL 直接获取这些的方法？如果网格真的很大，这意味着原始的左上角可能离我试图获取范围的像素块很远，误差不会累积在计算的顶点的地图坐标中吗？

Answer 1

您可以使用

ApplyGeoTransform

函数让 GDAL 执行此操作。我正在使用 Python API，但它在 C# 中应该是相同/相似的：

from osgeo import gdal

gt = (-124.5, 9.26e-5, 0, 50.0017, 0, -9.26e-5)

# upper left
ulx, uly = gdal.ApplyGeoTransform(gt, 0, 0)

# lower right
lrx, lry = gdal.ApplyGeoTransform(gt, 10, 20)

其中给出

(-124.5, 50.0017, -124.499074, 49.999848)

作为范围（对于

ulx, uly, lrx, lry

）。使用这个微小的分辨率，上面的值可能会发生一些舍入。

请注意，它给出了像素的左上角位置，而不是中心。要获得中心，您可以添加半个像素（例如

gdal.ApplyGeoTransform(gt, 10+0.5, 20+0.5)

。这与 NetCDF 等格式通常存储坐标的方式不同。

有一个附带的函数可以反转地理变换，它允许您计算反向（映射到像素）。

gt_inv = gdal.InvGeoTransform(gt)
ulx_px, uly_px = gdal.ApplyGeoTransform(gt_inv, lrx, lry)

使用相同的示例给我

(10.00000000023283, 20.0)

（对于

ulx_px, uly_px

）。该函数返回小数像素坐标，因此您可能需要根据应用程序将它们转换为整数。

如果您在地理变换中的分辨率不准确，它确实会累积错误。由于精度有限，这是不可避免的。因此，如果可能的话，谨慎选择范围和分辨率会有所帮助。如果您确实需要一个“干净”范围而不必依赖地理转换，则在元数据中存储已知范围可能是一种解决方法，但这将是非常规的。