如何在python gdal中通过经度和纬度切割dem？

问题描述

我想知道是否可以在 python 中将 srtm dem 数据（.tif）切割成更小的部分。

我们知道我们可以通过以下方式获得经度和纬度

# 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]

minx、miny 和 maxx 和 maxy 将给出经度和纬度的边界。

我的问题是：如果我们以minx、miny和maxx和maxy（即经度和纬度）的形式提供四个边界，我可以使用osgeo 来削减 dem 数据？

我可能没有把问题解释清楚，所以我在这里添加了一张图片。

如果我已经拥有美国以下地区的 DEM 数据，并且我想获得具有精确经纬度边界（红线）的红色阴影部分。我可以通过 osgeo.gdal 做到这一点吗？

解决方法

是的，您可以很容易地做到这一点。 唯一的要求是在整个过程中保持坐标系一致。 （例如，您的边界坐标应与您的 srtm 数据集位于同一坐标系中） 您可以使用 gdal 翻译

我还举了一个使用 GeoPandas 获取边界坐标的示例

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,miny]) # OR [ulx,uly,lrx,lry]
ds = None