[GIS] Get information from geotiff with lat/lon

Question

I have a geotiff file containing data like temperature, now I need to extract the temperature if a real-world lat/lon pair is given. I am using python and GDAL, following is a code I found and it worked.

However, when I used QGIS (I used value tool plugin) to validate the result, for a given lon/lat, the output value from QGIS is not consistent with the following code. The cell of my code is actually some 1- or 2-distance neighbors of what QGIS looks into. I don't see any problem of my code, is QGIS value tool has some precision issues that I should not trust?

from osgeo import gdal,ogr
from osgeo.gdalconst import *
import struct
import sys

lat = float(sys.argv[2])
lon = float(sys.argv[3])

def pt2fmt(pt):
    fmttypes = {
        GDT_Byte: 'B',
        GDT_Int16: 'h',
        GDT_UInt16: 'H',
        GDT_Int32: 'i',
        GDT_UInt32: 'I',
        GDT_Float32: 'f',
        GDT_Float64: 'f'
        }
    return fmttypes.get(pt, 'x')

ds = gdal.Open(sys.argv[1], GA_ReadOnly)
if ds is None:
    print 'Failed open file'
    sys.exit(1)

transf = ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount #1
band = ds.GetRasterBand(1)
bandtype = gdal.GetDataTypeName(band.DataType) #Int16
driver = ds.GetDriver().LongName #'GeoTIFF'
transfInv = gdal.InvGeoTransform(transf)
px = (lon-transf[0])/transf[1]
py = (lat-transf[3])/transf[5]
#px, py = gdal.ApplyGeoTransform(transfInv, lon, lat)
structval = band.ReadRaster(int(px), int(py), 1, 1,  buf_type = band.DataType )
fmt = pt2fmt(band.DataType)
intval = struct.unpack(fmt , structval)
print "value: ", intval[0]

Answer 1

I modified your code slightly to print more than one value for a very small (clipped) version of a real raster:

from osgeo import gdal,ogr
from osgeo.gdalconst import *
import struct
import sys

#lat = float(sys.argv[2])
#lon = float(sys.argv[3])

def pt2fmt(pt):
    fmttypes = {
        GDT_Byte: 'B',
        GDT_Int16: 'h',
        GDT_UInt16: 'H',
        GDT_Int32: 'i',
        GDT_UInt32: 'I',
        GDT_Float32: 'f',
        GDT_Float64: 'f'
        }
    return fmttypes.get(pt, 'x')

my_path = "/home/zeito/pyqgis_data/test_long_lat.tif"

#ds = gdal.Open(sys.argv[1], GA_ReadOnly)
ds = gdal.Open(my_path, GA_ReadOnly)

if ds is None:
    print 'Failed open file'
#    sys.exit(1)

transf = ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount #1
band = ds.GetRasterBand(1)
bandtype = gdal.GetDataTypeName(band.DataType) #Int16
driver = ds.GetDriver().LongName #'GeoTIFF'
transfInv = gdal.InvGeoTransform(transf)

#(3,4), (5,2), (6,8)
p1 = [(16.4851469606, 44.5051091348), 
      (16.4672665386, 44.4863346918),
      (16.5215038186, 44.4782885019)]

for item in p1:
    px = (item[0]-transf[0])/transf[1]
    py = (item[1]-transf[3])/transf[5]
    #px, py = gdal.ApplyGeoTransform(transfInv, lon, lat)
    structval = band.ReadRaster(int(px), int(py), 1, 1,  buf_type = band.DataType )
    fmt = pt2fmt(band.DataType)
    intval = struct.unpack(fmt , structval)
    print "value: ", intval[0]

After running the code at the Python Console of QGIS, it could be observed that it works well (see next image). I corroborated that because I have my own version of Value Tool Plugin (print row and column indices where first pixel is 0,0) and another one that uses indices where first pixel is (1,1). Values printed at Python Console were as expected in all three cases.