Why are my X and Y values always out of range of my raster file size?

Question

I have converted quite a few WGS84 coordinates, that I know exist within my raster data, to UTM and have plugged them into my program only to have it tell me they are out of range. My raster is 4695x9798 and I'm not sure why my coordinates keep falling outside of that window

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

Error Message: Access window out of range in RasterIO(). Requested (6126,1437) of size 1x1 on raster of 4695x9798.

Answer 1

The error description is very explicit. You are requesting a pixel that is outside of your raster extent. This may be related to the UTM coordinates you are supplying, or some aspect of the geotransform you are not taking into consideration (xskew or yskew). The more canonical way to get the row-col pixel indices is to use the inverse geotransform.

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