How can I cut a dem by longitude and latitude in python gdal?

Question

I am wondering if I can cut a srtm dem data(.tif) into smaller parts in python.

We know that we can get longitude and latitude by

# importing package
from osgeo import gdal

# load tiff data
dataset=gdal.Open("srtm_input.tif")

# transformation data
im_geotrans = dataset.GetGeoTransform()

# calcualte boundaries
minx = im_geotrans[0]
miny = im_geotrans[3] + im_width*im_geotrans[4] + im_height*im_geotrans[5]
maxx = im_geotrans[0] + im_width*im_geotrans[1] + im_height*im_geotrans[2]
maxy = im_geotrans[3]

The minx, miny and maxx and maxy will gives the borders in longitude and latitude.

My question is: If we provide four boundaries in forms of minx, miny and maxx and maxy (that is, longitude and latitude), can I use osgeo to cut the dem data?

I may not explaining the question clearly, so I am adding a image here.

If I already have the DEM data of the following area in U.S., and I want to get the red shaded part with exact longitude and latitude boundaries(the red lines). Can I manage to do that with osgeo.gdal?

Answer 1

Yes you can do this quite easily. Only requirement is that you keep the coordinate system consistent throughout the process. (e.g. your bounds coordinates should be in the same coordinate system than your srtm dataset) You can use gdal translate

I give an example using GeoPandas aswell to get the bounds coordinates

from osgeo import gdal
import geopandas as gpd

# e.g. srtm.tif coordinates epsg:4326
dataset = gdal.Open("srtm_input.tif")

# path to where you want the clipped raster
outputSrtm = "path/to/your/clipped/srtm.tif"

# e.g. shapefile.shp coordinates epsg:4326
minx, miny, maxx, maxy = gpd.read_file("path/to/shapefile.shp").unary_union.bounds

ds = gdal.Translate(outputSrtm , dataset, 
                    projWin = [minx, maxy,
                               maxx, miny]) # OR [ulx, uly, lrx, lry]
ds = None