Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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)
share|improve this question
up vote 3 down vote accepted

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]
gt = ds.GetGeoTransform()
xoff, yoff, xsize, ysize = extent_to_offset(xmin,ymin,xmax,ymax,gt)
array=band.ReadAsArray(xoff, yoff, xsize, ysize)
share|improve this answer
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
@user1738154 Thanks, typo fixed now ;) – Luke May 24 '13 at 21:51

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.