3

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]
  • I suggest to think if pixel count is zero or one based, and if you have considered that for GDAL the anchor point is in the top left corner of the pixel. They may be OK in your code, just considering what could go wrong. – user30184 Oct 11 '16 at 21:31
2

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.

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.