Plot——Filling

 

 

太极

 

# -*- coding: utf-8 -*- """ Python 一体化河流—流域—通量计算示例 1) 用 WhiteboxTools 填洼 2) 用 PySheds 计算流向 & 汇流累积 3) 用 GeoPandas 找出河—海交点,提取流域 4) 计算径流体积 & 加权累积 5) 采样出口点 & 计算通量 """ import os import rasterio import numpy as np from whitebox import WhiteboxTools from pysheds.grid import Grid import geopandas as gpd # —— 1. 数据路径 —— # DEM_PATH = r"D:\JW\BaiduSyncdisk\数据\DYW\GIS数据\高程DEM\高程WGS84\深圳市-谷歌高程18级\深圳市-谷歌高程18级_高程\L18\深圳市-谷歌高程18级.tif" RAIN_PATH = r"C:\Users\ASUS\AppData\Local\Temp\processing_JnOHDT\6f39273e4b524aa5b6480838ce7c80f0\OUTPUT.tif" RIVER_SHP = r"D:\JW\BaiduSyncdisk\数据\DYW\GIS数据\SHP\深圳市(WGS84)21下载20入库\深圳市_河流.shp" COAST_SHP = r"D:\JW\ArcGIS\liu_10m_coastline\10m.shp" OUT_DIR = r"D:\JW\BaiduSyncdisk\数据\SZML\全国站点逐日气象数据\气候数据GIS" os.makedirs(OUT_DIR, exist_ok=True) # —— 2. 用 Whitebox 填洼 —— # wbt = WhiteboxTools() filled_dem = os.path.join(OUT_DIR, "dem_filled.tif") # 使用 Whitebox 填洼算法 wbt.fill_depressions( dem=DEM_PATH, output=filled_dem ) print("DEM 填洼完成:", filled_dem) # —— 3. 用 PySheds 计算流向 & 汇流累积 —— # grid = Grid.from_raster(filled_dem, data_name='dem') # 计算流向:传入 Raster 对象 grid.flowdir(dem='dem', out_name='dir', dirmap='d8') grid.accumulation(direction='dir', out_name='acc') flow_dir_tif = os.path.join(OUT_DIR, "flow_direction.tif") flow_acc_tif = os.path.join(OUT_DIR, "flow_accumulation.tif") grid.to_raster(flow_dir_tif, data='dir') grid.to_raster(flow_acc_tif, data='acc') print("流向 & 累积输出:", flow_dir_tif, flow_acc_tif) # —— 4. 找河—海交点 & 提取流域 —— # rivers = gpd.read_file(RIVER_SHP).to_crs("EPSG:4326") coast = gpd.read_file(COAST_SHP).to_crs("EPSG:4326") intsct = gpd.overlay(rivers, coast, how='intersection') if intsct.empty: raise RuntimeError("未在河网与海岸线上找到交点,请检查 SHP。") outlet_pt = intsct.geometry.iloc[0] x0, y0 = outlet_pt.x, outlet_pt.y # 提取流域 basin = grid.catchment(x=x0, y=y0, data='dir', out_name='catch') basin_tif = os.path.join(OUT_DIR, "basin.tif") grid.to_raster(basin_tif, data='catch') print("流域提取完成:", basin_tif) # —— 5. 计算径流体积栅格 & 加权累积 —— # with rasterio.open(RAIN_PATH) as src: rain = src.read(1) meta = src.meta.copy() Cr = 0.5 # 径流系数 # 像元面积 (m²) pix_area = abs(src.transform.a * src.transform.e) # 体积 Vi = rain(mm)*1e-3(m/mm)*pix_area runoff_vol = rain * Cr * 1e-3 * pix_area # 加权累积:传入 Raster 对象和权重数组 grid.accumulation(grid.dir, out_name='flow_acc_vol', weights=runoff_vol) flow_acc_vol_tif = os.path.join(OUT_DIR, "flow_acc_volume.tif") grid.to_raster(flow_acc_vol_tif, data='flow_acc_vol') print("加权累积(体积)输出:", flow_acc_vol_tif) # —— 6. 采样出口 & 计算通量 —— # # 将出口坐标转换为像素索引 tuple_rc = grid.raster_to_pixel(x0, y0) row, col = tuple_rc V_total = grid.view('flow_acc_vol')[row, col] # 单位 m³ Q_m3s = V_total / 86400.0 print(f"出口累积总径流体积 V = {V_total:.2f} m³") print(f"日平均通量 Q = {Q_m3s:.3f} m³/s")
07-19
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值