基于DEM数据的河流提取

采用的是Windows窗体应用，本文仅提供主体代码。
算法思想：
1、读取数据，并采用BitMap格式进行可视化，高程越高，越接近白色
2、计算坡度值，并采用分层设色法可视化
3、河流提取
（1）首先进行填洼操作，本程序中采用了Moran和Vezina算法，首先用一极大高程水面将原始DEM数据表面淹没，然后通过迭代
去除 DEM上多余的水，最后得到的高程就是填洼后的高程。
（2）水流方向提取采用了D8算法
（3）汇水量计算，主要涉及到一个递归算法。

1、数据读取与可视化

        //文件菜单中打开文件
        private void 打开OToolStripMenuItem_Click(object sender, EventArgs e)
        {
            this.openFileDialog1.ShowDialog();
            if (this.openFileDialog1.FileName != "")
            {
                this.filename = this.openFileDialog1.FileName;
            }
         }
             //获取原始数据，并将DEM矩阵转换成位图BITMAP输出
        public void 原始数据绘图ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            //读取文件
            StreamReader streamReadfile = new StreamReader(this.filename);
            string data = streamReadfile.ReadToEnd();
            string[] strLine = data.Split('\n');//换行符分割字符串  
            NRows = Int32.Parse(this.textBox1.Text);
            NCols = Int32.Parse(this.textBox2.Text);

            //读取DEM数据，将高程值保存在二维数组中
            DEM = new double[NRows, NCols];
            string[] str1;
            int rgb;//寄存灰度
            int i; int j;
            int MaxBit = 2000;//最大高程值
            int MinBit = 0;   //最小高程值
            Bitmap Primary = new Bitmap(NRows, NCols, PixelFormat.Format32bppRgb);
            //文本数据转换成二维double数据
            {
                for (i = 0; i < NRows; i++)
                {
                    str1 = strLine[i].Split(',');
                    for (j = 0; j < NCols; j++)
                    {
                        DEM[i, j] = Convert.ToDouble(str1[j]);
                    }
                }
            }
           //转换成BMP，高程值越大，颜色越趋向于白色
            {
                for (i = 0; i < NRows; i++)
                {
                    for (j = 0; j < NCols; j++)
                    {
                        rgb = Convert.ToInt32(255 * (DEM[i, j] - MinBit) / (MaxBit - MinBit));
                        Color pixelColor = Color.FromArgb(rgb, rgb, rgb);
                        Primary.SetPixel(i, j, pixelColor);
                    }
                }
            }
            //输出位图
            {
                Graphics g = this.CreateGraphics();
                g.DrawImage(Primary, 50, 60);
                g.Dispose();
            }
        }
        

2、坡度的计算

 //计算坡度，并用分层设色法输出
        private void 坡度toolStripMenuItem_Click(object sender, EventArgs e)
        {
            int cellsize = Int32.Parse(this.textBox3.Text); //获取数据分辨率
            double[,] Slope = new double[NRows, NCols];    //坡度矩阵
            double fx, fy; //x方向和y方向
            int i, j;

            //计算坡度
            {
                for (i = 1; i < NRows - 1; i++)
                {
                    for (j = 1; j < NCols - 1; j++)
                    {

                        fx = (DEM[i - 1, j + 1] - DEM[i - 1, j - 1] + 2 * (DEM[i, j + 1] - DEM[i, j - 1]) + DEM[i + 1, j + 1] - DEM[i + 1, j - 1]) / (8 * cellsize);
                        fy = (DEM[i + 1, j - 1] - DEM[i - 1, j - 1] + 2 * (DEM[i + 1, j] - DEM[i - 1, j]) + DEM[i + 1, j + 1] - DEM[i - 1, j + 1]) / (8 * cellsize);
                        Slope[i, j] = Math.Atan(Math.Sqrt(fx * fx + fy * fy)) * 57.29578;//转换成角度制
                    }
                }
            }

            Bitmap slope = new Bitmap(NRows - 2, NCols - 2, PixelFormat.Format32bppRgb); //slope为坡度位图


            //转换成BMP,用不同的颜色表示不同的坡度
            int red; int green; int blue;//RGB色彩表示法
            {
                for (i = 1; i < NRows - 1; i++)
                {
                    for (j = 1; j < NCols - 1; j++)
                    {

                        if (Slope[i, j] < 2)
                        { red = 56; green = 168; blue = 0; }
                        else if (Slope[i, j] < 4)
                        { red = 139; green = 209; blue = 0; }
                        else if (Slope[i, j] < 6)
                        { red = 255; green = 255; blue = 0; }
                        else if (Slope[i, j] < 8)
                        { red = 255; green = 180; blue = 0; }
                        else if (Slope[i, j] < 10)
                        { red = 255; green = 90; blue = 0; }
                        else
                        { red = 255; green = 0; blue = 0; }
                        Color pixelColor = Color.FromArgb(red, green, blue);
                        slope.SetPixel(i - 1, j - 1, pixelColor);
                    }
                }
            }

            //输出坡度图
            {
                Graphics g = this.CreateGraphics();
                g.DrawImage(slope, 650, 60);
                g.Dispose();
            }
        }

3、河流提取

        private void 河流提取toolStripMenuItem_Click(object sender,EventArgs e)
        {
            int[,] river = Water(D8(Fill(DEM)));
            Bitmap River = new Bitmap(NRows, NCols, PixelFormat.Format32bppRgb);
            int Limit = Convert.ToInt32(this.toolStripComboBox1.Text);
            for (int i = 0; i < NRows; i++)
            {
                for (int j = 0; j < NCols; j++)
                {
                    if (river[i, j] > Limit)
                    {
                        Color pixelColor = Color.FromArgb(255, 255, 255);
                        River.SetPixel(i, j, pixelColor);
                    }
                    
                }
            }
            //输出位图
            {
                Graphics g = this.CreateGraphics();
                g.DrawImage(River, 300, 200);
                g.Dispose();
            }

        }


        //填充洼地算法，采用的是Moran和Vezina算法
        public double[,] Fill(double[,]DEM)
        {
            int H = 2000;//设置水面高度
            double[,] Fill = new double[NRows, NCols];

            //初始化
            for (int i = 0; i < NRows; i++)
            {
                for (int j = 0; j < NCols; j++)
                {
                    if (i == 0 || i == NRows - 1 || j == 0 || j == NCols - 1)
                    {
                        Fill[i, j] = DEM[i, j];
                    }
                    else
                    {
                        Fill[i, j] = H;
                    }
                }
            }

            bool Y_N;

            //采用迭代思想，去除多余的水
            do
            {
                Y_N = false;
                for (int i = 0; i < NRows; i++)
                {
                    for (int j = 0; j < NCols; j++)
                    {
                        if (Fill[i, j] > DEM[i, j])
                        {
                            for (int N = -1; N <= 1; N++)
                            {
                                for (int M = -1; M <= 1; M++)
                                {
                                    if ((N != 0) || (M != 0))
                                    {
                                        if (Fill[i, j] > Fill[i + N, j + M] + 0.00000001)
                                        {
                                            Fill[i, j] = Fill[i + N, j + M] + 0.00000001;
                                            Y_N = true;

                                        }
                                        if (DEM[i, j] >= Fill[i + N, j + M] + 0.00000001)
                                        {
                                            Fill[i, j] = DEM[i, j];
                                        }
                                    }
                                }
                            }
                        }
                    }
                }

            } while (Y_N == true);
             
            return Fill;
        }


        //D8算法
        public int[,] D8(double[,] Fill)
        {
            /*流向矩阵表示如下：
             6  7  8
             5  K  1
             4  3  2    */

            //Fill[,] 是经过无洼地处理后的带有高程值的矩阵,  Direction[,]是流向矩阵（用于最后输出）
            int[,] Direction = new int[NRows, NCols];  //水流方向矩阵
            //Direction、Water初始化为0
            for (int i = 0; i < NRows; i++)//初始化为0
            {
                for (int j = 0; j < NCols; j++)
                {
                    Direction[i, j] = 0;
                }
            }

            //解求Direction矩阵
            double L = Convert.ToDouble(this.textBox3.Text);

            for (int i = 1; i < NRows - 1; i++)
            {
                for (int j = 1; j < NCols - 1; j++)
                {

                    double[] direction = new double[8];    //寄存高差                                             
                                                           // 用高程差在List中的顺序表示数字指代的方向
                    direction[0] = (Fill[i, j] - Fill[i, j + 1]) / L; direction[1] = (Fill[i, j] - Fill[i + 1, j + 1]) / (L * Math.Sqrt(2)); //1、2方向
                    direction[2] = (Fill[i, j] - Fill[i + 1, j]) / L; direction[3] = (Fill[i, j] - Fill[i + 1, j - 1]) / (L * Math.Sqrt(2));//3、4方向
                    direction[4] = (Fill[i, j] - Fill[i, j - 1]) / L; direction[5] = (Fill[i, j] - Fill[i - 1, j - 1]) / (L * Math.Sqrt(2)); //5、6方向
                    direction[6] = (Fill[i, j] - Fill[i - 1, j]) / L; direction[7] = (Fill[i, j] - Fill[i - 1, j + 1]) / (L * Math.Sqrt(2));//7、8方向
                    double max = direction[0];
                    int NUM = 0;
                    for (int n = 7; n >= 0; n--)
                    {
                        if (direction[n] >= max)
                        {
                            max = direction[n];
                            NUM = n + 1;
                        }
                    }
                    Direction[i, j] = NUM;

                }
            }
            return Direction;


        }


        //水量汇聚算法
        public int[,] Water(int[,] Direction)
        {
            int[,] WATER = new int[NRows, NCols];
            for (int i = 0; i < NRows; i++)
            {
                for (int j = 0; j < NCols; j++)
                {
                    WATER[i, j] = 0;
                }
            }
            for (int i = 1; i < NRows - 1; i++)
            {
                for (int j = 1; j < NCols - 1; j++)
                {
                    WATER[i, j] = Cal(i, j, Direction, WATER);

                }
            }
            int Cal(int i, int j, int[,] dir, int[,] water)
            {

                if (water[i, j] == 0)
                {
                    water[i, j] = 1;
                    //中间点的情况,8方向
                    if ((i < NRows - 1) && (i > 0) && (j < NCols - 1) && (j > 0))
                    {
                        if (dir[i, j + 1] == 5)
                        { water[i, j] += Cal(i, j + 1, dir, water); }
                        if (dir[i + 1, j + 1] == 6)
                        { water[i, j] += Cal(i + 1, j + 1, dir, water); }
                        if (dir[i + 1, j - 1] == 8)
                        { water[i, j] += Cal(i + 1, j - 1, dir, water); }
                        if (dir[i, j - 1] == 1)
                        { water[i, j] += Cal(i, j - 1, dir, water); }
                        if (dir[i - 1, j - 1] == 2)
                        { water[i, j] += Cal(i - 1, j - 1, dir, water); }
                        if (dir[i - 1, j] == 3)
                        { water[i, j] += Cal(i - 1, j, dir, water); }
                        if (dir[i - 1, j + 1] == 4)
                        { water[i, j] += Cal(i - 1, j + 1, dir, water); }
                        if (dir[i + 1, j] == 7)
                        { water[i, j] += Cal(i + 1, j, dir, water); }
                    }
                }
                return water[i, j];
            }
            return WATER;
        }

完整代码：
https://download.csdn.net/download/qq_45767140/20323976

欢迎大家批评指正。