Почему мои значения X и Y всегда выходят за пределы размера моего растрового файла?

Я преобразовал довольно много координат 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(). Запрошено (6126,1437) размера 1x1 на растре 4695x9798.


python-2.7 arcgis gdal osgeo
person Jb11281992    schedule 21.12.2016    source источник
comment
Вы добавили тег arcgis, но какое отношение эта проблема имеет к ArcGIS? Вы, кажется, не используете его, основываясь на вашем импорте.   -  person Gary Sheppard    schedule 21.12.2016

Ответы (1)


arrow_upward
0
arrow_downward

Описание ошибки очень подробное. Вы запрашиваете пиксель, который находится за пределами экстента вашего растра. Это может быть связано с координатами 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]
person Logan Byers    schedule 22.12.2016