Grass Gis 是20世纪80年代初，美国军方建筑工程研究实验室（USA/CERL）的Bill Gorgan负责遴选一款具有土地管理、环境规划、环境评估的GIS软件。它有400多个地质分析模块，今天我们剖析一下r.lake算法。
1. 堰塞湖预测建模r.lake能做什么
如下图1左为一张地形栅格图，也是最常见的2D平面图。从图中可以看到许多的凹地盆地，当有山洪暴雨时就会形成堰塞湖。 Grass Gis 的r.lake就是能够预测堰塞湖形状大小的模块。例如，在图1左所示的红色圆圈位置开始注水当水位海拔高度为113.4米时，形成的堰塞湖的位置形状如图1右所示。Grass可以形成动画动态展示不断的提高水位后堰塞湖形态的变化。
我国长江三峡大坝水位设计海拔185米，蓄水库区淹没129座城镇，淹没耕地1.94万公顷。东起宜昌，西至重庆662.9Km沿长江两岸分水岭范围。在设计时就需要对185米高湖所能淹没的地区有非常精确的预测，如果沿岸有缺口造成计划外的地区被淹没就需要处理缺口。
Grass Gis 的r.lake 就可以完成这个事情，它的核心算法就是“湖心灌水洪泛算法“，它模拟了水在一点不断涌出时淹没了周边地区。
图1 堰塞湖预测模型演示
2. 数据降维简化模型
当一个点开始蓄水时，周边山岭的高度和形态都可以阻挡洪水的进一步蔓延，就像一个围栏一样。那怎样才能根据水位的海拔高度找出这些围栏的位置？首先想到的是三维建模以及空间关系计算等复杂的方式，其实不然。所有复杂的三维空间计算都要降维到二维平面计算，简化计算。
一个地区的遥感影像被处理成高程图也就是栅格图后，就变成了一个标准的M*N的数字矩阵。矩阵的每个点对于地图上的一个矩形区域，该点的值就是该区域的平均高度。图1就是这样的一个二维栅格图，只是每个像素被赋予不同的灰度颜色展现出来。一个复杂的三维地势图化简为矩阵后就变成了如图2a所示的二维数组，图中的数字为某块土地的海拔。
图2 栅格地图数据矩阵
当要求堰塞湖水位海拔9米时，满足条件的区域如图2b所示的三块区域。这样问题好像看起来非常简单了, 就是把数组遍历一遍找出<=9的区域不就行了？但是复杂的问题来了，当洪水爆发点在 6C 点，那9米的水只能在红色区域聚集成堰塞湖，黄色和蓝色都不应该被淹没，因为水被红色区域周围20多米高的地势挡住了。
现在的问题已经被简化成了一个简单的数据结构算法了。以 6C 点逐渐向往扩张，就像湖心有洪水不断涌出直到把周围的相邻8个格子中低于9米的区域淹没。6C为湖心的周边5B,5C,5D,6B,6D,7B,7C,7D被淹没了。接着以它们为湖心重复相同的操作，尝试淹没周围的8个邻居区域。不断重复直到没有新的可淹没区域被找到。这个算法很简单，就不给出代码了。
3. 以一个或多个堰塞湖区为洪水点
Grass Gis 的r.lake提供了更为复杂的建模。不是以某个点为洪水发生地，而是以某个湖区在涌出洪水，甚至多个湖区一起涌出。这样每次都以上一轮模拟输出的湖区为下次模拟计算的输入，不断增高湖水高度，达到动画效果输出。
这个算法思想和上面单点的算法相似，只是每轮需要扫描M*N个像素点，每次是旁边8个区域只要有一个有水的点就试图向自己淹没，直到一轮没有新的扩张点就停下来。
</pre><p><pre name="code" class="cpp">int pases = (int)(rows * cols) / 2;int water_level = .... // 输入的要求的湖面海拔高度int *in_terran; // row * col 输入的栅格地图，数值为地势的海拔高度int *in_wanter; // row * col 输入的已经有水的湖区矩阵. 0无水，>0 为湖水平面到湖底的深度int *out_water; // row * col 输出计算完毕的湖区矩阵. 初始为0无水，>0 为湖水平面到湖底的深度copy_seed_water(in_wanter, out_water); // 初始输出湖区为输入湖区int water_window[3][3]; // 某点为中心的3*3矩阵窗口，该窗口沿着地图每个点移动一遍为一趟for (pass = 0; pass < pases; pass++) { curcount = 0; /* 每个点扫描 */ for (row = 0; row < rows; row++) { for (col = 0; col < cols; col++) { //该点为中心去3*3矩阵 load_neibour_window_values(out_water, water_window, rows, cols, row, col); //该点周围8个邻居是否有水往自己这边淹没 if (is_near_water(water_window) == 1) { if (in_terran[row][col] < water_level) { // 该地势低，被淹没 out_water[row][col] = water_level - in_terran[row][col]; curcount++; } else { out_water[row][col] = 0; /* 该地势比湖水面高 */ } } } } if (curcount == lastcount) break; /* 没有发现新的淹没点，完成. */ lastcount = curcount; curcount = 0; }}return out_water;4. Grass Gis平台提供的服务支持
上面的算法没有看到任何Grass Gis 平台API的使用。其实栅格数据的输入，输出和图形化显示都是透明的交给了grass 平台了。向湖区的梯度渐变颜色表的建立就是调用API并绑定湖区地图的，但是开发人员并不知道该颜色表是怎样存储并绑定的。它提供的功能包括：
- 各种地质地图数据文件在内存中的直接访问
- 处理完成的地质数据向grass数据库的内存直接写
- 地质地图颜色表的生成并绑定
- 多个地质地图的三维叠加显示
- 统一管理用户设定的的数据空间范围
版权声明：本文为博主原创文章，未经博主允许不得转载。
版权声明：本文为博主原创文章，未经博主允许不得转载。
版权声明：本文为博主原创文章，未经博主允许不得转载。