VIC水文模型入门攻略（中）VIC的数据和输入文件准备

更正：上篇末尾关于运行汇流程序的部分，在“STEHE.stnloc”文件中需要修改的部分应改为NONE而非none。希望没有太多人因此造成困扰…

---

VIC的输入文件都有啥

其实VIC及其汇流程序真正的交互方式都是通过读入文件完成的。它们都要求有一个文件存储各种设定参数，包括运行的日期时间、运行模式以及所需的各种数据文件的位置。
上一篇中启动VIC的命令./vicNl -g /home/.../global_param.STEHE 实际上就是将VIC的全局参数文件的路径传递给VIC，然后VIC读入这些全局参数并执行。
VIC所需的必备输入文件有：

• 全局参数文件 也就是stehekin示例中的 global_param.STEHE 。

• 土壤参数文件 最核心的参数文件，保存每个网格的一些基本参数和各种土壤参数，包括网格编号，经纬坐标等，具体参考示例中的soil_param.STEHE。VIC就是从这里获得最基本的地理信息的。

• 植被参数文件 存储每个网格的土地覆盖情况，示例如veg_param.STEHE。

• 植被库文件 存储每种土地覆盖类型的一些参数，比如LAI，结构阻抗，零平面位移等。示例中为veglib.LDAS拿不到相关数据的情况下可以直接拿这个文件用。（难道就是从LDAS那里拿的？）

• 气象驱动数据文件 即meteorological forcing files（许多资料翻译作气象强迫、气象迫使，读之蛋疼，故我中二翻译为此）。存储每个网格的历时气象数据，一个网格一个文件，格式为forcing_纬度_经度。
  与上边全是文本文件形式的文件不同，这个可以为文本文件和二进制文件两种不同的类型。

其他一些可选的数据输入文件有高程带文件、雪线文件、湖泊/湿地文件以及一些基本数据的时间序列文件等。具体看官网相关内容。

---

VIC的数据准备

要生成以上相应的文件，你需要先弄到以下的数据：

• 研究区域内的多个气象站的观测数据。
  包括每个气象站的坐标以及各自的历时观测数据，精确到每日，至少包含降雨量，日平均气温或者日最高/最低气温以及风速。通常这些数据记录都不是完全的，有的数据项因为各种原因会在某些日子缺数据，甚至干脆却掉某个连续时间段的记录。于是就有必要用上数据处理工具。

• 率定用的数据。
  一些已经实际测得的VIC可输出的数据，比如某地的土壤含水量数据或者流域出口的流量数据。

• 研究区域的地理信息数据。
  首先研究区域范围的经纬度坐标表示以及DEM数据是必须的，最好还要有流域边界的shape文件，并且这些数据都能在某平台中打开和处理，比如ArcGIS。若实在没有流域边界的图形可以干脆直接开辟一块矩形区域来研究或者用高精度DEM提取盆域好了。
  DEM可以用SRTM计划的，对于VIC的尺度90m的分辨率应该够用了。链接在这里。
  有必要提醒下，VIC和汇流程序要求的坐标都是经纬度坐标，位置、网格尺寸什么的都要用经纬度表示和计算，因此在ArcGIS里面处理数据的时候数据框得设置为地理坐标系而非投影坐标系，获得的地理数据也要进行相应变换。

• 土壤数据。
  来源有很多，若手头没有的话可以使用HWSD，1km分辨率的土壤类型数据库，有全球和全国两个范围，可以从IIASA网站或者寒区旱区国家中心网站下载，优快云站内也有资源。
  不过直接从其得到的资料只有土壤类型的分类，砾、砂和粘粒的比例和一些理化性质等，所以还需要通过土壤转换函数换算为VIC需要的水力参数。若手头没有资源可以参照SPAW这个软件压缩包（本站也有资源）内附带的那篇论文里边给出的公式，或者在已发的VIC相关论文里找（记得有些论文带了土壤类型对应水力参数表）

• 土地覆盖数据（植被数据）。
  手头没数据可以用马里兰大学发布的全球1km分辨率土壤覆盖类型分布图，这是链接。点里边的Product Search -> AVHRR, Global Land Cover Product便是。

• 各种土地覆盖类型的参数数据。也就是各种土地覆盖类型的一些动力学参数以及不同月份LAI之类。其实这个通常不用特意准备，可以直接使用stehekin这个示例自带的植被库文件veglib.LDAS，其包含了常绿针叶林、常绿阔叶林等11种土地覆盖类型参数，并且与马里兰大学的土壤覆盖分布图的类型对应。
  （马里兰的覆盖分类一共有14种，其多出来的三种分别是水面、裸地和建筑用地，裸地在VIC中没有覆盖的地方默认为裸地，而建筑用地则比例太少可以忽略）

接着再来简要说说如何用ArcGIS对研究区域进行相应的网格划分。
（以下以ArcGIS10为示例，就当作你会操作ArcGIS）首先打开ArcMap，加入流域边界图形文件，打开数据框属性将坐标系设置为地理坐标系（比如1980西安坐标系），坐标表示为十进制度。流域边界也要用数据管理工具里的投影投影到相应坐标系。

接下来如果你刚用的是dem的话就可以设定经纬度范围之后重采样到你需要的网格大小。图形文件的话，为了确保之后生成的网格能尽量覆盖流域边界，使用制图工具里的掩膜工具里的要素轮廓线掩膜将边界图形扩大一点，边距通常为需生成的网格边长的0.3~0.5倍，得到一个新的稍微变胖了的图形要素。

接下来就可以使用转换工具里的转为栅格里的面转栅格将变胖的要素类转为你所需尺寸的尽可能覆盖流域边界网格。这样获得的DEM或者流向什么的都可以捕捉这个栅格进行处理了。

接着再用转换工具里的栅格转点将栅格转换成一个网格一个点的点要素，这样栅格信息的属性表就得到了。
打开点要素的属性表，添加记录坐标的两个字段“long”和“la”（其实叫什么随便了），右键点字段名选计算几何，分别得到点的在相应坐标系下的x和y坐标（单位十进制度）。

这样就可以获得一个属性表，包含网格的编号和其中点的经纬坐标。接下来就可以将其导出成数据库或者文本文件以便用R或者matlab什么的处理了。

---

各种输入文件的准备

关于全局参数文件（global parameter file）

全局参数文件可以直接copy示例stehekin的或者拿源码压缩包里带的samples改成自己需要的直接用。（当然文件名和文件路径就别照搬上边的XXX.STEHE了，那个STEHE指的就是stehekin项目，完全可以自己安排一个简洁清晰高效的文件名和组织形式…）

用文本编辑器（对于一些还在用UE甚至记事本的水文水利的同学我推荐下sublime text 3这个更轻更高效的文本编辑器）打开示例的全局参数文件之后可能刚开始会被大量的内容搞晕，但是其实大部分都是#号开头的注释，内容为一些说明文字以及一些可选的参数设置。
可以通过设置代码高亮（比如Perl、R等用#注释的语言）使之看起来清晰。

官网对各种参数的解释戳这。下面是一些下垫面相关基本参数的自翻：

• 基本参数部分

NLAYER      3   # 土壤层数，决定土壤参数文件中许多参数要弄几份
NODES       10  # 土壤热结点个数 
TIME_STEP   24  # 每个时间步的长度，单位小时。如果设为24那就是逐日。
SNOW_STEP   1   # 雪模型的时间步长度。如果上一项小于24就要等于上一项。
STARTYEAR   1949    # 模型运行的开始年份
STARTMONTH  01  # 模型计算时段的开始月份
STARTDAY    01  # 开始日期
STARTHOUR   00  # 从几点钟开始...
ENDYEAR     1951    # 结束年份
ENDMONTH    12  # 结束月份
ENDDAY      31  # 结束日期

FULL_ENERGY     FALSE   # 是否开启能量平衡模式，如果为FALSE就只开启水量平衡模式。
FROZEN_SOIL FALSE   # 是否开启冻土计算模式。

SOIL           /home/heheda/xriver/soil.duang   # 土壤参数的路径（这路径是我瞎编的）
ARC_SOIL        FALSE   # 土壤数据是不是ArcGIS的相关格式。通常用文本文件所以是FALSE。
BASEFLOW    ARNO    # ARNO表示土壤参数文件中第5到8列为VIC标准基流参数，为 NIJSSEN2001表示此处为Nijssen的基流参数（一般默认）
JULY_TAVG_SUPPLIED  FALSE   # 是否有在土壤参数最后添加一列六月份平均气温。如果没有程序就自己计算。
ORGANIC_FRACT   FALSE   # 是否进行土壤有机质相关的计算。如果为TRUE那么土壤参数就得额外再加三种参数：每层土壤里边有机质占总的比例，每层土壤的有机质的容积密度以及土粒密度。
VEGLIB          /home/heheda/xriver/veglib.LDAS # 植被类型库文件路径
VEGPARAM        /home/heheda/xriver/veg.duang   # 土地覆盖（也就是植被覆盖）文件路径（也是我瞎编的）
ROOT_ZONES      3   # 植被的根区数目（这个决定植被覆盖文件里边需要设定的相应参数的个数）
VEGPARAM_LAI    TRUE    # 植被覆盖文件是否自带LAI信息
LAI_SRC     LAI_FROM_VEGLIB    # 参数为LAI_FROM_VEGPARAM表示从植被参数文件里边读取LAI，参数为LAI_FROM_VEGLIB则表示从植被库文件里读取LAI。
SNOW_BAND   1   # snow bands（雪线？）的个数。求大牛解释其意思...我只知道如果大于1就要准备相应的文件了。

对于气象驱动文件输入设定和数据输出部分详见下面相关内容。

关于气象驱动文件（Meteorological Forcing FIle）

气象驱动文件的数据类别比较自由，可以在VIC提供的一个列表里选取可以提供的数据类别。
写此文时此链接markdown表格格式挂了故帖上一些过来：

气象
数据类别	是啥	描述
AIR_TEMP	平均气温	C (ALMA单位: K)
CATM	二氧化碳浓度	ppm
PREC	降水量（包括降雨、降雪）	mm (ALMA单位: mm/s)
DENSITY	空气密度	kg/m3
LONGWAVE	入射长波辐射量	W/m2
PREC	总降水量（包括降雨和降雪）	mm (ALMA单位: mm/s)
QAIR	比湿度	kg/kg
REL_HUMID	相对湿度	分数
SHORTWAVE	入射短波辐射量	W/m2
SNOWF	降雪量	mm (ALMA单位: mm/s)
TMAX	时段内最高气温	C (ALMA单位: K)
TMIN	时段内最低气温	C (ALMA单位: K)
TSKC	云量	分数
VP	水汽压	kPa (ALMA单位: Pa)
WIND	风速	m/s
WIND_E	风速在正东方向的分量	m/s
WIND_N	风速在正北方向的分量	m/s
植被
ALBEDO	地面反照率	分数
LAI_IN	叶面积指数	分数
VEGCOVER	植被覆盖率	分数
其他
SKIP	留出位置待将来放其他数据	没

其中降雨量PREC、平均气温AIR_TEMP或者日最高、最低气温TMAX和TMIN以及风速WIND三样为必选项。具体选啥还主要看手头上的气象数据有啥及其质量如何。

气象驱动数据输入文件形式为一个网格一个文件，全部放在同一个目录下，文件名格式示例data_12.3456_123.4560，前缀data_可以更改，但是经纬度后缀格式应按照示例，先纬度后经度，小数点后保留某位小数。
气象驱动数据的输入文件支持两种格式，文本文件（ASCII）格式和二进制文件（binary）格式，提供后者是为了优化数据的读取性能。

文本文件格式的气象驱动数据

如果研究规模比较小，对性能没有太大纠结的时候通常用比较方便的文本文件格式。其在全局参数文件中相应位置设置如下：

FORCING1    /home/hehehe/yriver/data_   # 气象驱动数据的路径格式。
N_TYPES     4   # 包含的气象数据类别数。
FORCE_TYPE  PREC    # 气象数据类别：降水量
FORCE_TYPE  TMAX    # 气象数据类别：最高气温
FORCE_TYPE  TMIN    # 气象数据类别：最低气温
FORCE_TYPE  WIND    # 气象数据类别：风速
FORCE_FORMAT    ASCII   # 输入文件格式：文本文件
FORCE_DT    24  # 数据时间步长度
FORCEYEAR   1950    # 数据记录的开始年份
FORCEMONTH  1   # 数据记录的开始月份
FORCEDAY    1   # 数据记录的开始日期
FORCEHOUR   0   # 数据记录的开始时间
GRID_DECIMAL    4   # 输入文件的经纬度小数点后保留几位小数。
WIND_H          10.0    # 测得风速的高度（m）。
MEASURE_H       2.0     # 测得湿度的高度（m）。
ALMA_INPUT  FALSE   # 数据是否使用ALMA单位制。

输入文件的格式也很简单，每一列是一种数据类别，从左到右按全局参数文件里设定的顺序；每一行则是逐时间步的数据，形式如：

0.000 -10.00 -2.00 2.00
0.000 -8.50 -1.10 2.00
5.000 -7.90 0.10 2.00
10.000 -8.20 -3.10 2.00