【ArcSWAT】确定水文站点控制面积及面降水量并绘制降水-径流过程图（含MATLAB全代码）

WW、forever

已于 2023-12-21 09:15:41 修改

阅读量1.2k

点赞数 2

分类专栏： # 水文模型文章标签： ArcGIS SWAT模型

于 2023-09-18 16:26:29 首次发布

本文链接：https://blog.csdn.net/qq_44246618/article/details/132983129

版权

水文模型专栏收录该内容

35 篇文章 174 订阅

订阅专栏

ArcSWAT径流模拟结果：绘制降水-径流过程图

参考

下面以洮河流域（甘肃省） 为例，水文站及气象站如下所示：
在这里插入图片描述
水文站信息如下：

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^3/s','FontSize',14);
set(get(AX(2),'Ylabel'),'string','\fontname{宋体}降雨量\fontname{Times New Roman}P/mm','FontSize',14);

% 率定期结果
text( nCali*12/2.8 , max(Qo),'率定期','FontSize',14);
text( nCali*12/2.1, max(Qo)*1.08,"\itR\rm^2 = "+round( R2( in, 2), 2 ) ,  'FontSize',12,'FontWeight','normal','FontName','Times New Roman');   
text( nCali*12/2.1, max(Qo),"\itNSE\rm = "+round( NSE( in, 2), 2 ),  'FontSize',13,'FontWeight','normal','FontName','Times New Roman');   
text( nCali*12/2.1, max(Qo)*0.92,"\itP_b_i_a_s\rm(%) = "+round( Pbias( in, 2), 2 ),  'FontSize',13,'FontWeight','normal','FontName','Times New Roman');   

% 检验期结果
text( nCali*12+nVali*12/9 , max(Qo),'检验期','FontSize',14);
text( nCali*12+nVali*12/2.4 , max(Qo)*1.08, "\itR\rm^2 = "+round( R2( in, 3), 2 ) ,  'FontSize',12,'FontWeight','normal','FontName','Times New Roman');   
text( nCali*12+nVali*12/2.4 , max(Qo),"\itNSE\rm = "+round( NSE( in, 3), 2 ),  'FontSize',13,'FontWeight','normal','FontName','Times New Roman');   
text( nCali*12+nVali*12/2.4 , max(Qo)*0.92,"\itP_b_i_a_s\rm(%) = "+round( Pbias( in, 3), 2 ),  'FontSize',13,'FontWeight','normal','FontName','Times New Roman');   

hl = legend(h([1 2 3]),'实测径流','模拟径流','降水量');           %设置图例
set(hl,'Box','off','NumColumns',3,'Location','none','position',[0.6,0.92, 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');

end


function [PCC , Pbias , NSE , KGE , R2 , RSR] = GetModelEfficience(Qs, Qo)
% 此函数用于计算水文模型评价指标

% 输入变量
% Qs    模拟径流
% Qo   实测径流

% 输出变量
% PCC       皮尔逊相关系数（Pearson’s correlation coefficient, PCC）
% Pbias     百分比偏差/相对误差（Percent bias, Pbias）
% NSE       纳什效率系数（Nash-Sutcliffe efficiency coefficient, NSE）
% KGE       克林-古普塔效率系数（Kling-Gupta efficiency coefficient, KGE）
% R2         决定/相关系数R2（Coefficient of Determination）
% RSR       均方根误差与观测值标准差的比率（RSR）
 
if length(Qs)==length(Qo)
    N = length(Qs);
    PCC = GetPCC(Qs, Qo);
    Pbias = GetPbias(Qs, Qo);
    NSE = GetNSE(Qs, Qo);
    KGE = GetKGE(Qs, Qo);
    R2 = GetR2(Qs, Qo);
    RSR = GetRSR(Qs, Qo);
else
    error("实测径流和模拟径流长度不等");
end

end


function PCC = GetPCC(Qs, Qo)
% 输入变量
% Qs    模拟径流
% Qo   实测径流
% 输出变量
% PCC       皮尔逊相关系数（Pearson’s correlation coefficient, PCC）

if length(Qs)==length(Qo)
    COV = cov(Qs, Qo);
    PCC = COV(1,2)/std(Qs)/std(Qo);
else
    error("实测径流和模拟径流长度不等");
end

end


function Pbias = GetPbias(Qs, Qo)
% 输入变量
% Qs    模拟径流
% Qo   实测径流
% 输出变量
% Pbias     百分比偏差/相对误差（Percent bias, Pbias）

if length(Qs)==length(Qo)
    Pbias = sum( (Qs-Qo)./Qo);
else
    error("实测径流和模拟径流长度不等");
end

end


function NSE = GetNSE(Qs, Qo)
% 输入变量
% Qs    模拟径流
% Qo   实测径流
% 输出变量
% NSE       纳什效率系数（Nash-Sutcliffe efficiency coefficient, NSE）

if length(Qs)==length(Qo)
    QoAve = mean(Qo);
    NSE = 1- sum( (Qs-Qo).^2./(Qo-QoAve).^2 );
else
    error("实测径流和模拟径流长度不等");
end

end


function KGE = GetKGE(Qs, Qo)
% 输入变量
% Qs    模拟径流
% Qo   实测径流
% 输出变量
% KGE       克林-古普塔效率系数（Kling-Gupta efficiency coefficient, KGE）

if length(Qs)==length(Qo)
    QsAve = mean(Qs);
    QoAve = mean(Qo);
    COV = cov(Qs, Qo);
    CC = COV(1,2)/std(Qs)/std(Qo);
    BR = QsAve/QoAve;
    RV = ( std(Qs)/QsAve )/ ( std(Qo)/QoAve );
    KGE = 1- ( (CC)^2+(BR+1)^2+(RV-1)^2 )^0.5;
else
    error("实测径流和模拟径流长度不等");
end

end


function R2 = GetR2(Qs, Qo)
% 输入变量
% Qs    模拟径流
% Qo   实测径流
% 输出变量
% R2         决定/相关系数R2（Coefficient of Determination）

if length(Qs)==length(Qo)
    QsAve = mean(Qs);
    QoAve = mean(Qo);
    R2 = ( sum( (Qo-QoAve).*(Qs-QsAve) )/ ( sum(Qo-QoAve).^2)^0.5/ ( sum(Qs-QsAve.^2)^0.5 )^2 );
else
    error("实测径流和模拟径流长度不等");
end

end


function RSR = GetRSR(Qs, Qo)
% 输入变量
% Qs    模拟径流
% Qo   实测径流
% 输出变量
% RSR       均方根误差与观测值标准差的比率（RSR）

if length(Qs)==length(Qo)
    QoAve = mean(Qo);
    RSR = ( sum(Qo-Qs).^2 )^0.5/( sum(Qo-QoAve).^2 )^0.5;
else
    error("实测径流和模拟径流长度不等");
end

end