兰伯特投影转经纬度算法

heywjiulh

求助，兰伯特投影转经纬度的算法？

luckycomcn

自己写的兰伯特投影反解

1.             /// <summary>
   
2.             /// 兰伯特投影转换为经纬投影，手动指定长半轴和短半轴
   
3.             /// </summary>
   
4.             /// <param name="LonRef">参考经度，为了方便，按照角度输入</param>
   
5.             /// <param name="LatRef">参考纬度，为了方便，按照角度输入</param>
   
6.             /// <param name="Slat">m=1时的南纬度，为了方便，按照角度输入</param>
   
7.             /// <param name="Nlat">m=1时的北纬度，为了方便，按照角度输入</param>
   
8.             /// <param name="deltaY">Y方向上的距离，以米为单位</param>
   
9.             /// <param name="deltaX">X方向上的距离，以米为单位</param>
   
10.             /// <param name="EarA">地球长半轴</param>
    
11.             /// <param name="EarB">地球短半轴</param>
    
12.             /// <returns></returns>
    
13.             public static double[] lccToLatLon(double LonRef, double LatRef, double Slat, double Nlat, double deltaY, double deltaX, double EarA, double EarB)
    
14.             {
    
15.                 //这个函数参考了European Petroleum Survey Group的Coordinate Conversions and Transformations including Formulas
    
16.                 double[] latAndLon = new double[2];
    
17.                 double L0, B0, B1, B2, Y, X;
    
18.                 double a, b, f, e;    //四个变量分别是地球长半轴，地球短半轴，扁率，第一偏心率
    
19.                 double mB1, mB2, tB0, tB1, tB2, n, F, r0, rd, td, Thetad, L, B, errf;
    
20.                 int num;
    
21.                 double pi = Math.PI; //PI，圆周率
    
22.                 a = EarA;    //地球长半轴
    
23.                 b = EarB;    //地球短半轴
    
24.                 L0 = LonRef;
    
25.                 B0 = LatRef;
    
26.                 B1 = Slat;
    
27.                 B2 = Nlat;
    
28.                 Y = deltaY;
    
29.                 X = deltaX;
    
30.                 //各种输入经纬度转化成弧度
    
31.                 L0 = L0 / 180 * pi;
    
32.                 B0 = B0 / 180 * pi;
    
33.                 B1 = B1 / 180 * pi;
    
34.                 B2 = B2 / 180 * pi;
    
35.                 f = (a - b) / a;
    
36.                 e = Math.Pow((2 * f - Math.Pow(f, 2)), 0.5);
    
37. 
    
38.                 mB1 = Math.Cos(B1) / Math.Pow((1 - Math.Pow(e, 2) * Math.Sin(B1) * Math.Sin(B1)), 0.5);
    
39.                 mB2 = Math.Cos(B2) / Math.Pow((1 - Math.Pow(e, 2) * Math.Sin(B2) * Math.Sin(B2)), 0.5);
    
40.                 tB0 = Math.Tan(pi / 4 - B0 / 2) / Math.Pow(((1 - e * Math.Sin(B0)) / (1 + e * Math.Sin(B0))), (e / 2));
    
41.                 tB1 = Math.Tan(pi / 4 - B1 / 2) / Math.Pow(((1 - e * Math.Sin(B1)) / (1 + e * Math.Sin(B1))), (e / 2));
    
42.                 tB2 = Math.Tan(pi / 4 - B2 / 2) / Math.Pow(((1 - e * Math.Sin(B2)) / (1 + e * Math.Sin(B2))), (e / 2));
    
43.                 n = Math.Log(mB1 / mB2) / Math.Log(tB1 / tB2);
    
44.                 F = mB1 / (n * Math.Pow(tB1, n));
    
45.                 r0 = a * F * Math.Pow(tB0, n);
    
46.                 rd = Math.Sign(n) * Math.Pow((Math.Pow(Y, 2) + Math.Pow((r0 - X), 2)), 0.5);
    
47.                 td = Math.Pow((rd / (a * F)), (1 / n));
    
48.                 Thetad = Math.Atan(Y / (r0 - X));
    
49.                 L = Thetad / n + L0;
    
50.                 B = B0;
    
51.                 errf = 1;
    
52.                 num = 0;
    
53. 
    
54.                 while (Math.Abs(errf) > (float)(1e-8) && num < 2000)
    
55.                 {
    
56.                     B1 = pi / 2 - 2 * Math.Atan(td * Math.Pow(((1 - e * Math.Sin(B)) / (1 + e * Math.Sin(B))), (e / 2)));
    
57.                     errf = B1 - B;
    
58.                     B = B1;
    
59.                     num = num + 1;
    
60.                 }
    
61.                 B = B * 180 / pi;
    
62.                 L = L * 180 / pi;
    
63.                 latAndLon[0] = L;
    
64.                 latAndLon[1] = B;
    
65.                 return latAndLon;
    
66.             }
    

            /// <summary><br />             /// 兰伯特投影转换为经纬投影，手动指定长半轴和短半轴<br />             /// </summary><br />             /// <param name="LonRef">参考经度，为了方便，按照角度输入</param><br />             /// <param name="LatRef">参考纬度，为了方便，按照角度输入</param><br />             /// <param name="Slat">m=1时的南纬度，为了方便，按照角度输入</param><br />             /// <param name="Nlat">m=1时的北纬度，为了方便，按照角度输入</param><br />             /// <param name="deltaY">Y方向上的距离，以米为单位</param><br />             /// <param name="deltaX">X方向上的距离，以米为单位</param><br />             /// <param name="EarA">地球长半轴</param><br />             /// <param name="EarB">地球短半轴</param><br />             /// <returns></returns><br />             public static double[] lccToLatLon(double LonRef, double LatRef, double Slat, double Nlat, double deltaY, double deltaX, double EarA, double EarB)<br />             {<br />                 //这个函数参考了European Petroleum Survey Group的Coordinate Conversions and Transformations including Formulas<br />                 double[] latAndLon = new double[2];<br />                 double L0, B0, B1, B2, Y, X;<br />                 double a, b, f, e;    //四个变量分别是地球长半轴，地球短半轴，扁率，第一偏心率<br />                 double mB1, mB2, tB0, tB1, tB2, n, F, r0, rd, td, Thetad, L, B, errf;<br />                 int num;<br />                 double pi = Math.PI; //PI，圆周率<br />                 a = EarA;    //地球长半轴<br />                 b = EarB;    //地球短半轴<br />                 L0 = LonRef;<br />                 B0 = LatRef;<br />                 B1 = Slat;<br />                 B2 = Nlat;<br />                 Y = deltaY;<br />                 X = deltaX;<br />                 //各种输入经纬度转化成弧度<br />                 L0 = L0 / 180 * pi;<br />                 B0 = B0 / 180 * pi;<br />                 B1 = B1 / 180 * pi;<br />                 B2 = B2 / 180 * pi;<br />                 f = (a - b) / a;<br />                 e = Math.Pow((2 * f - Math.Pow(f, 2)), 0.5);<br /> <br />                 mB1 = Math.Cos(B1) / Math.Pow((1 - Math.Pow(e, 2) * Math.Sin(B1) * Math.Sin(B1)), 0.5);<br />                 mB2 = Math.Cos(B2) / Math.Pow((1 - Math.Pow(e, 2) * Math.Sin(B2) * Math.Sin(B2)), 0.5);<br />                 tB0 = Math.Tan(pi / 4 - B0 / 2) / Math.Pow(((1 - e * Math.Sin(B0)) / (1 + e * Math.Sin(B0))), (e / 2));<br />                 tB1 = Math.Tan(pi / 4 - B1 / 2) / Math.Pow(((1 - e * Math.Sin(B1)) / (1 + e * Math.Sin(B1))), (e / 2));<br />                 tB2 = Math.Tan(pi / 4 - B2 / 2) / Math.Pow(((1 - e * Math.Sin(B2)) / (1 + e * Math.Sin(B2))), (e / 2));<br />                 n = Math.Log(mB1 / mB2) / Math.Log(tB1 / tB2);<br />                 F = mB1 / (n * Math.Pow(tB1, n));<br />                 r0 = a * F * Math.Pow(tB0, n);<br />                 rd = Math.Sign(n) * Math.Pow((Math.Pow(Y, 2) + Math.Pow((r0 - X), 2)), 0.5);<br />                 td = Math.Pow((rd / (a * F)), (1 / n));<br />                 Thetad = Math.Atan(Y / (r0 - X));<br />                 L = Thetad / n + L0;<br />                 B = B0;<br />                 errf = 1;<br />                 num = 0;<br /> <br />                 while (Math.Abs(errf) > (float)(1e-8) && num < 2000)<br />                 {<br />                     B1 = pi / 2 - 2 * Math.Atan(td * Math.Pow(((1 - e * Math.Sin(B)) / (1 + e * Math.Sin(B))), (e / 2)));<br />                     errf = B1 - B;<br />                     B = B1;<br />                     num = num + 1;<br />                 }<br />                 B = B * 180 / pi;<br />                 L = L * 180 / pi;<br />                 latAndLon[0] = L;<br />                 latAndLon[1] = B;<br />                 return latAndLon;<br />             }<br />

heywjiulh

兰伯特投影的反解公式，纬度B是如何实现迭代的

zqq0ew0

{:e311_b:}{:e311_b:}{:e311_b:}{:e311_b:}同求

ayzqs

luckycomcn 发表于 2012-12-10 20:14

                               
登录/注册后可看大图

自己写的兰伯特投影反解

谢谢分享！

黄河沙

很不错啊

kongfeng0824

luckycomcn 发表于 2012-12-10 20:14

                               
登录/注册后可看大图

自己写的兰伯特投影反解

感谢感谢！！！！

SEAMAN

学习学习，有问题再来请教

zitaxcy

可以请楼主分析下经纬度转兰伯特投影的公式和依据么..?

Liuhuamao

zitaxcy 发表于 2013-12-25 00:04

                               
登录/注册后可看大图

可以请楼主分析下经纬度转兰伯特投影的公式和依据么..?

学习中......,感谢分享