ghostkcleung

睇住條 Formula，
一知半解代返啲符號入去，
終於叫做轉到，應該係有啲 Roundup 上嘅誤差，
2,3 米範圍內，不過無所謂，反正 GPS 本身都自己都誤差。

測試結果：

lat: 22.48577143777777820
lon: 114.00029998111111240
n: 838478.464
e: 818098.341

1. #include <math.h>
   
2. #include <iostream>
   
3. #include <iomanip>
   
4. 
   
5. using namespace std;
   
6. 
   
7. double dms_to_dd(int deg, int min, double sec)
   
8. {
   
9.     return deg + (double)min / 60 + sec / 3600;
   
10. }
    
11. 
    
12. template <typename T>
    
13. double deg2rad(T deg)
    
14. {
    
15.     return (double)deg * M_PI / 180.0;
    
16. }
    
17. 
    
18. /*
    
19. The Formula is referenced from the following PDF document:
    
20. https://www.geodetic.gov.hk/common/data/pdf/explanatorynotes_c.pdf
    
21. */
    
22. 
    
23. const double e2 = 0.006722670022;
    
24. const double A0 = 1 - e2 / 4 - 3 * e2 * e2 / 64;
    
25. const double A2 = 3 * (e2 + e2 * e2 / 4) / 8;
    
26. const double A4 = 15 * e2 * e2 / 256;
    
27. const double a = 6378388;
    
28. const double N0 = 819069.80;
    
29. const double M0 = 2468395.728;
    
30. const double Vs = 6381480.502;
    
31. const double lon0 = deg2rad(dms_to_dd(114, 10, 42.80));
    
32. const double E0 = 836694.05;
    
33. const double iso_lat = 1.005792635;
    
34. 
    
35. void wgs84_to_hk1980(double lat, double lon, double &x, double &y)
    
36. {
    
37.     lat += 5.5 / 3600;
    
38.     lon -= 8.8 / 3600;
    
39. 
    
40.     lat = deg2rad(lat);
    
41.     lon = deg2rad(lon);
    
42. 
    
43.     double M = a * (A0 * lat - A2 * sin(2 * lat) + A4 * sin(4 * lat));
    
44.     double lonP = lon - lon0;
    
45. 
    
46.     x = N0 + M - M0 + Vs * sin(lat) * lonP * lonP / 2 * cos(lat);
    
47.     y = E0 + Vs * lonP * cos(lat) + Vs * pow(lonP, 3) / 6 * pow(cos(lat), 3) * (iso_lat - pow(tan(lat), 2));
    
48. }
    
49. 
    
50. int main()
    
51. {
    
52. 
    
53.     /*
    
54.     This example is referenced from the following PDF document:
    
55.     https://www.geodetic.gov.hk/common/data/parameter/7P_ITRF96_HK80_V1(Tra_Chi).pdf
    
56. 
    
57.     Lat: 22d, 29m, 8.777176s
    
58.     Lon: 114d, 00m, 1.079932s
    
59.     N: 838477.970m
    
60.     E: 818097.267m
    
61.     */
    
62. 
    
63.     double lat = dms_to_dd(22, 29, 8.777176);
    
64.     double lon = dms_to_dd(114, 00, 1.079932);
    
65. 
    
66.     double x, y;
    
67. 
    
68.     wgs84_to_hk1980(lat, lon, x, y);
    
69. 
    
70.     cout << fixed << setprecision(17);
    
71.     cout << "lat: " << lat << endl;
    
72.     cout << "lon: " << lon << endl;
    
73. 
    
74.     cout << fixed << setprecision(3);
    
75.     cout << "n: " << x << endl;
    
76.     cout << "e: " << y << endl;
    
77. 
    
78.     return 0;
    
79. }

複製代碼