[Gdal-dev] raster query

[Bryan Keith]

John,

I'm attaching a Python routine that queries a given band of a gdal 
raster at a given x and y and returns the interpolated value.  You'll 
have to modify it a bit if you want the value of the raster without 
interpolation.

I just updated it to work with rasters that have different cell size in 
i and j.  Hope it works and helps.

Bryan

Bryan Keith
GIS Specialist
Geomega, Inc.
Boulder, CO, USA

John Cartwright wrote:
> Hello All,
>
> does anyone know of a existing gdal tool for querying a raster?  I'd 
> like to get the raster value from a file for a given coordinate.
>
> Thanks!
>
> -- john
> _______________________________________________
> Gdal-dev mailing list
> Gdal-dev at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/gdal-dev
>
>
-------------- next part --------------
def PointValue(dataset, intBand, dblX, dblY):
  #this routine returns the cellvalue on the band at (dblX, dblY) using bilinear interpolation
  band = dataset.GetRasterBand(intBand)
  dblNDV = band.GetNoDataValue()
  lngXCount = band.XSize
  lngYCount = band.YSize
  geotransform = dataset.GetGeoTransform()
  dblXMin = geotransform[0]
  dblYMax = geotransform[3]
  dblCellSizeX = geotransform[1]
	dblCellSizeY = geotransform[5]
  dblXLoc = (dblX - dblXMin) / dblCellSizeX
  lngXLoc = int(dblXLoc)
  dblYLoc = (dblYMax - dblY) / dblCellSizeY
  lngYLoc = int(dblYLoc)
  if ((dblXLoc < 0.5) or (dblXLoc > lngXCount - 0.5) or (dblYLoc < 0.5) or \
    (dblYLoc > lngYCount - 0.5)):
    return dblNDV
  if ((dblXLoc - lngXLoc) < 0.5):
    lngXLoc = lngXLoc - 1
  if ((dblYLoc - lngYLoc) < 0.5):
    lngYLoc = lngYLoc - 1
  dblXPercent = ((dblXLoc - lngXLoc) + 0.5) - (int((dblXLoc - lngXLoc) + 0.5))
  dblYPercent = ((dblYLoc - lngYLoc) + 0.5) - (int((dblYLoc - lngYLoc) + 0.5))
  strRaster = band.ReadRaster(lngXLoc,lngYLoc,2,2,2,2,GDT_Float32)
  dblValue = struct.unpack('f' * 4,strRaster)
  if ((dblValue[0] == dblNDV) or (dblValue[1] == dblNDV) or (dblValue[2] == dblNDV) \
    or (dblValue[3] == dblNDV)):
    return dblNDV
  dblX1 = dblValue[0] + ((dblValue[1] - dblValue[0]) * dblXPercent)
  dblX2 = dblValue[2] + ((dblValue[3] - dblValue[2]) * dblXPercent)
  dblValue = ((dblX2 - dblX1) * dblYPercent) + dblX1
  #print dblX, dblY, dblValue
  return dblValue