学习记录十一：VIC模型—汇流

2023-09-01 17:35 作者:学习好难啊__ 0人读过 | 我要投稿

一、汇流原理

二、操作步骤（流向—产流比例—流程文件—水文站点位置—全局参数文件—汇流模型运行）

1.产流比例文件（fraction.txt）

文件内容：里面是边界网格被包含到实际的流域边界内的部分面积占该网格面积的比例，可以用ArcGIS制作。

（1）栅格转点+点转栅格：在ArcGIS中，导入研究区栅格数据"hhy_grid"，进行"栅格转点"+"点转栅格(pointed字段)"处理，得到分网格的栅格数据。

（2）栅格转面：取消勾选"简化面"，输出"HHY_rout_1.shp"

（3）栅格转面：导入"dem_one"数据，进行"栅格转面"处理，输出"HHY_rout_2.shp"

以上两个步骤输出数据的属性表里应有"面积"字段，可我输出的没有，重新运行，不修改输出路径，按GIS默认路径输出就没有这个问题，这是什么bug???

（4）相交：将"HHY_rout_1.shp"和"HHY_rout_2.shp"进行相交处理

（5）导出属性表：所需数据（红框内，数据表列数与教程中稍有出入，不影响结果）

（6）在RStudio中，打开代码"1_Fraction"，修改输入（上步属性表）输出路径后运行

（7）栅格转ASCII：在ArcGIS中，将研究区栅格数据"hhy_grid"进行栅格转ASCII处理

（8）在RStudio中，打开代码"2_precision"，修改输入输出路径、行列数后运行

Error1：

解决：创建ckgssj副本，命名为ckgssj1和ckgssj2，将ckgssj1表头删除，并将R中输入文件改为ckgssj1

Error2：同上

解决：修改行列数（写反了）

产流比例文件

2.流向文件（flowdir.txt）

文件内容：里面存放的是每个网格内水流的方向，用数宇1-8表示八个方位，这里一定不能用ArcGIS直接做流向，因为会丢失很多细节（重采样丢细节），而且容易造成封闭流向，所采用的方法是利用每个网格的最大汇流累积量代替高程，跑出来也基本上是对的。

（1）栅格转点+点转栅格：在ArcGIS中，导入研究区栅格数据"hhy_grid"，进行"栅格转点"+"点转栅格(pointed字段)"处理，得到分网格的栅格数据

（2）重采样：分辨率大小和dem的相同

（3）按掩膜提取

（4）栅格转点：对第二步"重采样"和第三步"按掩膜提取输出"的数据进行栅格转点处理（耗时较久，20min左右）

（5）添加连接：将第四步输出的数据用"pointid"字段进行"连接"

（6）导出属性表（10min）

（7）在RStudio中，打开代码"3_flowdir"，修改输入（上步属性表）输出路径后运行

（8）打开代码"4_precision"，修改输入输出路径、行列数后运行

（9）打开代码"5_flowdir2"，修改输入输出路径后运行

（10）将"F_div"中的表头拷贝至"hhy_F_dir1"

最终流向文件展示：

3.流程文件的制作（xmask）

（1）栅格转点+点转栅格+栅格转面：在ArcGIS中，导入研究区栅格数据"hhy_grid"，进行"栅格转点"+"点转栅格(pointed字段)"+"栅格转面"处理，分别输出xmask_1/2/3。

（2）在RStudio中，打开代码"7_Xmask"，先完成①②，后根据此修改代码输入输出路径等

①测量x_mask3中网格横、竖、对角长度（流向）

②拷贝flowdir输出结果"hhy_F_dir1"至x_mask文件夹，并删除其表头

③在"7_Xmask"中修改代码输入输出路径等，运行代码输出结果，给结果文件添加表头

4.确定水文站在模拟栅格中的位置

计算公式如下：根据已知水文点经纬度、xllcorner、yllcorner、cellsize，求RN(行数)和CN(列数)

其中xllcorner、yllcorner、cellsize（网格分辨率）值如下：

计算得出：

据此，制备水文站位置文件"F_F_stnloc.txt"

5.全局参数文件的制作

（1）将官网例子"stehekin\params\rout"中的全局参数文件"rout_input.STEHE"+"UH.all"和第1-4步制备的流向、产流比例、流程、水文站位置文件拷贝至"F:\VIC_modle\result\F_F\rout\param"中，将全局参数文件格式改为"txt"格式；将官网例子"route_1.0"文件夹拷贝至"F:\VIC_modle\result\F_F\rout"中

（2）修改全局参数文件输入输出路径

（3）将"rout_1.0"中的"Makefile"文件中第25行"FC=g77"修改为"FC=gfortran"

6.汇流模型的运行

（1）打开"Cygwin64 Terminal"软件，进行编译、运行，输出结果为"rout.exe"，将其拷贝至Cygwin安装软件目录bin里

（2）继续在"Cygwin64 Terminal"软件中，进行编译（全局参数）

Error1：

解决方法：参考https://blog.csdn.net/weixin_54697010/article/details/129163963?

修改：

Error2：

可能原因：产流（ArcGIS导出）、汇流（R）坐标不同

解决步骤：

①重采样：在ArcGIS中打开网格栅格数据，进行重采样

②栅格转ASCII

③修改汇流文件的坐标：用输出的ASCII文件的表头，替换以下三个文件的表头

④栅格转点：给转点数据的属性表添加两个字段"x""y"，并通过计算几何赋值，最后导出属性表

⑤重新制备forcing文件：打开制备气象驱动数据的最后一个代码"Conbination_PMNW"，修改经纬度参考坐标文件的路径，删除原有输出文件的数据，重新运行代码，将输出结果拷贝至"F:\VIC_modle\result\F_F\vic\forcing"

⑥修改土壤参数文件的坐标：打开"8modifySumSoilParam"代码，修改输入输出路径后运行

⑦重新制备产流文件：打开"cmd"命令提示符，运行vic模型

运行成功（比较耗时-1h）

⑧重新制备汇流文件：打开"cmd"命令提示符，运行rout模型

依然报错！！！

不要慌！不要急！仔细检查！继续修改！问题就是一个个慢慢解决的！！

Error3：水文站点位置问题

解决方法：经反复检查，在第四步"4.确定水文站在模拟栅格中的位置"中计算得出的数据没有错误，可能是数据稍有偏差，因此我尝试小幅度修改行列值，在修改为"52，40"时成功了

重新运行rout模型，水文站点位置问题成功解决！！

生成的过程线文件（TNH .uh_s）1686行，此前只有1行

为方便后续调参，可在"F_F_stnloc"中，将运行出的正确的过程线文件地址输入（替换"NONE"）

Error4（与Error2一样）：

自动识别出的经纬度坐标与径流输出经纬度不匹配（参数文件问题）制作的文件存在偏移

我通过对比报错经纬度和fluxes_原数据，发现：

①经度：0.0049——0.0050

②纬度：0.0016 —— 0.0017

解决方法：修改参数（流向、产流比例、流程文件表头的xllcorner和yllcorner）

（要确定此前用的同一个网格数据，否则调参无效，在Error2中已重做修改）

第一次修改：

失败原因：纬度小了0.0001（第一次改动幅度太大了，把小数点后面删了很多，吸取教训）

第二次修改：

失败原因：经度大了0.0001

第三次修改：

失败原因：经度小了0.0001

（本来以为这次成功了，前面非常丝滑，结果从第805个网格开始报错/(ㄒoㄒ)/~~）

第四次修改：

失败原因：和第三次一样，经度依然小了0.0001

第五次修改：

这次离成功很近了！！！只剩下25个not found ！！！

第六次修改......第n次修改......放弃......

我决定从根本上解决问题，因为每次都是那25行经度为_100.0417的问题。经检查，是经纬度原数据中"100.041650164299995"的问题，将其全部更改为"100.041640164299995"，即保证输出为_100.0416，然后重新制备产流文件（Error2中的第⑤⑥⑦步）

我想这次应该是没问题了，再次利用第五次修改中的坐标数据，进行汇流

没有报错！！！汇流成功！！！

7月中旬至今，经过间断性学习，历经重重报错，终于跑通了整个VIC模型，轻舟已过万重山，当然，后续还需要继续学习参数率定，以及需要从头来过再跑一遍VIC模型......因为

记得刚开始学习时，我有一个遗留问题，我搜索"VIC模型多站点"等关键词并没有得到答案，于是我天真的以为，VIC模型想输出多个站点数据，每个站点都要重新跑一遍，于是我只跑了第一个站点，直到完成汇流，我才在某篇文章中发现，汇流时可以输入多个站点数据......这意味着，这一个半月只是一次练习，在接下来的一周或者一个月，我得全部从头开始再跑一遍数据输出所有我需要的站点数据......哭了/(ㄒoㄒ)/~~......