[gdal-dev] query raster values by xy-location

[Bryan Keith]

Attached is some python code you can use to return the value at a cell
given an xy.  It does not interpolate.  I have a different function for
that if that's what you want.

Bryan

> Hi,
>
> 2008/9/26 <Andruit at gmx.de>
>
>>
>> I have an ESRI ascii raster and would like to get the raster values at
>> defined x y locations.
>> At the moment I am capable to load the Raster Dataset and can perform
>> several operations like querying the  Min/Max values, the origin, the
>> pixelsize and so on....
>>
>> But is there a function to get the cell value of a raster at a certain
>> location (x y)????
>
>
> You can work out the row and column of your pixel from your (x,y) location
> using the raster's geotransform. You can then read in the value of that
> particular pixel.
>
> Hope that helps,
> J
>
>
> --
> Centre for Terrestrial Carbon Dynamics
> Department of Geography, University College London
> Gower Street, London WC1E 6BT, UK
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
# returns the cellvalue on the band at (x, y)
def CellValue(dataset, band, x, y):
  #John Cartwright <John.C.Cartwright at noaa.gov> sent this via e-mail
  #3/23/2006 14:17
  #He modified the PointValue routine that I had sent him
  #print 'point (real-world coords): ',x,',',y
  band = dataset.GetRasterBand(band)
  # set a default NDV if none specified
  if (band.GetNoDataValue() == None):
    band.SetNoDataValue(-9999)
  ndv = band.GetNoDataValue()

  cols = band.XSize
  rows = band.YSize
  #print 'cols,rows: ',cols,',',rows

  geotransform = dataset.GetGeoTransform()
  cellSizeX = geotransform[1]
  # Y-cell resolution is reported as negative
  cellSizeY = -1 * geotransform[5]

  minx = geotransform[0]
  maxy = geotransform[3]
  maxx = minx + (cols * cellSizeX)
  miny = maxy - (rows * cellSizeY)
  #print 'bbox(real-world coords):',minx,',',miny,',',maxx,',',maxy
  if ((x < minx) or (x > maxx) or (y < miny) or (y > maxy)):
    #print 'given point does not fall within grid'
    return ndv

  # calc point location in pixels
  xLoc = (x - minx) / cellSizeX
  xLoc = int(xLoc)
  yLoc = (maxy - y) / cellSizeY
  yLoc = int(yLoc)
  #print 'point (pixels): ',xLoc,',',yLoc

  if ((xLoc < 0.5) or (xLoc > cols - 0.5)):
    #print 'xcoordinate out of bounds'
    return ndv

  if ((yLoc < 0.5) or (yLoc > rows - 0.5)):
    #print 'ycoordinate out of bounds'
    return ndv

  strRaster = band.ReadRaster(xLoc, yLoc, 1, 1, 1, 1, band.DataType)
  sDT = gdal.GetDataTypeName(band.DataType)
  if (sDT == 'Int16'):
    dblValue = struct.unpack('h', strRaster)
  elif (sDT == 'Float32'):
    dblValue = struct.unpack('f', strRaster)
  elif (sDT == 'Byte'):
    dblValue = struct.unpack('B', strRaster)
  else:
    print 'unrecognized DataType:', gdal.GetDataTypeName(band.DataType)
    return ndv

  return dblValue[0]