基于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;//设定死的每秒单个格子上升多少mdouble 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;//可识别原始点};
}

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

上证指数评论数据情感分析（随时间变化）有代码数据可直接运行
目录情感分析结果: 编辑首先是获取股票评论数据的网站: 程序:
从入门到入土：基于C语言实现并发Web服务器|父进程子进程|代码展示
此博客仅用于记录个人学习进度,学识浅薄,若有错误观点欢迎评论区指出.欢迎各位前来交流.(部分材料来源网络,若有侵权,立即删除) 本人博客所有文章纯属学习之用,不涉及商业利益.不合适引用,自当删除! 若 ...
基于DEM的降雨淹没算法
数字高程模型(Digital Elevation Model),简称DEM,是通过有限的地形高程数据实现对地面地形的数字化模拟(即地形表面形态的数字化表达),它是用一组有序数值阵列形式表示地面高程的一 ...
matlab读取电子海图,基于dem数据叠加的航海雷达回波模拟方法
基于dem数据叠加的航海雷达回波模拟方法 [技术领域] [0001] 本发明涉及航海雷达的回波模拟系统,具体是一种基于DEM(Digital Elevation Model,数字高程模型)数据叠加的航 ...
主题模型综述：短文本、细粒度、加入先验知识、作者写作偏好、主题内涵随时间的变迁、融入词嵌入特性、语言模型加持
原文链接:https://www.zhihu.com/question/34801598/answer/765580727 主题模型当然有用咯,谁用谁知道!这次我来展示下它的7个"变种&qu ...
基于数据智能的区域教育大平台建设与应用实践
点击上方蓝字关注我们基于数据智能的区域教育大平台建设与应用实践贺相春, 郭绍青西北师范大学教育技术学院,甘肃兰州 730070 摘要:数据智能引领是新时期区域教育大平台建设与应用的重要方向 ...
基于DEM的沟壑特征分析
您的位置:第十二章基于DEM的沟壑特征分析第十二章 1.本章主题编号专题序号专题名称子专题号子专题名称子专题主要内容实验内容备注 12 基于DEM的沟壑特征分析 1 基于DEM的沟 ...
基于Simulink模拟具有两个目标的双基地雷达（附源码）
目录一.示例二.发射机三.目标四.接收机五.结果和显示六.总结七.程序此示例演示如何仿真具有两个目标的双基地雷达系统.双基地雷达的发射器和接收器不位于同一位置,而是沿着不同的路径移动. ...
《Splunk智能运维实战》——3.11　制作折线图显示项目浏览量和购买量随时间的变化...
本节书摘来自华章计算机<Splunk智能运维实战>一书中的第3章,第3.11节,作者 [美]乔史·戴昆(Josh Diakun),保罗R.约翰逊(Paul R. Johnson),德莱克· ...

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

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

最新文章

热门文章