洪涝模拟仿真的实现方法主要有两种：一种是基于水动力学的洪水演进模型；还有一种是基于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 loadDataSet(string m_DEMFilePath)

{

//读取DEM数据

m_DEMDataSet = Gdal.Open(m_DEMFilePath, Access.GA_ReadOnly);

m_XSize = m_DEMDataSet.RasterXSize;

m_YSize = m_DEMDataSet.RasterYSize;

//获取影像坐标转换參数

m_DEMDataSet.GetGeoTransform(m_adfGeoTransform);

//读取DEM数据到内存中

Band m_DEMBand = m_DEMDataSet.GetRasterBand(1); //获取第一个波段

m_DEMdataBuffer = new int [m_XSize * m_YSize];

m_DEMBand.ReadRaster(0, 0, m_XSize, m_YSize, m_DEMdataBuffer, m_XSize, m_YSize, 0, 0);

//淹没范围填充缓冲区

m_FloodBuffer = new int[m_XSize * m_YSize];

IsFlood=new bool[m_XSize,m_YSize];

//初始化二维淹没区bool数组

for (int i = 0; i < m_XSize; i++)

{

for (int j = 0; j < m_YSize; j++)

{

IsFlood[i, j] = false;

}

}

//创建洪涝淹没范围影像

string m_FloodImagePath = System.IO.Path.GetDirectoryName(System.Windows.Forms.Application.ExecutablePath) + "\\FloodSimulation\\FloodedRegion.tif";

if (System.IO.File.Exists(m_FloodImagePath))

{

System.IO.File.Delete(m_FloodImagePath);

}

//在GDAL中创建影像,先须要明白待创建影像的格式,并获取到该影像格式的驱动

driver = Gdal.GetDriverByName("GTiff");

//调用Creat函数创建影像

// m_FloodSimulatedDataSet = driver.CreateCopy(m_FloodImagePath, m_DEMDataSet, 1, null, null, null);

m_FloodSimulatedDataSet = driver.Create(m_FloodImagePath, m_XSize, m_YSize, 1, DataType.GDT_Float32, null);

//设置影像属性

m_FloodSimulatedDataSet.SetGeoTransform(m_adfGeoTransform); //影像转换參数

m_FloodSimulatedDataSet.SetProjection(m_DEMDataSet.GetProjectionRef()); //投影

//m_FloodSimulatedDataSet.WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 1, null, 0, 0, 0);

//输出影像

m_FloodSimulatedDataSet.GetRasterBand(1).WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 0, 0);

m_FloodSimulatedDataSet.GetRasterBand(1).FlushCache();

m_FloodSimulatedDataSet.FlushCache();

}

///

/// 种子扩散算法淹没分析

///

/// 种子点纬度

/// 种子点经度

/// 淹没水位

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_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++;

}

}

}

}

///

/// 输出洪涝淹没范围图

///

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()

// {

// }

///

/// 获取第x行第y列栅格索引

///

/// 栅格点

/// 该点在DEM内存数组中的索引

private int getIndex(Point point)

{

return point.Y* m_XSize + point.X;

}

///

/// 获取高程值

///

/// 栅格点

/// 高程值

private int GetElevation(Point m_point)

{

return m_DEMdataBuffer[getIndex(m_point)];

}

///

/// 从像素空间转换到地理空间

///

/// 影像坐标变换參数

/// 像素所在行

/// 像素所在列

/// X

/// Y

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];

}

///

/// 从地理空间转换到像素空间

///

/// 影像坐标变化參数

/// X(经度)

/// Y(纬度)

/// 像素所在行

/// 像素所在列

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中的显示效果图：

欢迎大家留言交流。