【SWAT水文模型】SWAT水文模型建立及应用第四期：气象数据的准备（传统气象站）（含MATLAB数据整理代码）

WW、forever

已于 2023-08-01 22:44:47 修改

阅读量5.6k

点赞数 14

分类专栏： # 水文模型文章标签： ArcSWAT 气象数据

于 2023-04-26 21:51:45 首次发布

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

版权

水文模型专栏收录该内容

35 篇文章 181 订阅

订阅专栏

SWAT水文模型建立及应用：气象数据的准备

1 简介
- 1.1 数据来源
2 气象数据的准备（传统气象站）
3 SWAT模型中气象数据的导入
参考

本博客主要介绍气象数据的准备，分为数据下载和数据处理两部分。

1 简介

SWAT模型需要的气象数据主要包括流域的日降水量、最高/最低气温、太阳辐射、风速和相对湿度。这些数据可以是统计数据，也可以通过SWAT模型的天气模拟程序（Weather Generator）生成，或者是统计和模拟生产数据的结合。
理想情况下SWAT模型需要至少20年的天气资料，如果缺少必要的资料则可用天气发生器进行补充。
在时间尺度上，模型的模拟时间步长可以为年、月、日。

1.1 数据来源

降水、日最高／最低气温以及风速、相对湿度、太阳辐射等气象数据由中国气象数据共享服务网提供。

此案例研究区域为洮河流域（甘肃省），其周边气象站点分布如下：
在这里插入图片描述
数据系列长度为整个研究期（包括率定期和检验期）：1970-2020年。站点信息如下：

2 气象数据的准备（传统气象站）

气象数据对水文过程的重要性是不言而喻的。在SWAT 模型建立过程中有三个数据是模型所必须得，即天气发生器、降水数据、气温数据、相对湿度、风速和日照（可选），前者因其可以弥补气象数据的缺失，是SWAT 模型内置的，必须在建模之前提前建立好数据库信息，后两者可以从气象站点获取数据。

关键步骤：

天气发生器各参数的计算
降水及气温输入数据的准备

2.1 天气发生器各参数的计算

天气发生器可以根据多年逐月气象资料模拟生产逐日气象资料，但该数据库要求输入的参数较多，其主要输入数据有月平均最高气温、月平均最低气温、最高气温标准偏差、月平均降雨量、降雨量标准偏差、月内干日日数、露点温度、月平均太阳辐射量等。
天气发生器参数的计算公式如下：
在这里插入图片描述
天气发生器的详细介绍可参见另一博客-【SWAT水文模型】SwatWeather软件使用教程。

打开SWAT2012.mdb数据库，WGEN_user表，如下：
在这里插入图片描述
将其导出为excel表，并将研究区域内的站点资料整理进去。

随后，将整理好的表WGEN_user导入到SWAT2012.mdb数据库中的WGEN_user表中。在导入表时，选择【我自己选择主键】，如下图：

2.2 降水及气温输入数据的准备

在模型的气象输入文件中，降水与气温数据是必须的文件。降水所需尽可能的长时间序列实测数据，其格式如下图所示分为两列，时间序列必须为连续的，没有测值用-99代替，存储为.DBF格式。因为降水直接对径流有着重要的影响，虽然天气发生器可以模拟日降水数据，但强烈建议模型使用者直接输入实测的日尺度降水数据。
如果降水、气温的站点与天气发生器站点为同一个，可以由SwatWeather.exe直接生成所需格式的输入文件，如果有更多的站点，需要自行转化成相应格式的输入文件，建议降水与气温数据采用尽可能多的站点。

2.2.1 降水数据准备

SWAT模型气象数据的输入格式可以在软件安装后自带的案例中查看。
以降水为例，每一个站点的降水数据要形成一个单独的文本文件，文件命名为站点的名称，降水的索引文件和站点数据必须要在同一个文件夹下面，在SWAT输入中只需要输入索引文件，它会在同一个文件夹下面自动寻找站点数据。
在这里插入图片描述
具体的格式：每一种类型的数据需要有一个站点的索引文件，比如降水，需要有一个降水站点的索引文件，索引文件中需要包含站点的名称（具体表头的样子参考自带的案例），索引文件具体的样子如下：一定要严格按照格式输入，后三者为纬度、经度和高程。
每一个站点的数据要注意的点：

必须从模拟的第一天到最后一天，每天都要有一个数据，即便数据缺失也要用-99代替，一定要确保数据没有缺失，每个站点的数据长度都是相同的；
每个站点名称一定要与站点索引文件中的站点名称（NAME）字段对应起来；
注意异常值处理，网上下载的数据通常都有异常值，比如9999之类的，都必须替换为-99

PCPfork文件的MATLAB实现代码如下：

clc
close all
clear
%% 基本设置
pathFiles= '.\DataFiles\' ;

%% 导入数据

load('StationYT.mat');
nStation = length( StationYT );

%% 导出数据 ID NAME LATITUDE LONGITUDE ELEVATION
nVariable = 5;
Title = 'PCPfork';
fid = fopen([pathFiles,'PCPfork.txt'],'wt');
fprintf(fid,'%s\n','ID,NAME,LAT,LONG,ELEVATION');
for iStation=1:nStation
    for iVariable=1:nVariable
        switch iVariable
            case 1
                fprintf(fid,'%s',num2str( StationYT(iStation,1) ) );
                fprintf(fid,'%c',',');
            case 2
                fprintf(fid,'%s','PCP');
                fprintf(fid,'%s',num2str( StationYT(iStation,1) ) );
                fprintf(fid,'%c',',');
            case 3
                fprintf(fid,'%0.4f',StationYT(iStation,2));
                fprintf(fid,'%c',',');
            case 4
                fprintf(fid,'%0.4f',StationYT(iStation,3));
                fprintf(fid,'%c',',');
            case 5
                fprintf(fid,'%0.2f\n',StationYT(iStation,4));
        end
    end
end

输出文件格式如下：
在这里插入图片描述
各站点日降水数据的处理如下：

clc
close all
clear
%% 基本设置
pathFigure= '.\Figures\' ;
pathFiles= '.\DataFiles\' ;

%% 站点数据资料处理
load('StationID.mat');                                % 840个站点的站点信息                      
load('StationIDSSQ.mat');                          % 研究区内站点信息（区站号 X Y 位置）
nStation = length(StationIDSSQ);

load('Prec.mat');                   % 降水
Datamiss = -99;

% 截取实际采用数据（1970-2020年 共51年数据）18628天
% ------------------------------------------------------------------
nDay = datenum(2020,12,31)-datenum(1970,1,1)+1;
ii = datenum('01-Jan-1970');
jj = datenum('31-Dec-2020');
% 数据说明：分闰年/平年
if (jj-ii+1)~=nDay
    sprintf("数据有误！")
end

%% 研究区内降水数据处理 日降水资料PRE

% 根据研究区内站点ID得到相应降水数据
% ------------------------------------------------------------------
ID = StationIDSSQ;
PrecSSQ = zeros( length(ID) , length( Prec(1,:)) );
for ilength = 1:length( ID )
    for jlength = 1:length( StationID )
        if ID(ilength)==StationID(jlength)
            PrecSSQ(ilength,:) = Prec(jlength,:);
        end
    end
end

% 降水数据处理
% 32700	表示降水"微量"
% 32XXX	XXX为纯雾露霜
% 31XXX	XXX为雨和雪的总量
% 30XXX	XXX为雪量（仅包括雨夹雪，雪暴）
% ------------------------------------------------------------------
PSSQ = zeros( nStation ,  nDay );
for iStation=1:nStation
    for iday=1:nDay
        if isnan( PrecSSQ( iStation,iday) )
            PSSQ( iStation,iday) = Datamiss;
        elseif PrecSSQ( iStation,iday)==999990
            PSSQ( iStation,iday) = 0;
        elseif PrecSSQ( iStation,iday)==32700
            PSSQ( iStation,iday) = 0;
        elseif PrecSSQ( iStation,iday)>32000
            PSSQ( iStation,iday) = PrecSSQ( iStation,iday)-32000;
        elseif PrecSSQ( iStation,iday)>31000
            PSSQ( iStation,iday) = PrecSSQ( iStation,iday)-31000;
        elseif PrecSSQ( iStation,iday)>30000
            PSSQ( iStation,iday) = PrecSSQ( iStation,iday)-30000;
        else
            PSSQ( iStation,iday) = PrecSSQ( iStation,iday);
        end
    end
end

% txt格式数据导出   起始数据年月日  各日降水量（mm）
% ------------------------------------------------------------------

for iStation=1:nStation
    Title = ['PCP',num2str( StationIDSSQ(iStation,1) )];
    fid = fopen([pathFiles,'PCP\',Title,'.txt'],'wt');
    fprintf(fid,'%s\n','19700101');
    for iday=1:nDay
        fprintf(fid,'%0.2f\n', PSSQ( iStation,iday) );
    end
end

输出文件夹如下所示：
在这里插入图片描述

2.2.2 气温数据准备

气温数据输入格式如下：

TMPfork.txt文件储存各站点信息
各子文件夹存储气温数据

MATLAB处理代码与降水处理代码类似，如下：

clc
close all
clear
%% 基本设置
pathFiles= '.\DataFiles\' ;

%% 导入数据

load('StationYT.mat');
nStation = length( StationYT );


%% 导出数据 ID NAME LATITUDE LONGITUDE ELEVATION
nVariable = 5;
Title = 'TEMPfork';
fid = fopen([pathFiles,'TEMPfork.txt'],'wt');
fprintf(fid,'%s\n','ID,NAME,LAT,LONG,ELEVATION');
for iStation=1:nStation
    for iVariable=1:nVariable
        switch iVariable
            case 1
                fprintf(fid,'%s',num2str( StationYT(iStation,1) ) );
                fprintf(fid,'%c',',');
            case 2
                fprintf(fid,'%s','TEMP');
                fprintf(fid,'%s',num2str( StationYT(iStation,1) ) );
                fprintf(fid,'%c',',');
            case 3
                fprintf(fid,'%0.4f',StationYT(iStation,2));
                fprintf(fid,'%c',',');
            case 4
                fprintf(fid,'%0.4f',StationYT(iStation,3));
                fprintf(fid,'%c',',');
            case 5
                fprintf(fid,'%0.2f\n',StationYT(iStation,4));
        end
    end
end