Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

1 Answer 1

up vote 2 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]
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)
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

 
discard

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.