How to get the pixels from a geotiff file in gdal python for a given bbox?

Question

I have a GeoTIFF file in WSG 84 Projection reference (it's a GeoTIFF of SRTM). I can use python & gdal to read the whole lot of data into an array (with ReadAsArray()), but I want to only have a subset of the data, everything inside one bounding box (defined by a top latitude, bottom latitude, right longitude, left longitude). But I'm very new to gdal, and I can't figure out how to do it.

Here is some python code to confirm it is in WSG 84, etc.

>>> from osgeo import gdal
>>> gg = gdal.Open('srtm_34_02/srtm_34_02.tif')
>>> gg.GetGeoTransform()
(-15.000416327675339, 0.0008333333333333334, 0.0, 55.00041712672993, 0.0, -0.0008333333333333334)
>>> gg.GetProjectionRef()
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]'

It doesn't matter if the resultant array of data includes points that are just outside the bbox, due to rounding. Just so long as it doesn't exclude any points that are in the bbox, and it's as close enough as is practical.

Answer 1

I recommend the use of a vrt file, a virtual raster file that would point to your data with the bounding box that you want.

you can call it from python using

import subprocess
subprocess.call(["gdalbuildvrt", "-te", "xmin", "ymin", "xmax", "ymax", "output.vrt", "input.tif"]) 

then you use your vrt like an image.

If you want you can create the vrt with the binding, but this is notrecommended (see here)

Answer 2

A (late) python-only solution with gdal.

ReadAsArray allows to read subsets only. Consider that your bbox is defined by xmin, xmax, ymin, ymax. Then, something like that should work:

gdalRaster = gdal.Open(rasterFileName)

geoTransform = gdalRaster.GetGeoTransform()
x0 = geoTransform[0]
y0 = geoTransform[3]
dx = geoTransform[1]
dy = geoTransform[5]    
    
idXStart = int(math.floor((xmin-x0) / dx))
idYStart = int(math.floor((ymin-y0) / dy))
idXEnd = int(math.ceil((xmax-x0) / dx))
idYEnd = int(math.ceil((ymax-y0) / dy))

rasterBand = gdalRaster.GetRasterBand(1)  
subsetData = rasterBand.ReadAsArray(idXStart , idYStart , idXEnd - idXStart , idYEnd - idYStart )