1. 程式人生 > >洪澇有源淹沒演算法及淹沒結果分析

洪澇有源淹沒演算法及淹沒結果分析

洪澇模擬模擬的實現方法主要有兩種:一種是基於水動力學的洪水演進模型;另一種是基於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中的顯示效果圖:

123

--------------------- 
作者:召喚師峽谷 
來源:CSDN 
原文:https://blog.csdn.net/giser_whu/article/details/41288761 
版權宣告:本文為博主原創文章,轉載請附上博文連結!