DEM提取流域特征的代码

DEM（数字高程模型）用于地理信息系统（GIS）中，提取流域特征通常涉及到流线追踪和区域划分。以下是一个简单的Python示例，使用GDAL库（Geospatial Data Abstraction Library）和FlowDirector模块来提取流域：

from osgeo import gdal, ogr, gdal_array
import numpy as np

# 读取DEM数据
dem_raster = gdal.Open("path_to_your_DEM.tif")
band = dem_raster.GetRasterBand(1)
dem_data = band.ReadAsArray()

# 创建流网络
flow_dir = FlowDirectorSteepest(dem_data, -9999)  # 将未指定方向设为-9999 (通常是海洋或者无效区域)
flow_dir.run_one_step()

# 提取流域
basins = flow_dir.reach_code
catchments = basins.reshape((dem_data.shape[0], dem_data.shape[1]))

# 计算流域面积
cell_size = dem_raster.GetGeoTransform()[1]  # 获取栅格单元尺寸
areas = cell_size * cell_size * np.count_nonzero(catchments, axis=(0, 1))

# 输出流域特征（如最大流域、最小流域、平均流域面积）
max_area_index = np.argmax(areas)
min_area_index = np.argmin(areas)

print(f"最大流域面积: {areas[max_area_index]} 平方米")
print(f"最小流域面积: {areas[min_area_index]} 平方米")

# 存储流域分块
ogr_driver = ogr.GetDriverByName('ESRI Shapefile')
out_shp = ogr_driver.CreateDataSource('output.shp')
out_layer = out_shp.CreateLayer('catchments', geom_type=ogr.wkbPolygon)
field_defn = ogr.FieldDefn('Area', ogr.OFTReal)
out_layer.CreateField(field_defn)
for i, basin in enumerate(np.unique(catchments)):
    mask = catchments == basin
    polygon = ogr.Geometry(ogr.wkbPolygon)
    for ring in ogr.GeometryFromArray(mask).GetGeometryRef(0).ExportToWkt().split(";"):
        if "LINESTRING" not in ring:
            polygon.AddGeometryDirectly(ogr.CreateGeometryFromWkt(ring))
    feature = ogr.Feature(out_layer.GetLayerDefn())
    feature.SetField('Area', areas[i])
    feature.SetGeometry(polygon)
    out_layer.CreateFeature(feature)