I am trying to convert positions in a raster representing a .tif to the corresponding global coordinates. Converting the whole array to a tif and load it to QGIS everything is referenced fine, but using the below calculation method for single points there is a slight offset (to east-north-east in the resulting coordinates....
raster.tif uses ETRS 89 UTM Zone 32N
Does anyone has an idea?
from osgeo import ogr, gdal, osr
import numpy as np
raster = gdal.Open("rasters/raster.tif")
raster_array = np.array(raster.ReadAsArray())
def pixel2coord(x, y):
xoff, a, b, yoff, d, e = raster.GetGeoTransform()
xp = a * x + b * y + xoff
yp = d * x + e * y + yoff
return(xp, yp)
print(pixel2cood(500,598))
I think the issue might be that the xoff and yoff contain coordinates of the upper left corner of the most upper left pixel, and you need to calculate coordinates of the pixel's center.
def pixel2coord(x, y):
xoff, a, b, yoff, d, e = raster.GetGeoTransform()
xp = a * x + b * y + a * 0.5 + b * 0.5 + xoff
yp = d * x + e * y + d * 0.5 + e * 0.5 + yoff
return(xp, yp)