洪水淹没分析

背景:

    在使用skyline开发的过程中,有些功能并不满足实际需求,所以我们就需要自己重写设计以及做一些计算,比如模拟水淹过程,水淹洼地分析,以及水淹面积S和体积V计算等等.

需求整个计算过程:

    水淹算法本身并不难,只是种子填充算法的扩展. 通俗点讲,就是像素之间的连通域问题;  难的是在从鼠标选择区域->切割改区域的地形数据->将三维地形数据降维处理->应用算法处理之后->反映射原始数据->

三维数据面片化->构建棱柱体->计算共面非共面的不规则体积以及面积; 最后才能得到自己想要的功能以及实现,哦,对了,这个结果最后还要可视化;

  大家看这个过程,算法本身只是占了一小部分.剩下的用程序展现这个东西才是重点所在,大部分难题都是工程问题.以及开发过程中的种种异常以及流畅体验(异步多线程以及多进程等问题)

先展示一个结果,然后具体聊实现过程

这个结果我和其它主流的分析软件计算的结果对比过,相差不大,只不过其它的软件在计算水淹高度的过程中,可能设置了固定的采样间隔,我没有设置,

所以,可能计算的结果更逼近一些.(虽然这个意义感觉并不大而已)

实现:

数据获取我就不写了,可以用GDAL获取一个tif图像,转三维数据;或者直接从地球上裁剪;具体视情况而异

1.拿到数据之后,我们要将三维数据栅格化.就是

   public bool ElevationData2BitmP(ref List<ElevationPoint3D> ListFaces, int width, int height, double waterHeight, out Bitmap a)
        {
            a = new Bitmap(width, height);
            if (width < 1 && height < 1)
            {
                return false;
            }

            System.Drawing.Rectangle rect = new System.Drawing.Rectangle(0, 0, a.Width, a.Height);
            System.Drawing.Imaging.BitmapData bmpData = a.LockBits(rect, System.Drawing.Imaging.ImageLockMode.ReadWrite,
                                System.Drawing.Imaging.PixelFormat.Format24bppRgb);            //填充数据.
            int stride = bmpData.Stride;

            int Listnumber = ListFaces.Count;
            //查询的步长.
            int querySride = width * height / Listnumber;

            double lowWaterHeightNumber = 0;
            unsafe
            {
                byte* pIn = (byte*)bmpData.Scan0.ToPointer();
                byte* P;
                int R, G, B;
                double temp = 0;
                Random rm = new Random();
                for (int y = 0, k = 0; y < a.Height; y++)
                {
                    for (int x = 0; x < a.Width; x++, k += querySride)
                    {
                       // double dscal = 0;
                        P = pIn;
                        B = P[0];
                        G = P[1];
                        R = P[2];

                        if (k < Listnumber)
                        {
                            temp = ListFaces[k].PxyzConvertImg.Z;
                            if (temp >= waterHeight)  //背景色，以及过滤色
                            {

                                //dscal = 0;
                                P[2] = (byte)(73);//R
                                P[1] = (byte)(221);//G
                                P[0] = (byte)(93);//B
                            }
                            else
                            {
                                //dscal = 1;
                                P[2] = (byte)(0);//R
                                P[1] = (byte)(0);//G
                                P[0] = (byte)(255);//B
                                lowWaterHeightNumber++;
                            }
                        }

                        //P[2] = (byte)(255 * dscal);//R
                        //P[1] = (byte)(255 * dscal);//G
                        //P[0] = (byte)(255 * dscal);//B
                        pIn += 3;
                    }

                    pIn += stride - a.Width * 3;
                }

            }

            a.UnlockBits(bmpData);
            //得到当前蓝色点.
            return true;
        }

这里解释一下: 这一部分主要是将三维数据映射成二维uv数据,方便后续计算

通过这一步,就可以实现其数据的可视化. 这一步完成之后,就可以设置一不同高度,显示不同的计算结果

2.计算完成之后,将三维数据网格化

这样可以分析计算每一个小面片的具体参数; 注意其需要依据水淹高度, 在底层构建一个虚假的面.和这一层构成棱柱体.

类似这样.

3.计算,主要是利用海伦公式快速计算

 /// <summary>
        /// 依据传入的面,以及高. 求取拟柱体体积.
        /// </summary> 点的绘制按顺时针计算.
        /// <returns></returns>
        public double CalculateFitting(ref QuadPrism quadranglePrism)
        {
            double cylinderVolume = 0;

            //以点出发,计算每个点到水平面的距离.
            // double currentPlength

            //1/6*(S0 + 4S1 + S2)*H
            //注意这里的面积,一般都按照梯形面接计算.
            //点1 quadrangleFace.item0,  item1 ,item2 , item3.

            //  从L->         R->               RB->            LB 
            //  item0,        item1             item2           item3
            //  Distance[0]   Distance[1]       Distance[2]    Distance[3]
            Point3 currentP = quadranglePrism._bottomFace.pxyz0;
            Point3 RightcurrentP = quadranglePrism._bottomFace.pxyz1;
            Point3 RightBottomcurrentP = quadranglePrism._bottomFace.pxyz2;
            Point3 LeftBottomcurrentP = quadranglePrism._bottomFace.pxyz3;

            //SO面积.
            double S0h = ComputorFunction.calculateD2Distance(ref currentP, ref LeftBottomcurrentP);
            double S0Area = (quadranglePrism.PointToWaterLevelDistance[0] + quadranglePrism.PointToWaterLevelDistance[3]) / 2.0 * S0h;

            //S2面积.
            double S2h = ComputorFunction.calculateD2Distance(ref RightcurrentP, ref RightBottomcurrentP);
            double S2Area = (quadranglePrism.PointToWaterLevelDistance[1] + quadranglePrism.PointToWaterLevelDistance[2]) / 2.0 * S2h;

            //S1面积,按照中值划分,可以知道
            double S1_upBorder = 0.5 * (quadranglePrism.PointToWaterLevelDistance[0] + quadranglePrism.PointToWaterLevelDistance[1]);
            double S1_BottomBorder = 0.5 * (quadranglePrism.PointToWaterLevelDistance[2] + quadranglePrism.PointToWaterLevelDistance[3]);
            double S1h = ComputorFunction.calculateD2Distance(ref RightcurrentP, ref RightBottomcurrentP);
            double S1Area = (S1_upBorder + S1_BottomBorder) * S1h * 0.5;

            //整个体积的高.
            double Sh = ComputorFunction.calculateD3Distance(ref currentP, ref RightcurrentP);
            //幸普斯规则.数值积分进行拟柱体体积计算.
            cylinderVolume = (S0Area + 4 * S1Area + S2Area) / 6.0 * Sh;

            return cylinderVolume;
        }

4.最后使用循环体,累加计算就可以了

5.这里总结一下:写的过程中,一定要细致,否则特变容易出一些问题