5

Is there a way with GDAL (Python) to set a spatial filter (bounding-box) on a raster before I process it? I want to read it as an array.

BandReadAsArray(band, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None, buf_ysize=None, buf_obj=None)
5

You can do it by converting your bounding box to pixel offset and size:

from osgeo import gdal

def map_to_pixel(mx,my,gt):
    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    return px,py

def extent_to_offset(xmin,ymin,xmax,ymax,gt):
    pxmin,pymin = map_to_pixel(xmin,ymin,gt)
    pxmax,pymax = map_to_pixel(xmax,ymax,gt)

    return pxmin,pymin,pxmax-pxmin,pymax-pymin

xmin,ymin,xmax,ymax = [123,234,345,456]
ds=gdal.Open('some.tif')
gt = ds.GetGeoTransform()
xoff, yoff, xsize, ysize = extent_to_offset(xmin,ymin,xmax,ymax,gt)
band=ds.GetRasterBand(1)
array=band.ReadAsArray(xoff, yoff, xsize, ysize)
  • that look good. I'll give it a try. However in the bounding box you have twice 'ymax' – ustroetz May 24 '13 at 16:56

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.