[Gdal-dev] Gdal + python for map algebra script?

[Christopher Barker]

Python should work well for this, but I do have question:

> My basic idea:  (1) delimit a single band raster file with a vector
> polygon, 

What does it mean to delimit a raster file with a polygon?

Anyway, you'll probably want to use numpy to help with this kind of 
number crunching. I'm not sure the status of the gdal Python bindings 
and numpy, but you can get the raster data into a numpy array one way or 
another.

Here's some example code that might help:

This puts the RBG data in a numpy array
        ## XY is the coords of the corner of the part of the image you
        ## want, (0,0) if you want the whole thing
        ## WH is the size of the piece of the image you want:
        ##   (Width, Height)
        ## NewWH is the size you want the output image at:
        ##   (Width, Height)
        DataRect = N.zeros((NewWH[0],NewWH[1],3), N.uint8)
        for i in range(3):
             Band = self.Dataset.GetRasterBand(i+1)
             color = Band.ReadRaster( XY[0], # nXOff
                                      XY[1], # nYOff
                                      WH[0], # nXSize
                                      WH[1], # nYSize
                                      NewWH[0], #nBufXSize
                                      NewWH[1], #nBufYSize
                                      GDT_Byte, # data type
                                      )
             # put it in the numpy array:
             DataRect[:,:,i] = N.fromstring(color,N.uint8).reshape(NewWH)


-Chris

-- 
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov