Pourquoi mes valeurs X et Y sont-elles toujours hors de portée de ma taille de fichier raster?

Question

J'ai converti pas mal de coordonnées WGS84, que je connais dans mes données raster, en UTM et les ai branchées dans mon programme uniquement pour me dire qu'elles sont hors de portée. Ma trame est 4695x9798 et je ne sais pas pourquoi mes coordonnées continuent de tomber en dehors de cette fenêtre

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 

Message d'erreur: fenêtre d'accès hors de portée dans RasterIO(). Demandé (6126,1437) de taille 1x1 sur raster de 4695x9798.

Answer 1

La description de l'erreur est très explicite. Vous demandez un pixel en dehors de votre étendue raster. Cela peut être lié aux coordonnées UTM que vous fournissez, ou à certains aspects de la géotransformation que vous ne prenez pas en compte (xskew ou yskew). La façon la plus canonique d'obtenir les indices de pixel de rang-col est d'utiliser la géotransformation inverse.

#... 
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]