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.

2 Answers 2


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)


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 )
  • this should be accepted
    – cinqS
    Commented Jun 13, 2023 at 11:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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