Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

pixel/array position to lat long gdal Python

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))
like image 466
niaR Avatar asked Sep 21 '18 12:09

niaR


2 Answers

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)
like image 177
t_prz Avatar answered Sep 22 '22 04:09

t_prz


There is no reason to do this manually. You just need to install rasterio, a python library based on GDAL and numpy. Even if the image has an offset/rotation rasterio will take care of it.

You just have to do:

import rasterio

with rasterio.open('rasters/raster.tif') as map_layer:
    coords2pixels = map_layer.index(235059.32,810006.31) #input lon,lat
    pixels2coords = map_layer.xy(500,598)  #input px, py

Both functions return a tuple of the pixels and coordinates respectively.

like image 36
mobiuscreek Avatar answered Sep 19 '22 04:09

mobiuscreek