Flood active submergence algorithm and flood result analysis

There are two main ways to realize flood simulation: one is based on the hydrodynamic model of flood evolution, and the other is the flood inundation analysis based on Dem. The specific analysis is as follows:

I am a GIS practitioner, from our professional point of view, choose Dem-based flood analysis to do flood simulation. The method of flood inundation analysis based on DEM is mainly divided into active submerged and non-source submerged.

This blog uses an active flooding algorithm to achieve flood simulation, the algorithm is eight domain seed diffusion algorithm. The Floodsimulation class was written in C # version gdal, with all the source code:

  Class Floodsimulation {#region class member variable//point structure public struct-dot {public int          X          Line number public int Y;  Column number public int elevation; Pixel value (elevation value) public bool isflooded;        Flood Mark};                Private bool[,] isflood;  Flooded area labeled two-dimensional array, used to mark flooded raster private list<point> m_floodbufferlist;            Flood buffer stack public Dataset m_demdataset; Dem Data set public DataSet M_floodsimulateddataset;                     Flood flooded range DataSet public int m_xsize;                     Data x-directional raster number public int m_ysize;        Data y-Direction raster number public OSGeo.GDAL.Driver Driver;            Image format driver public int[] M_floodbuffer;          Fill buffer (flood flooded range) public int[] M_demdatabuffer;            Dem data (storage elevation value) public double m_areaflooded;            Surface area public double m_watervolume; Submerged water volume */* Here the Geotransform (image coordinate transformation parameter) is defined as: by the row of pixels is worth to its upper-left corner of the space coordinates of the operation parameterExample: The actual spatial coordinates of the upper-left corner of an image (P,l) point are: Xp = geotransform[0] + P * geotransform[1] + L * geotransform[2];                                                                     Yp = geotransform[3] + P * GEOTRANSFORM[4] + L * geotransform[5];           */public double[] m_adfgeotransform;            #endregion//Constructor public floodsimulation () {m_adfgeotransform = new double[6];                    M_floodbufferlist = new list<point> (); }///<summary>//Load the submerging range dem and create a flooded range image///</summary>//<param name= "M_demfil Epath ">dem file path </param>///<returns></returns> public void Loaddataset (string m_demfile            Path) {//read dem Data M_demdataset = Gdal.open (M_demfilepath, access.ga_readonly);            M_xsize = m_demdataset.rasterxsize;                        M_ysize = m_demdataset.rasterysize;        Get Image Coordinate transformation parameters    M_demdataset.getgeotransform (M_adfgeotransform); Read DEM data into memory Band M_demband = M_demdataset.getrasterband (1);            Gets the first band 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);            Flooded range fill Buffer m_floodbuffer = new int[m_xsize * M_ysize];            Isflood=new Bool[m_xsize,m_ysize]; Initializes a two-dimensional submerging range bool array for (int i = 0; i < m_xsize; i++) {for (int j = 0; J < M_ysiz E                J + +) {isflood[i, j] = false; }}//Create flood-flooded range image 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); }//Created in GdalTo create an image, you need to specify the format of the image to be created, and get the drive Driver = Gdal.getdriverbyname ("Gtiff") to the image format; Call the Creat function to create the image//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); Set image Properties M_floodsimulateddataset.setgeotransform (M_adfgeotransform); Image conversion Parameter m_floodsimulateddataset.setprojection (M_demdataset.getprojectionref ());  Projection//m_floodsimulateddataset.writeraster (0, 0, m_xsize, m_ysize, M_floodbuffer, M_xsize, m_ysize, 1, NULL, 0,                        0, 0); Output image 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> flood analysis of seed diffusion algorithm///</summary>//<param name= "M_lat" > Seed point Latitude </param>//<param name= "M_log" > Seed point Longitude </param&        Gt <param name= "M_floodlevel" > Flooded water level </param> public void Floodfill8direct (double m_lat,double M_log,doub            Le M_floodlevel) {//First get its row number according to the seed point latitude and longitude pfloodsourcepoint = new points ();            int x, y;            Geotoimagespace (M_adfgeotransform, M_log, M_lat, out X, out y);            Pfloodsourcepoint.x = X;                        Pfloodsourcepoint.y = Y;            Gets the seed point elevation value pfloodsourcepoint.elevation = getelevation (pfloodsourcepoint);            M_floodbufferlist.add (Pfloodsourcepoint);            Determine when there are no submerged points in the stack, and if there is a continuation search, no drowning analysis ends.                while (m_floodbufferlist.count!=0) {point pfloodsourcepoint_temp = m_floodbufferlist[0]; int rowx = pfloodsourcepoint_temp.                X int colmy = pfloodsourcepoint_temp.                               YThe mark can be submerged and removed from the flooded stack isflood[rowx, colmy] = true;                M_floodbuffer[getindex (pfloodsourcepoint_temp)] = 1;                                 M_floodbufferlist.removeat (0); Search for the connected domain for (int i = rowX-1; I < ROWX + 2; i++) to the 8 adjacent directions of the center grid cell. {FOR                    (int j = colmY-1; J < Colmy + 2; j + +)                        {//Determine if the grid boundary is reached 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); Searches for a raster cell that can be flooded and not marked if (temp_point. elevation<m_floodlevel| | Temp_point. Elevation <= pfloodsourcepoint_temp. Elevation) && isflood[temp_point. X, Temp_point.                      Y] = = False)      {//Add the eligible raster cells to the stack, mark as submerged, avoid repeating operations M_floodbufferlist.add (                                Temp_point); Isflood[temp_point. X, Temp_point.                                Y] = true;                            M_floodbuffer[getindex (temp_point)] = 1;            }                                                    }                                           }                }            }//Statistics flooded mesh number 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>//Output flooding range map///</sum mary> 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>///Get the x row y-column Raster Index </summary>//<param name= "point" > Grid points </param>//<returns> index of the point in the DEM memory array &L t;/returns> private int GetIndex (point point) {return point. y* m_xsize + point.        X }///<summary>//Get elevation values///</summary>//<param name= "M_point" > Grid point </p            aram>//<returns> elevation value </returns> private int getelevation (point m_point) {        Return M_demdatabuffer[getindex (M_point)];      }///<summary>//convert from pixel space to geo space  </summary>//<param name= "adfgeotransform" > Image coordinate Transformation parameters </param>//<param name= "Pi Xel "> Pixel row </param>//<param name=" line "> Pixel column </param>///<param name=" x ">X< /param>//<param name= "y" >Y</param> public void Imagetogeospace (double[] m_geotransform, in t pixel, int line, out double X, out double Y) {x = m_geotransform[0] + pixel * m_geotransform[1] + Lin            E * m_geotransform[2];        Y = m_geotransform[3] + pixel * M_geotransform[4] + line * m_geotransform[5]; }///<summary>//convert from geospatial to pixel space///</summary>//<param name= "adfgeotransfor M "> Image coordinate Change parameters </param>//<param name=" x ">x (longitude) </param>//<param name=" y ">y (latitude) &l t;/param>//<param name= "pixel" > Pixels </param>//<param name= "line" > pixels in the same column </param&        Gt 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]); }    }

The simulation results are displayed in Arcglobe:

Welcome to Message exchange.

Flood active submergence algorithm and flood result analysis