[FWTools] random access querying raster image - help

[Bryan Keith]

Kenton,

Below I pasted some Python code.  It will return the cell value of the 
cell where the xy is found.  It will not do any interpolation.  I have a 
separate routine for that.  I pulled this bit of code out of a bigger 
file.  I may be missing another import statement.

Bryan

from gdalconst import *
import struct
import gdal
# 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)
   if gdal.GetDataTypeName(band.DataType) == 'Int16':
     dblValue = struct.unpack('h',strRaster)
   elif gdal.GetDataTypeName(band.DataType) == 'Float32':
     dblValue = struct.unpack('f',strRaster)
   else:
     print 'unrecognized DataType: ',gdal.GetDataTypeName(band.DataType)
     return ndv

   return dblValue[0]

Kenton Williams wrote:
> Hello,
> 
> I am needing to write a simple (I think) program/script that takes a
> coordinate pair as input, queries a raster at those coordinates, and
> returns the raster's value at that point.  I would prefer to use
> Python, but I am open to better suggestions. Is this possible to do
> using any of FWtools or GDAL's utility programs?  If it's possible,
> could someone write a very brief list of the sequence of suggested
> steps?  Any help would be greatly appreciated.  Thanks in advance.
> 
> Kenton W.
> 
> _______________________________________________ FWTools mailing list 
> FWTools at lists.maptools.org 
> http://lists.maptools.org/mailman/listinfo/fwtools 
> http://fwtools.maptools.org/
> 
>