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 に答える