基于DEM模拟淹没区域随时间推演算法代码展示

之前写的一篇博客讲基于dem模拟淹没区域随时间推演的算法：https://blog.csdn.net/wqy248/article/details/81119550
有些朋友留言希望看一下实现源码，所以在这篇文章中展示其中主要实现代码，完整工程资源见：https://download.csdn.net/download/wqy248/10968410
但是我想实际上没有下载完整工程的必要，因为每个人看完这篇文章都可以写出适用于自己项目的程序了。

首先需要申明的是，数据处理的原dem文件因为涉密缘故不能发在网上，大家可以使用自己的dem文件（tif格式），并且需要自己记录出水口大坝的位置（符合对应dem的投影坐标），实际程序迭代过程需要用到。

本次实现的整体流程：
1.根据dem生成另一张图，这张图记录了所有被淹没点在时间轴上的位置，类似下图所示：（之所以呈扇形排布，是因为模拟水流入平地处呈圆形扩散开来的态势，在山地区域就不会看到明显的扇状）

2.使用arcgis的工具箱arctoolbox中的工具，把上面生成的tif图片生产成不同时间点的shp矢量图，方法是重分类，比如指定时间为1（这里的1并不是指1s，而是一个相对概念），则得到重分类的阈值为1，上图value值小于等于1的为一类，大于1的为另一类，把小于等于1的那一类由栅格转矢量就得到一个面状图形，这个面状图形代表时间为1时水淹没到哪里，然后给时间+10，继续重复以上步骤，得到另一个shp，代表时间为11时水淹没到哪里，依次类推直到最大的时间点，得到的shp是面积最大的淹没区域，可以根据实际需要确定每次增加时间数量，越小则约流畅平滑，但是生成的shp图层越多，越大则反之。
3.使用gisapi加载到应用项目中（按照时间顺序依次加载不同的shp图层，就可以看到动态的淹没区域演变过程了）。

本文主要介绍第一步，即如何根据一张dem 的tif格式地图生成标注淹没时间点的栅格图。

首先确定大坝位置，大坝位置是种子点、出水口，我们要实现的是有源淹没，必须有一个出水的位置。打开arcmap，打开dem地图，我们标出大坝的位置（相连的像素点）：

把鼠标移动到每个像素点记录下这些像素点的x、y坐标（这是在dem图的坐标系下的xy坐标）。

C#程序代码如下：

using OSGeo.GDAL;
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.IO;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using System.Collections.Generic;
using System.Collections;
using Oracle.DataAccess.Client;
using System.Threading;
namespace WindowsFormsApplication1
{
    public partial class Form1 : Form
    {
        int m_XSize, m_YSize;
        double[] m_adfGeoTransform = new double[6];
        string projection;
        double[] m_FloodBuffer;
        List<PoiD> dabapois = new List<PoiD>();
        double peradddem = 2;//设定死的每秒单个格子上升多少m
        double perextend = 1;//设定死一次扩充1个格子
        public Form1()
        {
            InitializeComponent();
            //此处是大坝位置坐标，有多少记录多少，这里是假数据
            dabapois.Add(new PoiD(1111, 2222));
            dabapois.Add(new PoiD(2222, 33333));
            dabapois.Add(new PoiD(2222, 33333));
            dabapois.Add(new PoiD(2222, 33333));
            dabapois.Add(new PoiD(2222, 33333));
            dabapois.Add(new PoiD(2222, 33333));
        }
        
        private void Form1_Load(object sender, EventArgs e)
        {
            try
            {
                OSGeo.GDAL.Gdal.AllRegister();
                OSGeo.GDAL.Dataset dataSet = OSGeo.GDAL.Gdal.Open(@"******zidingyi\fillabove0.tif", Access.GA_ReadOnly);    //dem tif图片所在文件位置           
                m_XSize = dataSet.RasterXSize;
                m_YSize = dataSet.RasterYSize;
                m_FloodBuffer = new double[m_XSize * m_YSize];
                dataSet.GetRasterBand(1).ReadRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 0, 0);               
                dataSet.GetGeoTransform(m_adfGeoTransform);
                projection=dataSet.GetProjection();
            }
            catch (Exception err)
            {
                Console.WriteLine(err.Message);
            } 
        }
 private void button1_Click(object sender, EventArgs e)
        {
            double floodheight = Convert.ToDouble(textBox1.Text);
            FloodFill8Direct(floodheight);           
            MessageBox.Show("创建完毕");
        }
         private int getIndex(PointEleva point)
        {
            return point.Y * m_XSize + point.X;
        }
        private double GetElevation(PointEleva m_point)
        {
            int idx = getIndex(m_point);
            if(idx>=m_FloodBuffer.Length||idx<0)
            {
                return 10000000;
            }
            return m_FloodBuffer[idx];
        }  
          /// <summary>  
        /// 种子扩散算法淹没分析  
        /// </summary>  
        /// <param name="m_Lat">种子点纬度</param>  
        /// <param name="m_Log">种子点经度</param>  
        /// <param name="m_FloodLevel">淹没水位</param>  
        public void FloodFill8Direct(double m_FloodLevel)
        {
            List<PointEleva> m_FloodBufferList = new List<PointEleva>();  //淹没缓冲区堆栈
            bool[,] IsFlood;                //淹没区域标记二维数组，用于标记淹没栅格  
            double[,] IsFloodLevel;
            double maxfloodlevel = 1;
            IsFlood = new bool[m_XSize, m_YSize];
            IsFloodLevel = new double[m_XSize, m_YSize];
            for(int i=0;i<m_YSize;i++)
            {
                for (int j = 0; j < m_XSize; j++)
                {
                    IsFlood[j, i] = false;                    
                }
            }
            foreach(PoiD daba in dabapois)
            {
                //首先根据种子点经纬度获取其所在行列号  
                PointEleva pFloodSourcePoint = new PointEleva();
                int x, y, idx;
                geoToImageSpace(m_adfGeoTransform, m_XSize, m_YSize, daba.X, daba.Y, out x, out y, out idx);
                pFloodSourcePoint.X = x;
                pFloodSourcePoint.Y = y;
                pFloodSourcePoint.IsDaba = true;
                pFloodSourcePoint.originPoi = pFloodSourcePoint;
                pFloodSourcePoint.FloodLevel = pFloodSourcePoint.PerDistance;
                //获取种子点高程值  
                pFloodSourcePoint.Elevation = m_FloodBuffer[m_XSize * y + x];
                m_FloodBufferList.Add(pFloodSourcePoint);
            }
            //判断堆栈中是否还有未被淹没点，如有继续搜索，没有则淹没分析结束。  
            while (m_FloodBufferList.Count != 0)
            {
                PointEleva pFloodSourcePoint_temp = m_FloodBufferList[0];
                int rowX = pFloodSourcePoint_temp.X;
                int colmY = pFloodSourcePoint_temp.Y;
                bool isdaba = pFloodSourcePoint_temp.IsDaba;
                //标记可淹没,并从淹没堆栈中移出  
                IsFlood[rowX, colmY] = true;
                IsFloodLevel[rowX, colmY] = pFloodSourcePoint_temp.FloodLevel;
                //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 (Math.Sqrt((i - rowX) * (i - rowX) + (j - colmY) * (j - colmY)) > 3) continue;
                        //if ((i == rowX - 2 || i == rowX + 2) && j != colmY) continue;
                        //if ((i == rowX - 1 || i == rowX + 1) && (j == colmY - 2 || j == colmY + 2)) continue;
                        if (isdaba && (i + j) <= (rowX + colmY))
                            continue;
                        //判断是否到达栅格边界  
                        if (i < m_XSize && i >= 0 && j < m_YSize && j >= 0)
                        {
                            PointEleva temp_point = new PointEleva();
                            temp_point.X = i;
                            temp_point.Y = j;
                            temp_point.Elevation = GetElevation(temp_point);
                            temp_point.IsDaba = false;
                            temp_point.parentPoint = pFloodSourcePoint_temp;
                            if (temp_point.Elevation <= 0)
                                continue;
                            bool isflood = IsFlood[temp_point.X, temp_point.Y];
                            //搜索可以淹没且未被标记的栅格单元  
                            if ((temp_point.Elevation <= m_FloodLevel || temp_point.Elevation <= pFloodSourcePoint_temp.Elevation) && IsFlood[temp_point.X, temp_point.Y] == false)
                           // if (temp_point.Elevation <= m_FloodLevel || temp_point.Elevation <= pFloodSourcePoint_temp.Elevation)
                            {
                                IsFlood[temp_point.X, temp_point.Y] = true;
                                temp_point.IsFlooded = true;
                                PointEleva flagpoi = pFloodSourcePoint_temp;                                
                                PointEleva newflagpoi = pFloodSourcePoint_temp;
                                bool flag = true;
                                while (flag && newflagpoi.FloodLevel > 1)
                                {
                                    flagpoi = newflagpoi;
                                    newflagpoi = flagpoi.originPoi;
                                    flag = PointsCanReturn(IsFlood, new Point(temp_point.X, temp_point.Y), new Point(newflagpoi.X, newflagpoi.Y));  
                                }
                                temp_point.originPoi = flagpoi;                                
                                temp_point.FloodLevel = temp_point.originPoi.FloodLevel + temp_point.PerDistance;
                                //将符合条件的栅格单元加入堆栈，标记为淹没，避免重复运算  
                                m_FloodBufferList.Add(temp_point);
                                IsFloodLevel[temp_point.X, temp_point.Y] = temp_point.FloodLevel;
                                if (IsFloodLevel[temp_point.X, temp_point.Y] > maxfloodlevel)
                                {
                                    maxfloodlevel = IsFloodLevel[temp_point.X, temp_point.Y];
                                }                               
                            }
                        }
                    }
                }             
            }
            double[] waterbuffer = new double[m_FloodBuffer.Length];
            for (int i = 0; i < m_XSize; i++)
                for (int j = 0; j < m_YSize; j++)
                {
                    waterbuffer[j * m_XSize + i] = IsFloodLevel[i, j];
                }
            writegeotif(System.IO.Path.GetDirectoryName(System.Windows.Forms.Application.ExecutablePath) + "\\FloodSimulation\\water" + "_" + m_FloodLevel.ToString() + ".tif",
                  waterbuffer, m_adfGeoTransform, m_XSize, m_YSize);          
        } 
         private bool PointsCanReturn(bool[,] isflood,Point poi1,Point poi2)
        {            
            int dx = poi2.X - poi1.X;
            int dy = poi2.Y - poi1.Y;
            List<Point> poi = new List<Point>();
            int dtx = (dx==0)?0:dx/Math.Abs(dx);
            int dty = (dy==0)?0:dy/Math.Abs(dy);
            if(Math.Abs(dx)>Math.Abs(dy))
            {
                int xmid = poi1.X+dtx;
                while(xmid!=poi2.X)
                {
                    int ymid = (xmid - poi1.X) * dy / dx + poi1.Y;
                    poi.Add(new Point(xmid, ymid));
                    xmid += dtx;
                }
            }
            else
            {
                int ymid = poi1.Y + dty;
                while(ymid!=poi2.Y)
                {
                    int xmid = (ymid - poi1.Y) * dx / dy + poi1.X;
                    poi.Add(new Point(xmid, ymid));
                    ymid += dty;
                }               
            }
            foreach(Point p in poi)
            {
                if(isflood[p.X,p.Y]==false)
                {
                    return false;
                }
            }
            return true;
        }
          private void writegeotif(string path,double[] data,double[] transform,int width,int height)
        {
            //在GDAL中创建影像,先需要明确待创建影像的格式,并获取到该影像格式的驱动
            OSGeo.GDAL.Driver driver = Gdal.GetDriverByName("GTiff");
            //调用Creat函数创建影像
            if (System.IO.File.Exists(path))
            {
                System.IO.File.Delete(path);
            }
            Dataset m_FloodSimulatedDataSet = driver.Create(path, width, height, 1, DataType.GDT_Float32, null);
            //设置影像属性
            m_FloodSimulatedDataSet.SetGeoTransform(transform); //影像转换参数
            m_FloodSimulatedDataSet.SetProjection(projection); //投影            
            m_FloodSimulatedDataSet.GetRasterBand(1).WriteRaster(0, 0, width, height, data, width, height, 0, 0);
            m_FloodSimulatedDataSet.GetRasterBand(1).FlushCache();
            m_FloodSimulatedDataSet.FlushCache();
        }
        /// <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,int sizex,int sizey, double x, double y, out int pixel, out int line,out int totalindx)
        {
            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]);
            totalindx = sizex * line + pixel;
        }
}
 public class PoiD
    {
        public double X;
        public double Y;
        public PoiD(double x, double y)
        {
            X = x;
            Y = y;
        }
    }   
    //点结构体  
    public class PointEleva
    {
        public int X;          //行号  
        public int Y;          //列号  
        public double Elevation;  //像素值(高程值)  
        public bool IsFlooded; //淹没标记  
        public bool IsDaba;//标记是否为大坝出水口
        public double PerDistance
        {
            get
            {
                return Math.Sqrt((X - originPoi.X) * (X - originPoi.X) + (Y - originPoi.Y) * (Y - originPoi.Y));
            }
        }
        public double FloodLevel; //淹没层级
        public PointEleva parentPoint;//父亲节点 
        public PointEleva originPoi;//可识别原始点
        
    };  
}