洪涝模拟仿真的实现方法主要有两种：一种是基于水动力学的洪水演进模型；还有一种是基于DEM的洪水淹没分析。详细分析例如以下：
我是GIS从业者，从我们的专业角度出发，选择基于DEM的洪水淹没分析来做洪涝的模拟仿真。而基于DEM的洪水淹没分析方法主要分为有源淹没和无源淹没。
本篇博客採用有源淹没算法实现洪涝的模拟，算法为八领域种子扩散算法。
採用C#版本号GDAL编写了FloodSimulation类，以下给出所有源码：
class FloodSimulation { #region 类成员变量 //点结构体 public struct Point { public int X; //行号 public int Y; //列号 public int Elevation; //像素值(高程值) public bool IsFlooded; //淹没标记 }; private bool[,] IsFlood; //淹没区域标记二维数组，用于标记淹没栅格 private List m_FloodBufferList; //淹没缓冲区堆栈 public Dataset m_DEMDataSet; //DEM数据集 public Dataset m_FloodSimulatedDataSet; //洪涝淹没范围数据集 public int m_XSize; //数据X方向栅格个数 public int m_YSize; //数据Y方向栅格个数 public OSGeo.GDAL.Driver driver; //影像格式驱动 public int[] m_FloodBuffer; //填充缓冲区（洪涝淹没范围） public int[] m_DEMdataBuffer; //DEM数据（存储高程值） public double m_AreaFlooded; //水面面积 public double m_WaterVolume; //淹没水体体积 /* 这里的GeoTransform(影像坐标变换參数)的定义是：通过像素所在的行列值得到其左上角点空间坐标的运算參数 比如:某图像上(P,L)点左上角的实际空间坐标为： Xp = GeoTransform[0] + P * GeoTransform[1] + L * GeoTransform[2]; Yp = GeoTransform[3] + P * GeoTransform[4] + L * GeoTransform[5]; */ public double[] m_adfGeoTransform; #endregion //构造函数 public FloodSimulation() { m_adfGeoTransform = new double[6]; m_FloodBufferList = new List (); } /// /// 载入淹没区DEM,并创建淹没范围影像 /// /// DEM文件路径 /// /// 种子扩散算法淹没分析 /// /// 种子点纬度 /// 种子点经度 /// 淹没水位 public void FloodFill8Direct(double m_Lat,double m_Log,double m_FloodLevel) { //首先依据种子点经纬度获取其所在行列号 Point pFloodSourcePoint = new Point(); int x, y; geoToImageSpace(m_adfGeoTransform, m_Log, m_Lat, out x, out y); pFloodSourcePoint.X = x; pFloodSourcePoint.Y = y; //获取种子点高程值 pFloodSourcePoint.Elevation = GetElevation(pFloodSourcePoint); m_FloodBufferList.Add(pFloodSourcePoint); //推断堆栈中时候还有未被淹没点，如有继续搜索。没有则淹没分析结束。 while (m_FloodBufferList.Count!=0) { Point pFloodSourcePoint_temp = m_FloodBufferList[0]; int rowX = pFloodSourcePoint_temp.X; int colmY = pFloodSourcePoint_temp.Y; //标记可淹没,并从淹没堆栈中移出 IsFlood[rowX, colmY] = true; m_FloodBuffer[getIndex(pFloodSourcePoint_temp)] = 1; m_FloodBufferList.RemoveAt(0); //向中心栅格单元的8个临近方向搜索连通域 for (int i = rowX - 1; i < rowX + 2; i++) { for (int j = colmY - 1; j < colmY + 2; j++) { //推断是否到达栅格边界 if (i<=m_XSize&&j<=m_YSize) { Point temp_point = new Point(); temp_point.X = i; temp_point.Y = j; temp_point.Elevation = GetElevation(temp_point); //搜索能够淹没且未被标记的栅格单元 if ((temp_point.Elevation<m_FloodLevel||temp_point.Elevation <= pFloodSourcePoint_temp.Elevation) && IsFlood[temp_point.X, temp_point.Y] == false) { //将符合条件的栅格单元增加堆栈，标记为淹没，避免反复运算 m_FloodBufferList.Add(temp_point); IsFlood[temp_point.X, temp_point.Y] = true; m_FloodBuffer[getIndex(temp_point)] = 1; } } } } } //统计淹没网格数 int count = 0; for (int i = 0; i < m_XSize; i++) { for (int j = 0; j < m_YSize; j++) { if (IsFlood[i,j]==true) { count++; } } } } /// <summary> /// 输出洪涝淹没范围图 /// </summary> public void OutPutFloodRegion() { m_FloodSimulatedDataSet.GetRasterBand(1).WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 0, 0); // m_FloodSimulatedDataSet.WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 1, null, 0, 0, 0); m_FloodSimulatedDataSet.GetRasterBand(1).FlushCache(); m_FloodSimulatedDataSet.FlushCache(); } // public void OutPutFloodedInfo() // { // } /// <summary> /// 获取第x行第y列栅格索引 /// </summary> /// <param name="point">栅格点</param> /// <returns>该点在DEM内存数组中的索引</returns> private int getIndex(Point point) { return point.Y* m_XSize + point.X; } /// <summary> /// 获取高程值 /// </summary> /// <param name="m_point">栅格点</param> /// <returns>高程值</returns> private int GetElevation(Point m_point) { return m_DEMdataBuffer[getIndex(m_point)]; } /// <summary> /// 从像素空间转换到地理空间 /// </summary> /// <param name="adfGeoTransform">影像坐标变换參数</param> /// <param name="pixel">像素所在行</param> /// <param name="line">像素所在列</param> /// <param name="x">X</param> /// <param name="y">Y</param> public void imageToGeoSpace(double[] m_GeoTransform, int pixel, int line, out double X, out double Y) { X = m_GeoTransform[0] + pixel * m_GeoTransform[1] + line * m_GeoTransform[2]; Y = m_GeoTransform[3] + pixel * m_GeoTransform[4] + line * m_GeoTransform[5]; } /// <summary> /// 从地理空间转换到像素空间 /// </summary> /// <param name="adfGeoTransform">影像坐标变化參数</param> /// <param name="x">X(经度)</param> /// <param name="y">Y（纬度）</param> /// <param name="pixel">像素所在行</param> /// <param name="line">像素所在列</param> public void geoToImageSpace(double[] m_GeoTransform, double x, double y, out int pixel, out int line) { line = (int)((y * m_GeoTransform[1] - x * m_GeoTransform[4] + m_GeoTransform[0] * m_GeoTransform[4] - m_GeoTransform[3] * m_GeoTransform[1]) / (m_GeoTransform[5] * m_GeoTransform[1] - m_GeoTransform[2] * m_GeoTransform[4])); pixel = (int)((x - m_GeoTransform[0] - line * m_GeoTransform[2]) / m_GeoTransform[1]); } }
模拟结果在ArcGlobe中的显示效果图：
欢迎大家留言交流。