首页 > 代码库 > 洪涝有源淹没算法及淹没结果分析
洪涝有源淹没算法及淹没结果分析
洪涝模拟仿真的实现方法主要有两种:一种是基于水动力学的洪水演进模型;另一种是基于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<Point> 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<Point>(); } /// <summary> /// 加载淹没区DEM,并创建淹没范围影像 /// </summary> /// <param name="m_DEMFilePath">DEM文件路径</param> /// <returns></returns> 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(); } /// <summary> /// 种子扩散算法淹没分析 /// </summary> /// <param name="m_Lat">种子点纬度</param> /// <param name="m_Log">种子点经度</param> /// <param name="m_FloodLevel">淹没水位</param> 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中的显示效果图:
欢迎大家留言交流。
洪涝有源淹没算法及淹没结果分析
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。