大地坐标系相互转换、测线规划和缓冲区生成源码介绍

大地坐标系相互转换源码介绍

基于一个小的需求,写的一个大地坐标系相互转换小demo,其中包括WGS84,北京54,西安80,和国家2000之间的相互转换,并支持UTM投影和高斯投影。

源码介绍:

1、数据读取

//加载点数据
        private void AddPointData_Click(object sender, EventArgs e)
        {
            OpenFileDialog p_dlg = new OpenFileDialog();
            p_dlg.Filter = "All files(*.*)|*.*|txt files (*.txt)|*.txt";
            p_dlg.FilterIndex = 2;
            p_dlg.RestoreDirectory = true;
            if (p_dlg.ShowDialog() == DialogResult.OK)
            {
                if (m_time != 0)
                {
                    m_time = 0;
                    m_GPSPoint54.Clear();
                    m_GPSPoint80.Clear();
                    m_GPSPoint84.Clear();

                    m_GPSBL54.Clear();
                    m_GPSBL80.Clear();
                    m_GPSBL84.Clear();
                }
                //读取文件类型
                if(radioButton4.Checked == true || radioButton7.Checked == true)
                {
                    OpenDataxyH(p_dlg.FileName);
                }
                else if (radioButton8.Checked == true || radioButton12.Checked == true)
                {
                    OpenDataBLH(p_dlg.FileName);
                }
                
            }
            
            
        }

读取投影坐标

private void OpenDataxyH(string p_filename)
        {
            //System.Windows.Forms.MessageBox.Show(m_Openfilename);
            m_time += 1; //打开文件次数加1
            StreamReader p_FRead = File.OpenText(p_filename);
            int p_num = 0;
            bool p_est = true;
            string p_s = p_FRead.ReadLine();
            while (p_est)
            {
                p_s = String_Replace(p_s, ' ');
                string[] p_ss = p_s.Split(' ');
                byte[] p_label  = new byte[10];
                p_label = System.Text.Encoding.Default.GetBytes(p_ss[0]);
                if (p_ss[0] == "")
                {

                }
                else if(p_label[0] =='P')
                {
                    break;
                }
                //if (p_ss[0] == "PT001")
                //{
                //    break;
                //}
                p_s = p_FRead.ReadLine();
            }
            //System.Windows.Forms.MessageBox.Show(p_s);
            while (p_s != null)
            {
                p_s = String_Replace(p_s, ' ');
                string[] p_ss = p_s.Split(' ');
                // string[] p_ss = System.Text.RegularExpressions.Regex.Split(p_s,@[*]+);
                double[] p_apoint = new double[2];
                p_apoint[0] = double.Parse(p_ss[1]);
                p_apoint[1] = double.Parse(p_ss[2]);

                GpsInfo p_PointGPS = new GpsInfo();
                p_PointGPS.Longitude = p_apoint[0];
                p_PointGPS.Latitude = p_apoint[1];

                if(radioButton1.Checked == true)
                {
                     m_GPSPoint54.Add(p_PointGPS);
                }
                else if(radioButton2.Checked == true)
                {
                    m_GPSPoint80.Add(p_PointGPS);
                }
                else if (radioButton3.Checked == true)
                {
                    m_GPSPoint84.Add(p_PointGPS);
                }
                else if (radioButton13.Checked == true)
                {
                    m_GPSPoint2000.Add(p_PointGPS);
                }
               
                p_num++;
                p_s = p_FRead.ReadLine();
                //System.Windows.Forms.MessageBox.Show(p_PointGPS.Latitude.ToString());
            }
            // 
        }

读取经纬度数据

 private void OpenDataBLH(string p_filename)
        { 
             m_time += 1; //打开文件次数加1
            StreamReader p_FRead = File.OpenText(p_filename);
            int p_num = 0;
            bool p_est = true;
            string p_s = p_FRead.ReadLine();
            while (p_est)
            {
                p_s = String_Replace(p_s, ' ');
                string[] p_ss = p_s.Split(' ');
                byte[] p_label  = new byte[10];
                p_label = System.Text.Encoding.Default.GetBytes(p_ss[0]);
                if (p_ss[0] == "")
                {

                }
                else if(p_label[0] =='L')
                {
                    p_s = p_FRead.ReadLine();
                    break;
                }
                //if (p_ss[0] == "PT001")
                //{
                //    break;
                //}
                p_s = p_FRead.ReadLine();
            }
            //System.Windows.Forms.MessageBox.Show(p_s);
            while (p_s != null)
            {
                p_s = String_Replace(p_s, ' ');
                string[] p_ss = p_s.Split(' ');
                // string[] p_ss = System.Text.RegularExpressions.Regex.Split(p_s,@[*]+);
                double[] p_apoint = new double[2];
                p_apoint[0] = double.Parse(p_ss[0]);
                p_apoint[1] = double.Parse(p_ss[1]);

                GpsInfo p_PointGPS = new GpsInfo();
                p_PointGPS.Longitude = p_apoint[0];
                p_PointGPS.Latitude = p_apoint[1];

                if (radioButton1.Checked == true)
                {
                    m_GPSBL54.Add(p_PointGPS);
                }
                else if (radioButton2.Checked == true)
                {
                    m_GPSBL80.Add(p_PointGPS);
                }
                else if (radioButton3.Checked == true)
                {
                    m_GPSBL84.Add(p_PointGPS);
                }
                else if (radioButton13.Checked == true)
                {
                    m_GPSBL2000.Add(p_PointGPS);
                }

                p_num++;
                p_s = p_FRead.ReadLine();
                //System.Windows.Forms.MessageBox.Show(p_PointGPS.Latitude.ToString());
            }
        }


2、坐标转换

坐标转换类CoordTran

各坐标系常量(大地坐标参数)

 public double m_54a = 6378245.0;
        public double m_54b = 6356863.0188;
        public double m_54c = 1.0 / 298.3;

        public double m_80a = 6378140;
        public double m_80b = 6356755.2881575;
        public double m_80c = 1.0 / 298.25722101;

        public double m_84a = 6378137.0;
        public double m_84b = 6356752.314245;
        public double m_84c = 1.0 / 298.257223565;

        public double m_2000a = 6378137.0;
        public double m_2000b = 6356752.314140;
        public double m_2000c = 1.0/298.257222101;

1)大地坐标与空间直角坐标系相互转换

 //大地坐标系换换成空间直角坐标系
        public void BLHToXYZ(double B, double L, double H, double aAxis, double bAxis)
        {
            double dblD2R = Math.PI / 180;
            double e1 = Math.Sqrt(Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / aAxis;
            
            double N = aAxis / Math.Sqrt(1.0 - Math.Pow(e1, 2) * Math.Pow(Math.Sin(B * dblD2R), 2));
            targetX = (N + H) * Math.Cos(B * dblD2R) * Math.Cos(L * dblD2R);
            targetY = (N + H) * Math.Cos(B * dblD2R) * Math.Sin(L * dblD2R);
            targetZ = (N * (1.0 - Math.Pow(e1, 2)) + H) * Math.Sin(B * dblD2R);

        }
 //七参数变换
        public void TranParameter(transParaSeven p_sevenPara)
        {
            OriX = p_sevenPara.m_dbXRotate / 180 * Math.PI;
            OriY = p_sevenPara.m_dbYRotate / 180 * Math.PI;
            OriZ = p_sevenPara.m_dbZRotate / 180 * Math.PI;
            NewX = p_sevenPara.m_dbXOffset + p_sevenPara.m_dbScale * targetX - targetY * OriZ + targetZ * OriY + targetX;
            NewY = p_sevenPara.m_dbYOffset + p_sevenPara.m_dbScale * targetY + targetX * OriZ - targetZ * OriX + targetY;
            NewZ = p_sevenPara.m_dbZOffset + p_sevenPara.m_dbScale * targetZ - targetX * OriY + targetY * OriX + targetZ;

        }

 //空间直接坐标系转换为大地坐标系
        public void XYZtoBLH(double X, double Y, double Z, double aAxis, double bAxis)
        { 
            double e1 = (Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(aAxis, 2);
            double e2 = (Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(bAxis, 2);
            double S = Math.Sqrt(Math.Pow(X, 2) + Math.Pow(Y, 2));
            double cosL = X / S;
            double B = 0;
            double L = 0;
            //System.Windows.Forms.MessageBox.Show(e1.ToString());
            L = Math.Acos(cosL);
            L = Math.Abs(L);

            double tanB = Z / S;
            B = Math.Atan(tanB);
            double c = aAxis * aAxis / bAxis;
            double preB0 = 0.0;
            double ll = 0.0;
            double N = 0.0;
            //迭代计算纬度
            do {
                preB0 = B;
                ll = Math.Pow(Math.Cos(B), 2) * e2;
                N = c / Math.Sqrt(1 + ll);

                tanB = (Z + N * e1 * Math.Sin(B)) / S;
                B = Math.Atan(tanB);
            }
            while (Math.Abs(preB0 - B) >= 0.0000000001);

            ll = Math.Pow(Math.Cos(B), 2) * e2;
            N = c / Math.Sqrt(1 + ll);

            targetZ = Z / Math.Sin(B) - N * (1 - e1);
            targetB = B * 180 / Math.PI;
            targetL = L * 180 / Math.PI;

         }

2)高斯投影坐标与大地坐标想换转换
//高斯投影逆变换
        public void GaussProjInvCal(double X, double Y, double a , double f,out double longitude, out double latitude)
        {
            int ProjNo; int ZoneWide; ////带宽   
            double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
            double e1, e2, ee, NN, T, C, M, D, R, u, fai, iPI;
            iPI = 0.0174532925199433; ////3.1415926535898/180.0;   
            //a = 6378245.0; f = 1.0 / 298.3; //54年北京坐标系参数   
            ////a=6378140.0; f=1/298.257; //80年西安坐标系参数   
            ZoneWide = 6; ////6度带宽   
            ProjNo = (int)(X / 1000000L); //查找带号  
            longitude0 = (ProjNo - 1) * ZoneWide + ZoneWide / 2;
            longitude0 = longitude0 * iPI; //中央经线  
            X0 = ProjNo * 1000000L + 500000L;
            Y0 = 0;
            xval = X - X0; yval = Y - Y0; //带内大地坐标  
            e2 = 2 * f - f * f;
            e1 = (1.0 - Math.Sqrt(1 - e2)) / (1.0 + Math.Sqrt(1 - e2));
            ee = e2 / (1 - e2);
            M = yval;
            u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
            fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.Sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.Sin(4 * u)
             + (151 * e1 * e1 * e1 / 96) * Math.Sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.Sin(8 * u);
            C = ee * Math.Cos(fai) * Math.Cos(fai);
            T = Math.Tan(fai) * Math.Tan(fai);
            NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(fai) * Math.Sin(fai)); //字串1
            R = a * (1 - e2) / Math.Sqrt((1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin
                     (fai) * Math.Sin(fai)));
            D = xval / NN;
            //计算经度(Longitude) 纬度(Latitude)  
            longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D
             * D * D * D * D / 120) / Math.Cos(fai);
            latitude1 = fai - (NN * Math.Tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
             + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
            //转换为度 DD  
            longitude = longitude1 / iPI;
            latitude = latitude1 / iPI;  

        }

 //高斯投影由经纬度(UnitD)反算大地坐标(含带号,Unit:Metres) 
        public void GaussProjCal(double longitude, double latitude, double a , double f,ref double X, ref double Y) 
        {
            int ProjNo=0; int ZoneWide; ////带宽 
            double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
            double e2,ee, NN, T,C,A, M, iPI;
            iPI = 0.0174532925199433; ////3.1415926535898/180.0; 
            ZoneWide = 6; ////6度带宽 
            //a=6378245.0; f=1.0/298.3; //54年北京坐标系参数 
            ////a=6378140.0; f=1/298.257; //80年西安坐标系参数 
            ProjNo = (int)(longitude / ZoneWide) ; 
            longitude0 = ProjNo * ZoneWide + ZoneWide / 2; 
            longitude0 = longitude0 * iPI ;
            latitude0=0; 
            longitude1 = longitude * iPI ; //经度转换为弧度
            latitude1 = latitude * iPI ; //纬度转换为弧度
            e2=2*f-f*f;
            ee=e2*(1.0-e2);
            NN=a/Math.Sqrt(1.0-e2*Math.Sin(latitude1)*Math.Sin(latitude1));
            T = Math.Tan(latitude1) * Math.Tan(latitude1);
            C = ee * Math.Cos(latitude1) * Math.Cos(latitude1);
            A = (longitude1 - longitude0) * Math.Cos(latitude1);
            M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(2 * latitude1)
               + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.Sin(6 * latitude1));
            xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
            yval = M + NN * Math.Tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);
           
            // X0 = 1000000L * (ProjNo + 1) + 500000L;//加带号
            X0 =  500000L;
            Y0 = 0;
            xval = xval+X0; yval = yval+Y0; 
            X = xval;
            Y = yval;
         }


3)各坐标系分类型转换函数

//高斯投影平面坐标转大地经纬度坐标(不同坐标系下)
        public void xyHtoBLH(double x, double y, double H, double p_addL, transParaSeven p_ParaSeven, double p_xa, double p_xb, double p_xc,double p_Ba, double p_Bb)
        {
            //高斯投影逆变换
            GaussProjInvCal(x, y, p_xa,p_xc,out targetL, out targetB);
            targetL += p_addL;//加上带号

            BLHToXYZ(targetB, targetL, targetZ, p_xa, p_xb);
            TranParameter(p_ParaSeven);
            XYZtoBLH(NewX, NewY, NewZ, p_Ba, p_Bb);
           // targetL += p_addL;
            //直接调用targetL、targetB即可
        }
        //不同坐标系的高斯平面直角坐标互转
        public void xyHtoxyH(double x, double y, double H, double p_addL, transParaSeven p_ParaSeven, double p_xa, double p_xb, double p_xc,double p_Ba, double p_Bb,double p_Bc)
        {
            //高斯投影逆变换
            GaussProjInvCal(x, y, p_xa,p_xc,out targetL, out targetB);
            targetL += p_addL;//加上带号
           
            BLHToXYZ(targetB, targetL, targetZ, p_xa, p_xb);
            TranParameter(p_ParaSeven);
            XYZtoBLH(NewX, NewY, NewZ, p_Ba, p_Bb);
            GaussProjCal(targetL, targetB, p_Ba,p_Bc,ref targetx, ref targety);
            //直接调用targetx、targety即可
        }
        //不同经纬度坐标互转
        public void BLHtoBLH(double p_B, double p_L, double p_Z, double p_addL, transParaSeven p_ParaSeven, double p_xa, double p_xb, double p_Ba, double p_Bb)
        {
           // p_Z = 100;
            BLHToXYZ(p_B, p_L, p_Z, p_xa, p_xb);

            TranParameter(p_ParaSeven);
            XYZtoBLH(NewX, NewY, NewZ, p_Ba, p_Bb);
            //targetL += p_addL;
            //直接调用targetL、targetB即可
        }
        //大地经纬度坐标转高斯平面直角坐标
        public void BLHtoxyH(double p_B, double p_L, double p_Z, transParaSeven p_ParaSeven, double p_xa, double p_xb, double p_Ba, double p_Bb,double p_Bc)
        {
            BLHToXYZ(p_B, p_L, p_Z, p_xa, p_xb);
            TranParameter(p_ParaSeven);
            XYZtoBLH(NewX, NewY, NewZ, p_Ba, p_Bb);
            GaussProjCal(targetL, targetB, p_Ba, p_Bc,ref targetx, ref targety);
            //直接调用targetx、targety即可
        }

        //高斯投影转UTM投影
        public void GStoUTM(double x, double y, double H, out double O_x, out double O_y, out double O_H) 
        {
            O_x = TransformFactor * x;
            O_y = TransformFactor * y;
            O_H = H;
        }

        //UTM投影转高斯投影
        public void UTMtoGS(double x, double y, double H, out double O_x, out double O_y, out double O_H)
        {
            O_x = x / TransformFactor;
            O_y = y / TransformFactor;
            O_H = H;

        }

3、规划航迹

 //规划航迹
        private void PlanTrack_Click(object sender, EventArgs e)
        {
            List<double> p_Dis = new List<double>();
            List<double> p_CDis = new List<double>();
            CoordTran p_CoordTran = new CoordTran();
            PointProcess p_pproces = new PointProcess();

            BufferArea p_buffer = new BufferArea();
            double p_laL = double.Parse(textBox1.Text);
            double p_laR = -p_laL;
            //生成航迹
            if(radioButton9.Checked == true )
            {
                List<GpsInfo> m_GPSNew54_L = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew54_R = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew54 = new List<GpsInfo>();
                //保存距离
               
                p_pproces.DistanceCalculateBL(m_GPSBL54,p_Dis,p_CDis,p_CoordTran,p_CoordTran.m_54a,p_CoordTran.m_54c);
                SavePointDistance(p_CDis, p_Dis, "北京54坐标");

                //生成浅剖面
                p_buffer.CreateBufferArea2(m_GPSBL54, m_GPSNew54_L, p_laL, 0, p_CoordTran, m_proLine, p_CoordTran.m_54a, p_CoordTran.m_54c);               
                p_buffer.CreateBufferArea2(m_GPSBL54, m_GPSNew54_R, p_laR, 0, p_CoordTran, m_proLine, p_CoordTran.m_54a, p_CoordTran.m_54c);

                p_buffer.LineDissect(m_GPSNew54_L, m_GPSNew54_R, m_GPSNew54);

                string gpxName = System.Environment.CurrentDirectory + "\\" + "北京54"  + "浅剖"+p_laL.ToString()+".gpx";
                m_GpxWriter = new GpxWriter(gpxName);
                int p_num = m_GPSNew54.Count;
                for (int p_i = 0; p_i < p_num; p_i++)
                {
                    m_GpxWriter.AddGpsInfo(m_GPSNew54[p_i]);
                    //System.Windows.Forms.MessageBox.Show(m_GPSPoint84[p_i].Longitude.ToString());
                }
                m_GpxWriter.WriteToGpx();
            }
            else if(radioButton10.Checked == true)
            {
                List<GpsInfo> m_GPSNew80_L = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew80_R = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew80 = new List<GpsInfo>();
                
                p_pproces.DistanceCalculateBL(m_GPSBL80, p_Dis, p_CDis, p_CoordTran, p_CoordTran.m_80a, p_CoordTran.m_80c);
                SavePointDistance(p_CDis, p_Dis, "西安80坐标");

                //生成浅剖面
                p_buffer.CreateBufferArea2(m_GPSBL80, m_GPSNew80_L, p_laL, 0, p_CoordTran, m_proLine, p_CoordTran.m_80a, p_CoordTran.m_80c);
                p_buffer.CreateBufferArea2(m_GPSBL80, m_GPSNew80_R, p_laR, 0, p_CoordTran, m_proLine, p_CoordTran.m_80a, p_CoordTran.m_80c);

                p_buffer.LineDissect(m_GPSNew80_L, m_GPSNew80_R, m_GPSNew80);
                string gpxName = System.Environment.CurrentDirectory + "\\" + "西安80" + "浅剖" + p_laL.ToString() + ".gpx";
                m_GpxWriter = new GpxWriter(gpxName);
                int p_num = m_GPSNew80.Count;
                for (int p_i = 0; p_i < p_num; p_i++)
                {
                    m_GpxWriter.AddGpsInfo(m_GPSNew80[p_i]);
                    //System.Windows.Forms.MessageBox.Show(m_GPSPoint84[p_i].Longitude.ToString());
                }
                m_GpxWriter.WriteToGpx();
            }
            else if (radioButton11.Checked == true)
            {
                List<GpsInfo> m_GPSNew84_L = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew84_R = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew84 = new List<GpsInfo>();

                //int p_num = m_GPSBL84.Count;
                p_pproces.DistanceCalculateBL(m_GPSBL84, p_Dis, p_CDis, p_CoordTran, p_CoordTran.m_84a, p_CoordTran.m_84c);
                SavePointDistance(p_CDis, p_Dis, "WGS84坐标");

                //生成浅剖面
                p_buffer.CreateBufferArea2(m_GPSBL84, m_GPSNew84_L, p_laL, 0, p_CoordTran, m_proLine, p_CoordTran.m_84a, p_CoordTran.m_84c);
                p_buffer.CreateBufferArea2(m_GPSBL84, m_GPSNew84_R, p_laR, 0, p_CoordTran, m_proLine, p_CoordTran.m_84a, p_CoordTran.m_84c);

                p_buffer.LineDissect(m_GPSNew84_L, m_GPSNew84_R, m_GPSNew84);
                string gpxName = System.Environment.CurrentDirectory + "\\" + "WGS84" + "浅剖" + p_laL.ToString() + ".gpx";
                m_GpxWriter = new GpxWriter(gpxName);
                int p_num = m_GPSNew84.Count;
                for (int p_i = 0; p_i < p_num; p_i++)
                {
                    m_GpxWriter.AddGpsInfo(m_GPSNew84[p_i]);
                    //System.Windows.Forms.MessageBox.Show(m_GPSPoint84[p_i].Longitude.ToString());
                }
                m_GpxWriter.WriteToGpx();
            }
            else if (radioButton16.Checked == true)
            {
                List<GpsInfo> m_GPSNew2000_L = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew2000_R = new List<GpsInfo>();
                List<GpsInfo> m_GPSNew2000 = new List<GpsInfo>();

                //int p_num = m_GPSBL84.Count;
                p_pproces.DistanceCalculateBL(m_GPSBL2000, p_Dis, p_CDis, p_CoordTran, p_CoordTran.m_2000a, p_CoordTran.m_2000c);
                SavePointDistance(p_CDis, p_Dis, "CGCS2000坐标");

                //生成浅剖面
                p_buffer.CreateBufferArea2(m_GPSBL2000, m_GPSNew2000_L, p_laL, 0, p_CoordTran, m_proLine, p_CoordTran.m_84a, p_CoordTran.m_84c);
                p_buffer.CreateBufferArea2(m_GPSBL2000, m_GPSNew2000_R, p_laR, 0, p_CoordTran, m_proLine, p_CoordTran.m_84a, p_CoordTran.m_84c);

                p_buffer.LineDissect(m_GPSNew2000_L, m_GPSNew2000_R, m_GPSNew2000);
                string gpxName = System.Environment.CurrentDirectory + "\\" + "WGS84" + "浅剖" + p_laL.ToString() + ".gpx";
                m_GpxWriter = new GpxWriter(gpxName);
                int p_num = m_GPSNew2000.Count;
                for (int p_i = 0; p_i < p_num; p_i++)
                {
                    m_GpxWriter.AddGpsInfo(m_GPSNew2000[p_i]);
                    //System.Windows.Forms.MessageBox.Show(m_GPSPoint84[p_i].Longitude.ToString());
                }
                m_GpxWriter.WriteToGpx();
            }
            

        }

4、航迹缓冲区生成

根据航迹点生成其对应的缓冲区(可以应用在航迹规划缓冲区,道路缓冲区等一系列线性点的缓冲区生成)

类:BufferArea

public double bufferWidth;       //缓冲区宽度
        public double bufferHeight;      //缓冲区的高度
        public double LineOrizontal;     //缓冲区为水平方向
        public double LineVertical;      //缓冲区为垂直方向

        //List<GpsInfo>
        //点转换
        public void PointTran(GpsInfo p_OldGPS, CoordTran p_CoordTran, double p_xa, double p_xc,ref double p_X,ref double p_Y)
        {
            double p_B, p_L, p_H;
            p_L = p_OldGPS.Longitude;
            p_B = p_OldGPS.Latitude;
            p_H = 0.0;
            p_CoordTran.GaussProjCal(p_L, p_B, p_xa, p_xc, ref p_X, ref p_Y);//转换
        }
        //点转换
        public void PointTranBL(double p_X, double p_Y, CoordTran p_CoordTran, double p_AddL, double p_xa, double p_xc, ref double p_B, ref double p_L)
        {            
            double p_Z = 0.0;
            double p_lon = 0.0;
            double p_lat = 0.0;
            p_CoordTran.GaussProjInvCal(p_X, p_Y, p_xa, p_xc, out p_lon, out p_lat);//转换
            p_L = p_lon + p_AddL;
            p_B = p_lat;
        }
        public void CreateBufferArea(List<GpsInfo> p_OldGPS, List<GpsInfo> p_NewGPS, double p_width, int p_lable, CoordTran p_CoordTran, double p_AddL, double p_xa, double p_xc)
        {
            int p_pointNum = p_OldGPS.Count;
            double point1_x = 0.0;
            double point1_y = 0.0;
            double point2_x = 0.0;
            double point2_y = 0.0;
            PointTran(p_OldGPS[0], p_CoordTran, p_xa, p_xc, ref point1_x, ref point1_y);
            PointTran(p_OldGPS[1], p_CoordTran, p_xa, p_xc, ref point2_x, ref point2_y);

            double p_a1 = (point1_y - point2_y) / (point1_x - point2_x);
            double p_a2 ,p_b2;
            if(p_width == 0.0)
            {
                for (int p_i = 0; p_i < p_pointNum;p_i++ )
                {
                    p_NewGPS.Add(p_OldGPS[p_i]);
                }
            }
            else
            {
                if (p_a1 == 0)   //生成垂直方向缓冲区
                {
                    //暂不处理
                }
                else
                {
                    p_a2 = -1.0 / p_a1;   //求直线斜率
                    p_b2 = point1_y - (p_a2 * point1_x);    //求直线偏移量

                    //采用直接计算方法
                    double p_d = Math.Pow(p_width, 2) /(Math.Pow(p_a2, 2) + 1.0);
                    double p_X1 = point1_x + Math.Sqrt(p_d);
                    double p_X2 = point1_x - Math.Sqrt(p_d);
                    double p_Y1 = p_a2 * p_X1 + p_b2;
                    double p_Y2 = p_a2 * p_X2 + p_b2;
                    double p_B1 = 0.0;
                    double p_L1 = 0.0;
                    double p_B2 = 0.0;
                    double p_L2 = 0.0; 
                    PointTranBL(p_X1, p_Y1, p_CoordTran, p_AddL, p_xa, p_xc, ref p_B1, ref p_L1);
                    PointTranBL(p_X2, p_Y2, p_CoordTran, p_AddL, p_xa, p_xc, ref p_B2, ref p_L2);

                    if (p_lable == 0)
                    {
                        for (int p_i = 0; p_i < p_pointNum; p_i++)
                        {
                            GpsInfo p_pinfo = new GpsInfo();
                            if (p_i == 0) //第一个点
                            {
                                if(p_width > 0 )
                                {
                                    p_pinfo.Latitude = p_B1;
                                    p_pinfo.Longitude = p_L1;
                                    p_pinfo.Altitude = p_i + 1;
                                    p_pinfo.Speed = "unistrong:104";
                                    p_NewGPS.Add(p_pinfo);
                                }
                                else
                                {
                                    p_pinfo.Latitude = p_B2;
                                    p_pinfo.Longitude = p_L2;
                                    p_pinfo.Altitude = p_i + 1;
                                    p_pinfo.Speed = "unistrong:104";
                                    p_NewGPS.Add(p_pinfo);
                                }
                               

                            }//第一个点结束//
                            else
                            {
                                double p_lon = p_OldGPS[p_i].Longitude - p_OldGPS[p_i - 1].Longitude;
                                double p_lat = p_OldGPS[p_i].Latitude - p_OldGPS[p_i - 1].Latitude;
                                p_pinfo.Longitude = p_NewGPS[p_i - 1].Longitude + p_lon;
                                p_pinfo.Latitude = p_NewGPS[p_i - 1].Latitude + p_lat;
                                p_pinfo.Altitude = p_i + 1;
                                p_pinfo.Speed = "unistrong:104";
                                p_NewGPS.Add(p_pinfo);
                            }
                        }//for结束
                    }
                    else if(p_lable == 1)
                    {
                        for (int p_i = 0; p_i < p_pointNum; p_i++)
                        {
                            GpsInfo p_pinfo = new GpsInfo();
                            if (p_i == 0) //第一个点
                            {
                                p_pinfo.Latitude = p_Y2;
                                p_pinfo.Longitude = p_X2;
                                p_pinfo.Altitude = p_i + 1;
                                p_pinfo.Speed = "unistrong:104";
                                p_NewGPS.Add(p_pinfo);

                            }//第一个点结束//
                            else
                            {
                                double p_lon = p_OldGPS[p_i].Longitude - p_OldGPS[p_i - 1].Longitude;
                                double p_lat = p_OldGPS[p_i].Latitude - p_OldGPS[p_i - 1].Latitude;
                                p_pinfo.Longitude = p_NewGPS[p_i - 1].Longitude + p_lon;
                                p_pinfo.Latitude = p_NewGPS[p_i - 1].Latitude + p_lat;
                                p_pinfo.Altitude = p_i + 1;
                                p_pinfo.Speed = "unistrong:104";
                                p_NewGPS.Add(p_pinfo);
                            }
                        }//for结束
                    }//if结束

                }//缓冲区判断
            }
        }
        public void CreateBufferArea2(List<GpsInfo> p_OldGPS, List<GpsInfo> p_NewGPS, double p_width, int p_lable, CoordTran p_CoordTran, double p_AddL, double p_xa, double p_xc)
        {
            int p_pointNum = p_OldGPS.Count;
            double point1_x = 0.0;
            double point1_y = 0.0;
            double point2_x = 0.0;
            double point2_y = 0.0;

            double p_a, p_b, p_d;
            p_a = p_b = p_d = 0.0;
            for (int p_ii = 0; p_ii < p_pointNum ;p_ii++ )
            {
                if (p_ii < p_pointNum - 1)
                {
                    PointTran(p_OldGPS[p_ii], p_CoordTran, p_xa, p_xc, ref point1_x, ref point1_y);
                    PointTran(p_OldGPS[p_ii + 1], p_CoordTran, p_xa, p_xc, ref point2_x, ref point2_y);

                     double p_a1 = (point1_y - point2_y) / (point1_x - point2_x);
                     double p_a2, p_b2, p_d2;//记录直线信息
                     p_a2 = p_b2 = p_d2 = 0.0;
                     if (p_width == 0.0)
                     {
                         for (int p_i = 0; p_i < p_pointNum; p_i++)
                         {
                             p_NewGPS.Add(p_OldGPS[p_i]);
                         }
                     }
                     else
                     {
                         if (p_a1 == 0)   //生成垂直方向缓冲区
                         {
                             //暂不处理
                         }
                         else
                         {
                             p_a2 = -1.0 / p_a1;   //求直线斜率
                             p_b2 = point1_y - (p_a2 * point1_x);    //求直线偏移量

                             //采用直接计算方法
                             p_d2 = Math.Pow(p_width, 2) / (Math.Pow(p_a2, 2) + 1.0);
                             double p_X1 = point1_x + Math.Sqrt(p_d2);
                             double p_X2 = point1_x - Math.Sqrt(p_d2);
                             double p_Y1 = p_a2 * p_X1 + p_b2;
                             double p_Y2 = p_a2 * p_X2 + p_b2;
                             double p_B1 = 0.0;
                             double p_L1 = 0.0;
                             double p_B2 = 0.0;
                             double p_L2 = 0.0;
                             PointTranBL(p_X1, p_Y1, p_CoordTran, p_AddL, p_xa, p_xc, ref p_B1, ref p_L1);
                             PointTranBL(p_X2, p_Y2, p_CoordTran, p_AddL, p_xa, p_xc, ref p_B2, ref p_L2);
                             GpsInfo p_pinfo = new GpsInfo();
                             if ((p_width > 0 && p_a2 > 0) || (p_width < 0 && p_a2 < 0))
                             {
                                 p_pinfo.Latitude = p_B1;
                                 p_pinfo.Longitude = p_L1;
                                 p_pinfo.Altitude = p_ii + 1;
                                 p_pinfo.Speed = "unistrong:104";
                                 p_NewGPS.Add(p_pinfo);
                             }
                             else
                             {
                                 p_pinfo.Latitude = p_B2;
                                 p_pinfo.Longitude = p_L2;
                                 p_pinfo.Altitude = p_ii + 1;
                                 p_pinfo.Speed = "unistrong:104";
                                 p_NewGPS.Add(p_pinfo);
                             }
                             

                         }
                     }
                     p_a = p_a2;
                     p_b = p_b2;
                     p_d = p_d2; 
                }
                else if(p_ii == p_pointNum -1)  //最后一个点
                {
                    //PointTran(p_OldGPS[p_ii], p_CoordTran, p_xa, p_xc, ref point1_x, ref point1_y);                   

                    //double p_X1 = point1_x + Math.Sqrt(p_d);
                    //double p_X2 = point1_x - Math.Sqrt(p_d);
                    //double p_Y1 = p_a * p_X1 + p_b;
                    //double p_Y2 = p_a * p_X2 + p_b;
                    //double p_B1 = 0.0;
                    //double p_L1 = 0.0;
                    //double p_B2 = 0.0;
                    //double p_L2 = 0.0;
                    //PointTranBL(p_X1, p_Y1, p_CoordTran, p_AddL, p_xa, p_xc, ref p_B1, ref p_L1);
                    //PointTranBL(p_X2, p_Y2, p_CoordTran, p_AddL, p_xa, p_xc, ref p_B2, ref p_L2);
                    GpsInfo p_pinfo = new GpsInfo();
                    double p_lon = p_OldGPS[p_ii].Longitude - p_OldGPS[p_ii - 1].Longitude;
                    double p_lat = p_OldGPS[p_ii].Latitude - p_OldGPS[p_ii - 1].Latitude;
                    p_pinfo.Longitude = p_NewGPS[p_ii - 1].Longitude + p_lon;
                    p_pinfo.Latitude = p_NewGPS[p_ii - 1].Latitude + p_lat;
                    p_pinfo.Altitude = p_ii + 1;
                    p_pinfo.Speed = "unistrong:104";
                    p_NewGPS.Add(p_pinfo);
                   
                    //if (p_width > 0)
                    //{
                    //    p_pinfo.Latitude = p_B1;
                    //    p_pinfo.Longitude = p_L1;
                    //    p_pinfo.Altitude = p_ii + 1;
                    //    p_pinfo.Speed = "unistrong:104";
                    //    p_NewGPS.Add(p_pinfo);
                    //}
                    //else
                    //{
                    //    p_pinfo.Latitude = p_B2;
                    //    p_pinfo.Longitude = p_L2;
                    //    p_pinfo.Altitude = p_ii + 1;
                    //    p_pinfo.Speed = "unistrong:104";
                    //    p_NewGPS.Add(p_pinfo);
                    //}
                             
                }
               


            }
        }
        //生成浅剖/
        public void LineDissect(List<GpsInfo> p_LGPS, List<GpsInfo> p_RGPS, List<GpsInfo> p_NewGPS)
        {
            int p_num1 = p_LGPS.Count;
            int p_num2 = p_RGPS.Count;
            int p_Label, p_Label_R, p_i, p_j;
            p_Label = p_Label_R = 0;
            int p_PointNum = 0; // 
            for (p_i = 0; p_i < p_num1; p_i++)
            {
                p_PointNum += 1; //点数+1              
                p_NewGPS.Add(p_LGPS[p_i]);
                p_NewGPS[p_PointNum - 1].Altitude = p_PointNum;
                //加载右侧的点
                while (p_Label < 2 && p_Label_R < p_num2)
                {
                    p_NewGPS.Add(p_RGPS[p_Label_R]);
                    p_Label_R += 1;
                    p_PointNum++; //点数+1
                    p_NewGPS[p_PointNum - 1].Altitude = p_PointNum;
                    p_Label++;
                }
                p_Label = 0; //初始化
                //加载左侧的点
                p_i += 1;
                if (p_i < p_num1)
                {
                    p_NewGPS.Add(p_LGPS[p_i]);
                    p_PointNum++; //点数+1
                    p_NewGPS[p_PointNum - 1].Altitude = p_PointNum;
                }


            }

缓冲区生成

//缓冲区生成
        private void BufferArea_Click(object sender, EventArgs e)
        {
            CoordTran p_CoordTran = new CoordTran();
            p_CoordTran.targetX = 0.0;
            p_CoordTran.targetY = 0.0;
            p_CoordTran.targetZ = 0.0;

            double p_interval = double.Parse(textBox12.Text); //间隔大小
            double p_low = double.Parse(textBox10.Text);
            double p_high = double.Parse(textBox11.Text); //数值上下限 
            BufferArea p_buffer = new BufferArea();
            int p_N =(int)((p_high - p_low) / p_interval);

            if (radioButton9.Checked == true)
            {
                
                for (int p_i = 0; p_i <= p_N; p_i++)
                {
                    List<GpsInfo> m_GPSNew54 = new List<GpsInfo>();
                    double p_la = p_low + p_i * p_interval;
                    string gpxName = System.Environment.CurrentDirectory + "\\" + "北京54" + p_la.ToString() + "m缓冲.gpx";
                    m_GpxWriter = new GpxWriter(gpxName);
                    int p_num = m_GPSBL54.Count;
                    p_buffer.CreateBufferArea2(m_GPSBL54, m_GPSNew54, p_la, 0, p_CoordTran, m_proLine, p_CoordTran.m_54a, p_CoordTran.m_54c);
                    for (int p_ii = 0; p_ii < p_num; p_ii++)
                    {
                        m_GpxWriter.AddGpsInfo(m_GPSNew54[p_ii]);
                        //System.Windows.Forms.MessageBox.Show(m_GPSPoint84[p_i].Longitude.ToString());
                    }
                    m_GpxWriter.WriteToGpx();

                }
            }
            else if (radioButton10.Checked == true)
            {
                
                for (int p_i = 0; p_i <= p_N; p_i++)
                {
                    List<GpsInfo> m_GPSNew80 = new List<GpsInfo>();
                    double p_la = p_low + p_i * p_interval;
                    string gpxName = System.Environment.CurrentDirectory + "\\" + "西安80" + p_la.ToString() + "m缓冲.gpx";
                    m_GpxWriter = new GpxWriter(gpxName);

                    int p_num = m_GPSBL80.Count;
                    p_buffer.CreateBufferArea2(m_GPSBL80, m_GPSNew80, p_la, 0, p_CoordTran, m_proLine, p_CoordTran.m_80a, p_CoordTran.m_80c);
                    for (int p_ii = 0; p_ii < p_num; p_ii++)
                    {
                        m_GpxWriter.AddGpsInfo(m_GPSNew80[p_ii]);
                        //System.Windows.Forms.MessageBox.Show(m_GPSPoint84[p_i].Longitude.ToString());
                    }
                    m_GpxWriter.WriteToGpx();
                }
            }
            else if (radioButton11.Checked == true)
            {
               
                for (int p_i = 0; p_i <= p_N; p_i++)
                {
                    List<GpsInfo> m_GPSNew84 = new List<GpsInfo>();
                    double p_la = p_low + p_i * p_interval;
                    string gpxName = System.Environment.CurrentDirectory + "\\" + "WGS84" + p_la.ToString() + "m缓冲.gpx";
                    m_GpxWriter = new GpxWriter(gpxName);
                    int p_num = m_GPSBL84.Count;
                    p_buffer.CreateBufferArea2(m_GPSBL84, m_GPSNew84, p_la, 0, p_CoordTran, m_proLine, p_CoordTran.m_84a, p_CoordTran.m_84c);
                    for (int p_ii = 0; p_ii < p_num; p_ii++)
                    {
                        m_GpxWriter.AddGpsInfo(m_GPSNew84[p_ii]);
                        //System.Windows.Forms.MessageBox.Show(m_GPSPoint84[p_i].Longitude.ToString());
                    }
                    m_GpxWriter.WriteToGpx();

                }
            }
            
           
        }

源码位置:http://download.csdn.net/download/u011326478/10246916

猜你喜欢

转载自blog.csdn.net/u011326478/article/details/79296161