洪涝有源淹没算法及淹没结果分析
洪涝模拟仿真的实现方法主要有两种:一种是基于水动力学的洪水演进模型;另一种是基于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中的显示效果图:
欢迎大家留言交流。
洪涝有源淹没算法及淹没结果分析相关推荐
- 有源淹没分析arcgis_洪涝有源淹没算法及淹没结果分析
洪涝模拟仿真的实现方法主要有两种:一种是基于水动力学的洪水演进模型:还有一种是基于DEM的洪水淹没分析.详细分析例如以下: 我是GIS从业者,从我们的专业角度出发,选择基于DEM的洪水淹没分析来做洪涝 ...
- 洪涝有源淹没算法及淹没结果分析【转】
http://blog.csdn.net/giser_whu/article/details/41288761 洪涝模拟仿真的实现方法主要有两种:一种是基于水动力学的洪水演进模型:另一种是基于DEM的 ...
- 基于DEM的降雨淹没算法
数字高程模型(Digital Elevation Model),简称DEM,是通过有限的地形高程数据实现对地面地形的数字化模拟(即地形表面形态的数字化表达),它是用一组有序数值阵列形式表示地面高程的一 ...
- 图像处理中,SIFT,FAST,MSER,STAR等特征提取算法的比较与分析(利用openCV实现)
图像处理中,SIFT,FAST,MSER,STAR等特征提取算法的比较与分析(利用openCV实现) 本文实验为自己原创,转载请注明出处. 本人为研究生,最近的研究方向是物体识别.所以就将常用的几种特 ...
- SURF算法与源码分析、下
FROM: http://www.cnblogs.com/ronny/p/4048213.html 上一篇文章 SURF算法与源码分析.上 中主要分析的是SURF特征点定位的算法原理与相关OpenCV ...
- 动图图解C语言插入排序算法,含代码分析
C语言文章更新目录 C语言学习资源汇总,史上最全面总结,没有之一 C/C++学习资源(百度云盘链接) 计算机二级资料(过级专用) C语言学习路线(从入门到实战) 编写C语言程序的7个步骤和编程机制 C ...
- 动图图解C语言选择排序算法,含代码分析
C语言文章更新目录 C语言学习资源汇总,史上最全面总结,没有之一 C/C++学习资源(百度云盘链接) 计算机二级资料(过级专用) C语言学习路线(从入门到实战) 编写C语言程序的7个步骤和编程机制 C ...
- c语言将十进制转化为二进制算法_base64算法初探即逆向分析
算法分析 虽说base64严格意义上来说并不能算是加密算法,但的确应用方面来说还算是比较广,在CTF的算法逆向中Base系列算是也比较常见的,萌新刚开始学算法,就以base64为例,对该算法进行一个简 ...
- 数据结构与算法--代码鲁棒性案例分析
代码鲁棒性 鲁棒是robust的音译,就是健壮性.指程序能够判断输入是否符合规范,对不合要求的输入能够给出合理的结果. 容错性是鲁棒的一个重要体现.不鲁棒的代码发生异常的时候,会出现不可预测的异常,或 ...
最新文章
- 合并果子(贪心,优先队列)
- MCMC+马尔科夫链蒙特卡罗
- Libidn 简介 对国际化域名进行编码和解码
- 重新想象 Windows 8 Store Apps (49) - 输入: 获取输入设备信息, 虚拟键盘, Tab 导航, Pointer, Tap, Drag, Drop...
- Docker 1.10版本发布
- 小师妹学JVM之:cache line对代码性能的影响
- java永久冻结_Java如何解决脆弱基类(基类被冻结)问题
- 电脑机器人_视频|电话积分换平板电脑和扫地机器人?女子拿回家后……-
- srsLTE源码学习:生成多播信道表gen_mch_tables
- [附源码]java毕业设计社区医院电子病历系统
- 【MMDetection3D】基于单目(Monocular)的3D目标检测入门实战
- 云计算 | 浅议云计算发展趋势
- OLED之U8g2中文库使用
- 深入浅出学习Linux(基础知识一)
- 区块链是什么通俗解释?
- 神奇葩! 史上最牛的博士论文答辩
- mysql表情符存储设置
- 教你在微信头像上加口号,很实用!
- Mysql数据库经验总结
- 海康威视摄像头使用网线连网输入IP地址跳转显示“网站处于联机状态,但未对联机尝试做出反应”
热门文章
- 富人越富,穷人越穷,我为什么反对PoS
- SpringBoot项目中使用set方法后,自动保存问题
- python爬虫学习笔记分析Ajax爬取果壳网文章
- TCP网络编程之chat聊天室
- 计算机应用研究、计算机工程与应用、计算机科学与探索投稿经验
- 国产猫粮高端化难题不少,网易天成拿什么出众?
- 【工具介绍】fastcopy的下载与使用方法,可用于硬盘对拷
- 计算机网络MOOC期末考试答案与解析
- No qualifying bean of type 'com.xxx.xx.service.xxService' available: expected at leas
- python期末试题汇总