0

我已经将很多 WGS84 坐标（我知道它们存在于我的栅格数据中）转换为 UTM 并将它们插入我的程序只是让它告诉我它们超出了范围。我的栅格是 4695x9798，我不确定为什么我的坐标总是落在那个窗口之外

import numpy as np
from osgeo import gdal,ogr
import struct


gdata = gdal.Open('sinusoidal.tif')
geot = gdata.GetGeoTransform()


x = (284905 - geot[0])/geot[1]
y = (5936117 - geot[3])/(geot[5])

myarray = np.array(gdata.GetRasterBand(1).ReadAsArray())


print gdata.RasterXSize
print gdata.RasterYSize

rb = gdata.GetRasterBand(1)
intval = rb.ReadAsArray(x,y,1,1)
print intval

错误消息： RasterIO() 中的访问窗口超出范围。在 4695x9798 的栅格上请求大小为 1x1 的 (6126,1437)。

4

1 回答 1

1

错误描述非常明确。您正在请求一个超出栅格范围的像素。这可能与您提供的 UTM 坐标有关，或者与您未考虑的地理变换的某些方面有关（xskew 或 yskew）。获取行列像素索引的更规范方法是使用逆地理变换。

#...
rb = gdata.GetRasterBand(1)
geot = gdata.GetGeoTransform() # maps i,j to x,y
# try-except block to handle different output of InvGeoTransform with gdal versions
try:
    inv_gt_success, inverse_gt = gdal.InvGeoTransform(geot) # maps x,y to i,j
except:
    inverse_gt = gdal.InvGeoTransform(geot) # maps x,y to i,j
x_utm = 284905
y_utm = 5936117
pix_x = int(inverse_gt[0] + inverse_gt[1] * x_utm +
            inverse_gt[2] * y_utm)
pix_y = int(inverse_gt[3] + inverse_gt[4] * y_utm +
            inverse_gt[5] * y_utm)
val = rb.ReadAsArray(pix_x, pix_y, 1, 1)[0,0]
于 2016-12-22T16:08:23.567 回答