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)
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)
