圣维南方程python代码 圣维南方程组

汇流程序的数据准备

前面说过VIC作为分布式水文模型，其功能只有计算产流量而无法给出汇流量，于是要进行汇流操作得用Dag Lohmann研发的一个汇流模型程序进行，其名为Routing，编程语言为Fortran，使用单位线进行坡面汇流以及线性圣维南方程进行河网汇流，输入数据为VIC的输出文件，流域的水系流向信息，以及你想得到模拟的径流量数据的若干个流量站或者流域出口的位置，输出为逐日和逐月的此站点模拟的径流量数据。

汇流程序运行的原理

汇流程序用法跟VIC一样，也是需要准备一些文本文件格式的数据文件和一个输入参数文件来记录这些数据文件的位置和一些参数，然后也是先make编译，然后输入./rout后面带上输入参数文件的路径来运行。具体必须的输入文件有：

• 输入参数文件。相当于VIC的全局参数文件。
• 流向文件。这是一个栅格文件，也是最核心最重要的文件，记录流域内每个网格所流向的下一个网格。其实不单单是流向，汇流程序需要的流域的网格大小、经纬坐标、流域形状等信息都是从这儿获得的。
• 站点文件。记录流量站的名称和位置。
• 坡面汇流单位线文件。这个可以直接拿示例stehekin的来用。

然后这些是可选文件：

• 流速文件。也是栅格文件。不过流速之类的参数都不好弄到，所以可以只用一个单一值代替。这个值和下面的水力扩散值以及坡面汇流单位线是可以率定的。不过根据官网文档关于率定的部分所说，对于大尺度、逐月的汇流计算，由于要求精度较低一般不需要率定，根据前人的研究此时可以采用一个单一值1.50。
• 水力扩散文件。同上，可用单一值800。
• Xmask文件。栅格文件。其实不太清楚为啥叫这个名字，官网文档的一行小字解释这其实就是流程数据，记录的是每个网格中心到其流向的下一个网格中心的距离，对精度要求不高设定为一个单一值也可以。我将其解释为流程文件。
• fraction文件。记录的是每个网格实际产流量的比例，详细解释看下文。我将其解释为产流比例文件。
• 站点单位线文件。后缀名.uh_s，一个站点一个，记录的是流域内每个网格对流域出口也就是站点位置的单位线。刚开始用汇流程序的时候不需要这个文件，因为汇流程序可以生成，详见下文。

程序运行的流程是这样的：

1. 汇流程序首先会读入输入数据文件获取各项参数（废话）；
2. 然后读入各种栅格文件的数据和坡面汇流单位线数据；
3. 然后打开站点文件，读入第一个站点的信息，比如站点名、位置；
4. 根据流向文件的数据找出最终会流入这个站点的所有网格（即站点控制的流域内的所有网格）并记录到列表；
5. 如果该站点信息的站点单位线文件路径不为“NONE”则到第6步，否则到第7步；
6. 根据此路径读取站点单位线数据，然后到第8步；
7. 利用流向、流速、水力扩散和流程数据和坡面汇流单位线数据使用Dag Lohmann的方程计算出站点单位线数据，并保存到文件以供下次直接使用；
8. 从第4步得出的网格列表打开第一个网格对应的VIC输出文件，读入产流数据；
9. 根据产流比例数据和站点单位线数据计算出此站点对流域出口的径流量并加到站点的总径流量上；
10. 如果还有没计算完的网格则回到第8步拿下一个网格继续；
11. 输出该站点的径流量序列（包括逐日和逐月的径流量（立方英尺每秒）和径流深（mm）表示的数据）；
12. 如果站点文件还有下一个站点则读入下一个站点的信息然后回到第4步；
13. 收工。

说白了就是假设每个网格的产流对流域出口径流的贡献互相独立，先是生成一套网格产流对出口径流的单位线，然后拿这套单位线和网格产流量做汇流计算，而且这两个步骤是可以分开的。

关于输入参数文件

内容比较简单，不过需要注意的是汇流程序的健壮性实在不咋地，最好是复制示例的输入参数文件rout_input.STEHE对原文进行修改，信息的顺序也不要改，甚至注释也不要增删。因为我看到源码就是简单粗暴的跳过一行读取一行来读取信息的… 内容自翻如下：

# 这行不要删。   
# 流向文件的路径。   
/mnt/aww/ohio/route/rout_data/oh.flow   
# 流速文件。如果有就改为.true.然后下一行改为文件路径，没有就如下设定一个单一值。
.false.   
1.5   
# 水力扩散文件，同上。   
.false.   
800   
# 流程文件Xmask file，同上。 
.true.   
/mnt/aww/ohio/route/gen_mask/kanawha.xmask
# 产流比例文件fraction file，同上。   
.true.   
/mnt/aww/ohio/route/rout_data/oh.fract
# 站点文件。   
/mnt/aww/ohio/route/run/oh.stnloc.kanaw
# VIC输出文件的路径和格式，以及文件名坐标保留的小数位数。这个应跟VIC的设定对应。   
/mnt/aww/ohio/vic_ohio/results/fluxes_   
4   
# 输出文件的路径。
/mnt/aww/ohio/route/rout_out/kanaw/
# 第一行是进行汇流计算的时间段，第二行是实际输出到文件里的数据的时间段。
1990 1 1996 12
1990 1 1996 12
# 坡面汇流单位线文件。可以用stehekin提供的。
/mnt/aww/ohio/route/rout_data/oh.uh_all

关于流向文件

把官网的样例搬过来：

ncols           22
nrows           20
xllcorner       -97.000
yllcorner       38.000
cellsize        0.50
NODATA_value    0
0 0 0 5 5 5 6 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 4 0 4 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 5 3 5 5 6 7 5 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 5 4 0 5 5 5 5 6 7 5 5 5 7 5 0 0 0 0 0 0
0 0 5 4 0 5 4 5 5 7 7 5 6 7 5 6 0 0 0 0 0 0
4 0 5 7 4 0 3 5 5 7 5 5 7 4 5 6 0 0 0 0 0 0
4 0 3 4 0 4 0 4 3 5 6 7 7 4 5 6 0 0 0 0 0 0
0 4 0 4 0 3 1 2 4 0 4 5 7 5 5 0 0 0 0 0 0 0
0 0 3 4 3 1 8 4 5 4 0 5 5 3 5 6 5 0 0 0 0 0
0 0 0 3 4 5 5 4 5 4 3 5 6 7 7 5 5 7 0 0 0 0
0 0 0 4 3 5 5 3 4 0 4 3 4 3 5 5 7 7 0 0 0 0
0 0 0 4 0 5 4 0 4 3 4 0 4 5 3 5 7 7 0 0 0 0
0 0 0 0 4 0 4 0 4 0 4 5 3 6 7 7 5 5 6 0 5 7
0 0 0 0 0 4 0 3 4 3 5 3 6 7 6 5 7 7 7 7 7 7
0 0 0 0 0 0 0 0 4 0 4 5 7 5 6 7 1 7 1 7 0 0
0 0 0 0 0 0 0 0 3 4 0 5 5 6 7 7 7 6 0 0 0 0
0 0 0 0 0 0 0 0 0 4 0 4 5 7 1 7 7 0 0 0 0 0
0 0 0 0 0 0 0 0 0 3 2 4 5 7 7 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 4 0 5 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

如果你玩过ArcGIS的从栅格导出到ASCII工具你应该会发现这其实就是标准的导出成ASCII文件的ArcGIS的栅格数据的格式。其中第一二行记录的是栅格的列数行数，第三四行记录的是栅格左下角（不是最左下角的网格的中心而是最左下角的网格的左下角）的坐标，第五行定义的是栅格的尺寸，第六行是栅格中无数据栅格的表示形式。 然后下面一行行一列列的就是每个网格的值。 汇流程序使用D8流向定义方式，用整数1~8定义8个不同的方向，表示各网格周围的8个水流可以流向的下一个网格。 整数1到8对应的方向分别为北、东北、东、东南、南、西南、西和西北，即说八个方向从北方向开始顺时针用整数1到8编码。 （这里跟ArcGIS的编码方式不同，ArcGIS用的是整数2n，其中n为整数0到7，编码方向也是从东方向开始顺时针，也就是说1、2、4、8、16、32、64、128分别对应东、东南、南、西南、西、西北、北和东北。）

至于流向文件的准备，不要直接拿ArcGIS的流向工具生成。不要直接拿ArcGIS的流向工具生成。不要直接拿ArcGIS的流向工具生成。重要的事情说三遍。 因为ArcGIS的这个工具用的是以相同分辨率的DEM为依据，假定每个网格流向其周边最低的网格这种简单粗暴的方式生成流向的，VIC需要的这么大的网格尺度，经过重采样之后的DEM肯定会抹掉网格内地形的细节，到时候河网到底怎么流都不知道了，弄出来的流向肯定是错的。 所以最好的方式还是用高精度的DEM提取水系再想办法搞。这是个比较蛋疼的问题。实在不行可以参考我的略二百五的做法：拿90m的DEM用ArcGIS的流量工具提取水系，然后记录每个网格范围内的流量的最大值，再根据水量累积的原理假设每个网格流向其邻近的流量最大的网格或者流向边界（如果边界是海岸线的话）。这种方法好处就是简单粗暴，划分出来的流向也大致是对的，缺点就是在一些地形奇葩的地区可能会形成回路，直接放进汇流程序跑会死循环，得用一个程序检测排除。

关于站点文件

站点文件就简单多了。形式如下：

1 UPMIS               14  2  -9999
NONE
1 WISCR               14 11  -9998
WISCR.uh_s
1 ANOKA                8 15  -9998
store/ANOKA.uh_s
...

每个站点占用两行。第一行从左到右依次是：是否运行这个站点，是则为1，否为0；站名；站点在流向文件定义的栅格中位于左到右位于第几列；站点位于从下到上第几行；站点控制的流域的面积（用不上所以设成一串负数）。

第二行则是站点的单位线文件的路径，也就是记录站点控制的流域内各网格对站点位置的单位线数据的文件，这个不用自己准备。因为前面说过，汇流程序假设每个网格的汇流互相独立，当设定为NONE的时候，汇流程序会根据一套算法计算出单位线，然后保存到汇流程序所在的目录中生成.uh_s文件，接下来程序才会读入VIC输出的各网格产流数据并拿这些单位线计算径流过程。

所以下次计算的时候就可以将程序生成的单位线文件直接拿来用，设置好单位线文件的路径，而不用再进行单位线计算了。 这个机制主要是由于计算单位线的时间通常要比卷积的时间长不少，而且一次算完可多次利用，以提高效率。

关于产流比例文件Fraction file

再次搬运官网内容：

ncols           22
nrows           20
xllcorner       -97.000
yllcorner       38.000
cellsize        0.50
NODATA_value    0.00
0.00 0.00 0.00 0.11 0.38 0.63 0.65 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.78 1.00 1.00 1.00 0.83 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.64 1.00 1.00 1.00 1.00 0.87 0.35 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.83 1.00 1.00 1.00 1.00 1.00 0.93 0.73 0.90 0.64 0.39 0.49 0.44 0.33 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.91 0.00 0.00 0.00 0.00 0.00 0.00
1.00 0.99 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.75 0.00 0.00 0.00 0.00 0.00 0.00
0.52 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.29 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.57 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.98 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.57 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.73 0.21 0.40 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.26 0.82 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.56 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.41 0.99 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.82 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.34 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.84 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.57 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.96 0.15 0.00 0.28 0.28
0.00 0.00 0.00 0.00 0.00 0.48 0.97 0.99 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.88 0.83 0.65 0.43
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.53 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.91 0.99 0.65 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.13 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.42 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.93 1.00 1.00 1.00 1.00 1.00 1.00 0.42 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.52 0.86 1.00 1.00 0.87 0.56 0.30 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.31 0.57 0.05 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

可能你直接一眼就看出这个文件其实是干嘛用的了。 这个栅格文件记录的是每个网格被包含到实际的流域边界内的部分的面积占网格面积的比例。

现实中的流域边界往往都要打位于流域边缘的网格中间经过，因而实际上位于流域边界的网格的产流量也不会全部流到流域出口，而是有一部分会流到隔壁流域。也就是说每个网格只有被包含在实际的流域边界的部分的产流量才会流入到流域，最终汇聚到流域出口。因而为了提高精度，汇流程序可以接受记录每个网格实际产流比例的数据，在进行汇流计算的时候将VIC输出的产流数据乘上对应网格的产流比例。

而汇流程序在这里假定VIC每个网格的产流量空间上均匀分布。所以每个网格的产流比例，实际就是网格被包含到实际流域边界内的面积比例与网格总面积的比例。

而这个文件在ArcGIS中的获得方法可将栅格转为图形要素（就是将网格变成一个个正方形图形的形式），然后和流域边界图形要素用分析工具的相交工具做交操作，计算面积后除以完整网格的面积即可得到产流比例。之后再从图形转回栅格（栅格值为产流比例值）并用转换工具中的栅格转ASCII即可得到此文件。

当然要是你的网格精度够高，可以忽略流域边界与网格不重合的误差，也可以全部设个单一值1.0。

注意事项

根据我的发现许多人的汇流程序运行的成功率要比VIC小不少。这有一个重要原因是官网给出的汇流程序的健壮性实在是不怎么样。 这主要是因为编程语言的本身性质问题。Routing这个程序由Fortran所写，这是一种非常古老的高级语言，比C和BASIC都要老，基本上再往前就是汇编了。通常这种早期的语言的内存管理方式都很原始，数据结构都要预先定义所需要的空间大小，许多简单的功能都难以实现。也或许作者已经在建模上投入了太多精力，无暇顾及实际运行中的一些细节。 至于这个程序的C版本，官方表示在搞着但是没搞好… 于是没办法，只好由人迁就模型。目前发现的一些要遵守的细节如下：

• 首先各种输入文件最好从示例照搬，别随便填删改注释。
• 如果你建立的网格规模较大，网格的行数或列数超过250，模拟的时长超过50年，或者站点控制的流域内的网格数超过10000，请将源文件rout.f中的相应限制值变量改为你需要的值（分别为位于32、33和44行的NROW、NCOL、NYR、PMAX）。
• 文件路径长度要控制，最大为72个字符，如实在不行请将源码中相应的定义字符串长度的值改为你所需要的值。
• VIC输出文件的设定最好不要更改，否则其中需要包含这两样数据：OUT_RUNOFF和OUT_BASEFLOW，即径流量和基流量，且在输出文件中应位于除日期列以外的第3和第4列，参数PRT_HEADER应设为FALSE，即不要带带注释符号“#”的表头，因为汇流程序认不出注释。
• 程序输出的站点名称为固定5个字符，多了裁掉后面的，少了用空格补充，所以不要用同样的一串长于5个字符的字符串后面加个编号比如OUTTST01，OUTTST02来作为站点名，因为在程序里一处理通通都变成OUTTS，然后输出数据的时候就互相把输出文件覆盖掉了。

顺便贴个我的针对这些问题自己用C++撸的汇流程序CRout：https://github.com/Sibada/CRout 汇流机理和公式照搬原版，只是在程序结构上有所修改。 此外增加了指定产流数据在VIC输出文件中的列数等选项。 用法：除了输入参数文件和站点单位线文件有修改其他照原版。 欢迎debug。

一些想补充的东西

MinGW的配置

其实对于MinGW和Cygwin的性能区别我是没深入研究的，主要是靠MinGW是直接编译成win32API的程序所以性能会高些。当然还是不太推荐在Windows下运行VIC，因为速度还是比不上linux的。 这里就给出下MinGW的配置方法，各位水文人员可以在数据量不大的情况下在Windows下玩玩。

因为使用官方的在线安装程序一来下载速度不稳定，二来安装选项可能会使不太懂GNU的人眼花缭乱，所以自己搞了个傻瓜包共享：http://pan.baidu.com/s/1i4jwCnz

先将里面的MinGW.zip下载下来，其中的MinGW文件夹解压到某地，通常推荐是C盘；

然后右键点我的电脑->属性，再如下操作配置环境变量：

在变量值的框里加上;C:\MinGW\bin，也就是英文分号加上你放MinGW的路径后边再加个\bin。

然后确定确定确定。再点开开始菜单输入CMD回车，输入make回车显示这个就说明能用MinGW了：

然后就可以像linux那样编译运行vic和rout了。

土地覆盖（植被覆盖）文件的一种生成方法

因为植被覆盖文件格式比较蛋疼，不像土壤文件那样是规规矩矩的二维表，于是分享下我用的比较简单粗暴的方法。 灵感来自我老师。 并欢迎其他大牛提出更有效率的方法。 本方法得出的植被覆盖文件并不包括LAI数据，所以使用时全局参数文件中VEGPARAM_LAI应设为FALSE。

接上一篇，首先你有了研究区域的栅格和栅格转成的点图层：（点忘了在图上显示出来…）

然后使用转换工具的点转栅格，其中值字段选pointid，像元大小还跟你的设定；

然后就得到了栅格值为你的网格ID的栅格。

然后载入马里兰大学的全球1km植被覆盖文件：

使用空间分析工具->地图代数->栅格计算器工具，表达式为网格ID的栅格乘100加上马大的植被覆盖栅格；

先别点确定，点取消右边的环境…，然后进行如下设置：

上面设置为网格ID的栅格，作用是框定生成的栅格的界限，

下面设置为马大的植被覆盖栅格，作用是将生成的栅格尺寸定位马大植被覆盖栅格的尺寸。然后得到这样一个东西。

然后是再从栅格转为点，用栅格转点工具，就会生成一个密密麻麻的点要素图层。

打开这个点要素图层的属性表。

其中GRID_CODE列就是我们要的数据。可发现这个数据后两位为植被类型编号，前若干位则为其所属的网格ID。这就是乘上100的原因。然后导出属性表数据，保存为文本文件，得到如下内容：

接下来就是弄一段程序来统计各网格内植被类型的比例。 比如这个我用R写的。

# 上一步导出的文本文件的路径
veg_stat = read.table('d:/grid_veg_point.txt',header = T,sep = ',')[,3]

root_pm = c('0.10 0.05 1.00 0.45 5.00 0.50',
            '0.10 0.05 1.00 0.45 5.00 0.50',
            '0.10 0.05 1.00 0.45 5.00 0.50',
            '0.10 0.05 1.00 0.45 5.00 0.50',
            '0.10 0.05 1.00 0.45 5.00 0.50',
            '0.10 0.10 1.00 0.65 1.00 0.25',
            '0.10 0.10 1.00 0.65 1.00 0.25',
            '0.10 0.10 1.00 0.65 0.50 0.25',
            '0.10 0.10 1.00 0.65 0.50 0.25',
            '0.10 0.10 1.00 0.70 0.50 0.20',
            '0.10 0.10 0.75 0.60 0.50 0.30')
            # 各种植被的根区数据，此处为ROOT_ZONES=3的情况。

id = round(veg_stat/100)
veg_type = veg_stat - id*100

veg_stat = data.frame(id,veg_type)
veg_stat = veg_stat[order(veg_stat$id),]

vn = 0
tl = nrow(veg_stat)

type_stat = rep(0,11)

out_text = c()
r = 1

aban = 0.02 # 面积比例少于百分之二的植被类型不登记

for(i in 1:tl){
  current_id = veg_stat$id[i]
  tp = veg_stat$veg_type[i]

  vn = vn + 1

  if(tp %in% c(1:11))
    type_stat[tp] = type_stat[tp] + 1

  if(i == tl || veg_stat$id[i + 1] != current_id){
    type_stat = type_stat / vn

    type_sum = length(type_stat[type_stat > aban])
    out_text[r] = paste(current_id,type_sum, sep = '\t')
    r = r + 1

    for(t in 1:11){
      if(type_stat[t] > aban){
        out_text[r] = paste(t,round(type_stat[t],5),root_pm[t], sep = '\t')
        r = r + 1
      }
    }
    type_stat = rep(0,11)
    vn = 0
  }
}

# 输出的文件
write.table(out_text,'d:/veg_param.txt', row.names = F,col.names = F,quote = F)

跑完程序就可以去找输出的植被覆盖文件了。 这个是效果：

1   5
2   0.12    0.10 0.05 1.00 0.45 5.00 0.50
6   0.27556 0.10 0.10 1.00 0.65 1.00 0.25
7   0.34222 0.10 0.10 1.00 0.65 1.00 0.25
10  0.04444 0.10 0.10 1.00 0.70 0.50 0.20
11  0.19556 0.10 0.10 0.75 0.60 0.50 0.30
2   6
2   0.04    0.10 0.05 1.00 0.45 5.00 0.50
4   0.03556 0.10 0.05 1.00 0.45 5.00 0.50
6   0.33778 0.10 0.10 1.00 0.65 1.00 0.25
7   0.39556 0.10 0.10 1.00 0.65 1.00 0.25
10  0.03556 0.10 0.10 1.00 0.70 0.50 0.20
11  0.15556 0.10 0.10 0.75 0.60 0.50 0.30
3   4
6   0.07111 0.10 0.10 1.00 0.65 1.00 0.25
7   0.23111 0.10 0.10 1.00 0.65 1.00 0.25
10  0.26667 0.10 0.10 1.00 0.70 0.50 0.20
11  0.42667 0.10 0.10 0.75 0.60 0.50 0.30
4   4
6   0.04889 0.10 0.10 1.00 0.65 1.00 0.25
7   0.36444 0.10 0.10 1.00 0.65 1.00 0.25
10  0.02667 0.10 0.10 1.00 0.70 0.50 0.20
11  0.56    0.10 0.10 0.75 0.60 0.50 0.30
5   4
6   0.49333 0.10 0.10 1.00 0.65 1.00 0.25
7   0.37333 0.10 0.10 1.00 0.65 1.00 0.25
10  0.04889 0.10 0.10 1.00 0.70 0.50 0.20
11  0.06667 0.10 0.10 0.75 0.60 0.50 0.30
...

关于气象驱动数据的生成

这里放下我的用C++写的基于IDW的插值程序：https://github.com/Sibada/ForcingMaker。 欢迎debug。

在这里你需要的数据有：

• 各个站点的坐标数据。
• 各个站点的各种气象数据的时间序列（一个文件一种数据，一行一个时间点，一列一个站点）。
• VIC的各个网格的坐标数据。

可以选择输出二进制文件或者文本文件。 实测i7-4790K输出一万八个网格20年逐日的5类气象数据大概用时8到10分钟。

还想说的一些

要是还是遇到问题的话： 先把ArcGIS（或者其他GIS软件）上手 先把一门脚本语言上手 多看官方文档 多百度（谷歌） 就酱…