洪涝有源淹没算法及淹没结果分析【转】

http://blog.csdn.net/giser_whu/article/details/41288761

洪涝模拟仿真的实现方法主要有两种：一种是基于水动力学的洪水演进模型；另一种是基于DEM的洪水淹没分析。具体分析如下：

我是GIS从业者，从我们的专业角度出发，选择基于DEM的洪水淹没分析来做洪涝的模拟仿真。而基于DEM的洪水淹没分析方法主要分为有源淹没和无源淹没。

本篇博客采用有源淹没算法实现洪涝的模拟，算法为八领域种子扩散算法。采用C#版本GDAL编写了FloodSimulation类，下面给出全部源代码：

[csharp] view plain copy
  1.   class FloodSimulation  
  2.     {  
  3.         #region 类成员变量  
  4.   
  5.         //点结构体  
  6.         public struct Point  
  7.         {  
  8.             public int X;          //行号  
  9.             public int Y;          //列号  
  10.             public int Elevation;  //像素值(高程值)  
  11.             public bool IsFlooded; //淹没标记  
  12.   
  13.         };  
  14.         private bool[,] IsFlood;                //淹没区域标记二维数组，用于标记淹没栅格  
  15.         private List<Point> m_FloodBufferList;  //淹没缓冲区堆栈  
  16.         
  17.         public Dataset m_DEMDataSet;            //DEM数据集  
  18.         public Dataset m_FloodSimulatedDataSet; //洪涝淹没范围数据集  
  19.         public int m_XSize;                     //数据X方向栅格个数  
  20.         public int m_YSize;                     //数据Y方向栅格个数  
  21.         public OSGeo.GDAL.Driver driver;        //影像格式驱动  
  22.         public int[] m_FloodBuffer;            //填充缓冲区（洪涝淹没范围）  
  23.         public int[] m_DEMdataBuffer;          //DEM数据（存储高程值）   
  24.   
  25.         public double m_AreaFlooded;            //水面面积  
  26.         public double m_WaterVolume;            //淹没水体体积  
  27.         /* 这里的GeoTransform(影像坐标变换参数)的定义是：通过像素所在的行列值得到其左上角点空间坐标的运算参数 
  28.             例如:某图像上(P,L)点左上角的实际空间坐标为： 
  29.             Xp = GeoTransform[0] + P * GeoTransform[1] + L * GeoTransform[2]; 
  30.             Yp = GeoTransform[3] + P * GeoTransform[4] + L * GeoTransform[5];                                                                     */  
  31.         public double[] m_adfGeoTransform;     
  32.  
  33.         #endregion  
  34.           
  35.         //构造函数  
  36.         public FloodSimulation()  
  37.         {  
  38.             m_adfGeoTransform = new double[6];  
  39.             m_FloodBufferList = new List<Point>();  
  40.               
  41.         }  
  42.   
  43.         /// <summary>  
  44.         /// 加载淹没区DEM,并创建淹没范围影像  
  45.         /// </summary>  
  46.         /// <param name="m_DEMFilePath">DEM文件路径</param>  
  47.         /// <returns></returns>  
  48.         public void loadDataSet(string m_DEMFilePath)  
  49.         {  
  50.             //读取DEM数据  
  51.             m_DEMDataSet = Gdal.Open(m_DEMFilePath, Access.GA_ReadOnly);  
  52.             m_XSize = m_DEMDataSet.RasterXSize;  
  53.             m_YSize = m_DEMDataSet.RasterYSize;  
  54.               
  55.             //获取影像坐标转换参数  
  56.             m_DEMDataSet.GetGeoTransform(m_adfGeoTransform);   
  57.             //读取DEM数据到内存中  
  58.             Band m_DEMBand = m_DEMDataSet.GetRasterBand(1); //获取第一个波段  
  59.             m_DEMdataBuffer = new int [m_XSize * m_YSize];  
  60.             m_DEMBand.ReadRaster(0, 0, m_XSize, m_YSize, m_DEMdataBuffer, m_XSize, m_YSize, 0, 0);  
  61.   
  62.             //淹没范围填充缓冲区  
  63.             m_FloodBuffer = new int[m_XSize * m_YSize];  
  64.             IsFlood=new bool[m_XSize,m_YSize];  
  65.   
  66.             //初始化二维淹没区bool数组  
  67.             for (int i = 0; i < m_XSize; i++)  
  68.             {  
  69.                 for (int j = 0; j < m_YSize; j++)  
  70.                 {  
  71.                     IsFlood[i, j] = false;  
  72.                 }  
  73.             }  
  74.   
  75.             //创建洪涝淹没范围影像  
  76.             string m_FloodImagePath = System.IO.Path.GetDirectoryName(System.Windows.Forms.Application.ExecutablePath) + "\FloodSimulation\FloodedRegion.tif";  
  77.             if (System.IO.File.Exists(m_FloodImagePath))  
  78.             {  
  79.                 System.IO.File.Delete(m_FloodImagePath);  
  80.             }  
  81.   
  82.             //在GDAL中创建影像,先需要明确待创建影像的格式,并获取到该影像格式的驱动  
  83.             driver = Gdal.GetDriverByName("GTiff");  
  84.             //调用Creat函数创建影像  
  85.             // m_FloodSimulatedDataSet = driver.CreateCopy(m_FloodImagePath, m_DEMDataSet, 1, null, null, null);  
  86.             m_FloodSimulatedDataSet = driver.Create(m_FloodImagePath, m_XSize, m_YSize, 1, DataType.GDT_Float32, null);  
  87.             //设置影像属性  
  88.             m_FloodSimulatedDataSet.SetGeoTransform(m_adfGeoTransform); //影像转换参数  
  89.             m_FloodSimulatedDataSet.SetProjection(m_DEMDataSet.GetProjectionRef()); //投影  
  90.             //m_FloodSimulatedDataSet.WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 1, null, 0, 0, 0);  
  91.               
  92.             //输出影像  
  93.             m_FloodSimulatedDataSet.GetRasterBand(1).WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 0, 0);  
  94.             m_FloodSimulatedDataSet.GetRasterBand(1).FlushCache();  
  95.             m_FloodSimulatedDataSet.FlushCache();  
  96.   
  97.         }  
  98.   
  99.         /// <summary>  
  100.         /// 种子扩散算法淹没分析  
  101.         /// </summary>  
  102.         /// <param name="m_Lat">种子点纬度</param>  
  103.         /// <param name="m_Log">种子点经度</param>  
  104.         /// <param name="m_FloodLevel">淹没水位</param>  
  105.         public void FloodFill8Direct(double m_Lat,double m_Log,double m_FloodLevel)  
  106.         {  
  107.             //首先根据种子点经纬度获取其所在行列号  
  108.             Point pFloodSourcePoint = new Point();  
  109.             int x, y;  
  110.             geoToImageSpace(m_adfGeoTransform, m_Log, m_Lat, out x, out y);  
  111.             pFloodSourcePoint.X = x;  
  112.             pFloodSourcePoint.Y = y;  
  113.               
  114.             //获取种子点高程值  
  115.             pFloodSourcePoint.Elevation = GetElevation(pFloodSourcePoint);  
  116.             m_FloodBufferList.Add(pFloodSourcePoint);  
  117.   
  118.             //判断堆栈中时候还有未被淹没点，如有继续搜索，没有则淹没分析结束。  
  119.             while (m_FloodBufferList.Count!=0)  
  120.             {  
  121.                 Point pFloodSourcePoint_temp = m_FloodBufferList[0];  
  122.                 int rowX = pFloodSourcePoint_temp.X;  
  123.                 int colmY = pFloodSourcePoint_temp.Y;  
  124.                  
  125.                 //标记可淹没,并从淹没堆栈中移出  
  126.                 IsFlood[rowX, colmY] = true;  
  127.                 m_FloodBuffer[getIndex(pFloodSourcePoint_temp)] = 1;  
  128.                 m_FloodBufferList.RemoveAt(0);   
  129.                   
  130.                 //向中心栅格单元的8个临近方向搜索连通域  
  131.                 for (int i = rowX - 1; i < rowX + 2; i++)  
  132.                 {  
  133.                     for (int j = colmY - 1; j < colmY + 2; j++)  
  134.                     {  
  135.                         //判断是否到达栅格边界  
  136.                         if (i<=m_XSize&&j<=m_YSize)  
  137.                         {  
  138.                             Point temp_point = new Point();  
  139.                             temp_point.X = i;  
  140.                             temp_point.Y = j;  
  141.                             temp_point.Elevation = GetElevation(temp_point);  
  142.                             //搜索可以淹没且未被标记的栅格单元  
  143.                             if ((temp_point.Elevation<m_FloodLevel||temp_point.Elevation <= pFloodSourcePoint_temp.Elevation) && IsFlood[temp_point.X, temp_point.Y] == false)  
  144.                             {  
  145.                                 //将符合条件的栅格单元加入堆栈，标记为淹没，避免重复运算  
  146.                                 m_FloodBufferList.Add(temp_point);  
  147.                                 IsFlood[temp_point.X, temp_point.Y] = true;  
  148.                                 m_FloodBuffer[getIndex(temp_point)] = 1;  
  149.                             }  
  150.                               
  151.                         }  
  152.                          
  153.                     }  
  154.                 }  
  155.             }  
  156.             //统计淹没网格数  
  157.             int count = 0;  
  158.             for (int i = 0; i < m_XSize; i++)  
  159.             {  
  160.                 for (int j = 0; j < m_YSize; j++)  
  161.                 {  
  162.                     if (IsFlood[i,j]==true)  
  163.                     {  
  164.                         count++;  
  165.                     }  
  166.                 }  
  167.             }  
  168.         
  169.         }  
  170.           
  171.         /// <summary>  
  172.         /// 输出洪涝淹没范围图  
  173.         /// </summary>  
  174.         public void OutPutFloodRegion()  
  175.         {  
  176.             m_FloodSimulatedDataSet.GetRasterBand(1).WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 0, 0);  
  177.             // m_FloodSimulatedDataSet.WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 1, null, 0, 0, 0);  
  178.             m_FloodSimulatedDataSet.GetRasterBand(1).FlushCache();  
  179.             m_FloodSimulatedDataSet.FlushCache();  
  180.         }  
  181.   
  182. //         public void OutPutFloodedInfo()  
  183. //         {  
  184. //         }  
  185.         /// <summary>  
  186.         /// 获取第x行第y列栅格索引  
  187.         /// </summary>  
  188.         /// <param name="point">栅格点</param>  
  189.         /// <returns>该点在DEM内存数组中的索引</returns>  
  190.         private int getIndex(Point point)  
  191.         {  
  192.             return  point.Y* m_XSize + point.X;  
  193.         }  
  194.   
  195.         /// <summary>  
  196.         /// 获取高程值  
  197.         /// </summary>  
  198.         /// <param name="m_point">栅格点</param>  
  199.         /// <returns>高程值</returns>  
  200.         private int GetElevation(Point m_point)  
  201.         {  
  202.             return m_DEMdataBuffer[getIndex(m_point)];  
  203.   
  204.         }  
  205.   
  206.         /// <summary>  
  207.         /// 从像素空间转换到地理空间  
  208.         /// </summary>  
  209.         /// <param name="adfGeoTransform">影像坐标变换参数</param>  
  210.         /// <param name="pixel">像素所在行</param>  
  211.         /// <param name="line">像素所在列</param>  
  212.         /// <param name="x">X</param>  
  213.         /// <param name="y">Y</param>  
  214.         public void imageToGeoSpace(double[] m_GeoTransform, int pixel, int line, out double X, out double Y)  
  215.         {  
  216.             X = m_GeoTransform[0] + pixel * m_GeoTransform[1] + line * m_GeoTransform[2];  
  217.             Y = m_GeoTransform[3] + pixel * m_GeoTransform[4] + line * m_GeoTransform[5];  
  218.         }  
  219.   
  220.         /// <summary>  
  221.         /// 从地理空间转换到像素空间  
  222.         /// </summary>  
  223.         /// <param name="adfGeoTransform">影像坐标变化参数</param>  
  224.         /// <param name="x">X(经度)</param>  
  225.         /// <param name="y">Y（纬度）</param>  
  226.         /// <param name="pixel">像素所在行</param>  
  227.         /// <param name="line">像素所在列</param>  
  228.         public void geoToImageSpace(double[] m_GeoTransform, double x, double y, out int pixel, out int line)  
  229.         {  
  230.             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]));  
  231.             pixel = (int)((x - m_GeoTransform[0] - line * m_GeoTransform[2]) / m_GeoTransform[1]);  
  232.         }  
  233.   
  234.   
  235.     }  

模拟结果在ArcGlobe中的显示效果图：

欢迎大家留言交流。

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/mazhenyu/p/8119019.html