[ArcSWAT] 确定水文站的控制区和地表降水量，绘制降水-径流过程图（含完整的 MATLAB 代码）
最编程
2024-06-30 11:35:18
...
ArcSWAT径流模拟结果：绘制降水-径流过程图
- 1 确定水文站控制/集水面积
- 2 根据泰森多边形插值得到面降水量
- 3 绘制降水-径流过程线（MATLAB全代码）
- 参考
下面以洮河流域（甘肃省） 为例，水文站及气象站如下所示：
水文站信息如下：
1 确定水文站控制/集水面积
在SWAT CUP分上下流域率定时，需要查找某水文站上游子流域。
SWAT划分子流域如下所示：
根据各子流域河流流向，确定各水文站集水面积。
2 根据泰森多边形插值得到面降水量
根据泰森多边形插值得到面降水量的详细步骤可参见另一博客-【ArcGIS】基于泰森多边形求流域面降水量。
- 首先根据泰森多边形对整个流域进行插值，结果如下：
- 再根据各水文站控制面积裁剪，得到各水文站控制面积下不同水文站占比，
3 绘制降水-径流过程线（MATLAB全代码）
降水-径流过程线如下所示：（数据为随机拟定，仅供图形展示参考）
模拟径流与实测径流散点图如下：（数据为随机拟定，仅供图形展示参考）
MATLAB代码如下所示：
clc;
clear;
close all;
%% 函数说明-绘制径流模拟过程线
%% 基础设置
pathFigure= '.\Figures\' ;
pathFiles= '.\DataFiles\' ;
%% 导入数据
load('Qo.mat');
[nMonth , nStation] = size(Qo);
load('PmonthYT.mat');
P = PmonthYT( end - nMonth+1 :end , :);
% 模型模拟时间
yearStart = 1976;
yearVali = 2009;
yearEnd = 2020;
nYear = yearEnd - yearStart + 1;
nMonthCali = (yearVali-yearStart)*12;
nMonthVali = (yearEnd-yearVali+1)*12;
%% 计算径流模拟评价结果
% 指标1：相关系数R2（Coefficient of Determination）
% 指标2：纳什效率系数（Nash-Sutcliffe efficiency coefficient, NSE）
% 指标3：偏差百分比（Percent bias, Pbias）
% 整个时期 率定期 验证期
PCC = zeros( nStation, 3);
Pbias = zeros( nStation, 3);
NSE = zeros( nStation, 3);
KGE = zeros( nStation, 3);
R2 = zeros( nStation, 3);
RSR = zeros( nStation, 3);
for in=1:nStation
% 整个时期
% ------------------------------------------
QQo = Qo(:,in);
QQs = Qo(:,in)*0.98;
[PCC( in, 1) , Pbias( in, 1) , NSE( in, 1) , KGE( in, 1) , R2( in, 1) , RSR( in, 1)] = GetModelEfficience( QQs, QQo);
% 率定期
% ------------------------------------------
Qo_cali = Qo(:,in);
Qs_cali = Qo(:,in)*0.98;
[PCC( in, 2) , Pbias( in, 2) , NSE( in, 2) , KGE( in, 2) , R2( in, 2) , RSR( in, 2)] = GetModelEfficience(Qs_cali, Qo_cali);
% 检验期
% ------------------------------------------
Qo_vali = Qo(:,in);
Qs_vali = Qo(:,in)*0.98;
[PCC( in, 3) , Pbias( in, 3) , NSE( in, 3) , KGE( in, 3) , R2( in, 3) , RSR( in, 3)] = GetModelEfficience(Qs_vali, Qo_vali);
end
%% 绘制降水-实测径流/模拟图
for in=1:nStation
QQo = Qo(:,in);
QQm =Qo(:,in)*0.98;
PP = P(:,in);
plotQQP( in, QQm , QQo , PP , yearStart, yearVali , yearEnd, Pbias , NSE , R2);
str= strcat(pathFigure, " ("+Title(in)+") "+Basins(in)+"径流过程线" , '.tiff');
print(gcf, '-dtiff', '-r600', str);
end
%% 绘制模拟径流-观测径流 散点图
for in=1:nStation
QQm = Qo(:,in);
QQo = Qo(:,in);
nData = length( QQm ) ;
XYmax = 1.1*max( QQo );
figure(nStation+in)
hold on;box on;grid off;
h(1) = plot( QQm(1:nMonthCali), QQo(1:nMonthCali),'.','MarkerSize',8,'Color', [0.28235 0.87843 0.8]); % 模拟径流
h(2) = plot( QQm((nMonthCali+1):nData), QQo((nMonthCali+1):nData),'.','MarkerSize',8,'Color', [1 0.25 0]); % 实测径流
h(3) = plot( [ 0.1*max( QQo ) max( QQo ) ], [ 0.1*max( QQo ) max( QQo )] ,'k--', 'linewidth',1.2 );
set(gcf,'color','w'); % 设置当前figure的背景颜色
set(gca,'FontSize',12,'Fontname', 'Times New Roman');
xlabel("\fontname{宋体}模拟径流\fontname{Times New Roman}m^3/s",'FontSize',14);
ylabel("\fontname{宋体}观测径流\fontname{Times New Roman}m^3/s",'FontSize',14);
set(gca,'XLim',[0 XYmax],'YLim',[0 XYmax]); % X轴的数据显示范围
% 整个时期结果
text( 'string',"\itR\rm^2 = "+round( R2( in, 1), 2 ) , 'Units','normalized','position',[0.03,0.92], 'FontSize',12,'FontWeight','normal','FontName','Times New Roman');
text( 'string',"\itNSE\rm = "+round( NSE( in, 1), 2 ), 'Units','normalized','position',[0.03,0.86], 'FontSize',12,'FontWeight','normal','FontName','Times New Roman');
text( 'string',"\itP_b_i_a_s\rm(%) = "+round( Pbias( in, 1), 2 ), 'Units','normalized','position',[0.03,0.80], 'FontSize',12,'FontWeight','normal','FontName','Times New Roman');
hl = legend(h([1 2]),'率定期','检验期'); %设置图例
set(hl,'Box','off','NumColumns',3,'Location','none','position',[0.6,0.923, 0.1, 0.05],'FontSize',14,'Fontname', '宋体');
text( 'string',"\fontname{Times New Roman}("+Title(in)+") \fontname{宋体}"+Basins (in), 'Units','normalized','position',[0.02,1.03], 'FontSize',14,'FontWeight','Bold','FontName','Times New Roman');
set(gca,'Layer','top');
str= strcat(pathFigure, " ("+Title(in)+") "+Basins(in)+"散点图" , '.tiff');
print(gcf, '-dtiff', '-r600', str);
end
%% 调用函数
function plotQQP( in, Qm , Qo , P , yearStart, yearVali , yearEnd, Pbias , NSE , R2)
% 绘制月尺度降水-实测径流/模拟图
% 输入数据
% Qmodel 模拟径流（线表示）
% Qobserve 实测径流（点表示）
% Precipitation 降水
global Title
global Basins
figureUnits = 'centimeters';
figureWidth = 24;
figureHeight = 12;
if length(Qm)==length(P)
nData = length(Qm);
else
error("降水和径流数据长度不等！");
end
nCali = yearVali-yearStart;
nVali = yearEnd-yearVali+1;
% 根据实测径流确定Y轴（流量）范围及间距
s = ceil( log10( 1.2*max(Qo) ) );
switch s
case 1
% 实测径流0-10
internalQ = 1;
case 2
% 实测径流10-100
if 1.2*max(Qo)<50
internalQ = 5;
else
internalQ = 10;
end
case 3
% 实测径流100-1000
if 1.2*max(Qo)<200
internalQ = 40;
elseif 1.2*max(Qo)<500
internalQ = 80;
else
internalQ = 150;
end
case 4
% 实测径流1000-10000
if 1.2*max(Qo)<2000
internalQ = 300;
else
internalQ = 500;
end
end
% 根据实测降水确定Y轴（降水）范围及间距
if max(P)<150
internalP = 30;
elseif max(P)<200
internalP = 50;
else
internalP = 100;
end
figure(in)
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight]);
hold on;box on;grid off;
[AX, h(2) ,h(3) ] = plotyy(1:nData, Qm,1:nData, P,'plot','bar'); % 画双轴，AX(1)左轴，AX（2）右轴，H为曲线本身
set(AX,'xlim',[0 nData+1]);
set(AX(2),'YDir','reverse','Ylim',[0, max(P)*5],'YTick',[0: internalP:max(P)*1.5],'FontSize',12,'Fontname', 'Times New Roman'); % 设置右边轴为倒立
set(gca,'box','off','Ytick',[])
set(AX(1),'YLim',[0,1.4*max(Qo)],'YTick',[0: internalQ :1.4*max(Qo)],'Fontsize',12,'YColor','k');
h(1) = plot(1:nData,Qo,'r.','LineWidth',1);
h(4) = plot([12*nCali+1 12*nCali+1],[0 1.5*max(Qo)],'k--','LineWidth', 1.2);
%设置figuer中线，柱的属性
set(h(3), 'BarWidth', 0.8 ,'facecolor',[0.2 0.4 1]);
set(h(2), 'LineWidth', 1 ,'Color', [0.28235 0.87843 0.8]); % 设置H1的曲线宽为4，颜色为g绿色,
set(h(1), 'Markersize', 8 ,'Color', [1 0.25 0]);
set(gca,'xtick',[ 1:12*5:nData ],'xticklabel', [1976:5:2020],'XMinorTick','on');
ax = gca;
ax.XAxis.MinorTickValues = 1:1*12:nData;
set(gcf,'color','w'); %设置当前figure的背景颜色
set(gca,'FontSize',14,'Fontname', 'Times New Roman');
%设置坐标轴的标题
set(get(AX(1),'Xlabel'),'String','年份','FontSize',14,'Fontname', '宋体');
set(get(AX(1),'Ylabel'),'String','\fontname{宋体}流量\fontname{Times New Roman}Q/m^