source: subversion/applications/lib/ccoord/LatLng.h @ 18255

Last change on this file since 18255 was 2097, checked in by nick, 13 years ago

updated libosm

File size: 4.8 KB
Line 
1#ifndef LATLNG_H
2#define LATLNG_H
3
4#include "RefEll.h"
5#include "util.h"
6#include <cmath>
7
8// Based on Jonathan Stott's JCoord.
9// Licenced under the GNU GPL.
10
11class OSRef;
12
13/**
14 * Class to represent a latitude/longitude pair.
15 *
16 * (c) 2006 Jonathan Stott
17 *
18 * Created on 11-02-2006
19 *
20 * @author Jonathan Stott
21 * @version 1.0
22 * @since 1.0
23 */
24class LatLng {
25
26private:
27  /**
28   * Latitude in degrees.
29   */
30  double lat;
31
32  /**
33   * Longitude in degrees.
34   */
35  double lng;
36
37public:
38  /**
39   * Create a new LatLng object to represent a latitude/longitude pair.
40   *
41   * @param lat
42   *          the latitude in degrees
43   * @param lng
44   *          the longitude in degrees
45   * @since 1.0
46   */
47 
48   LatLng(double lat, double lng) {
49    this->lat = lat;
50    this->lng = lng;
51  }
52
53
54  /**
55   * Convert this->latitude and longitude into an OSGB (Ordnance Survey of Great
56   * Britain) grid reference.
57   *
58   * @return the converted OSGB grid reference
59   * @since 1.0
60   */
61  OSRef toOSRef(); 
62
63
64  /**
65   * Convert this->LatLng from the OSGB36 datum to the WGS84 datum using an
66   * approximate Helmert transformation.
67   *
68   * @since 1.0
69   */
70  void toWGS84() {
71         RefEll AIRY_1830 (6377563.396, 6356256.909);
72    double a = AIRY_1830.getMaj();
73    double eSquared = AIRY_1830.getEcc();
74    double phi = toRadians(lat);
75    double lambda = toRadians(lng);
76    double v = a / (sqrt(1 - eSquared * sinSquared(phi)));
77    double H = 0; // height
78    double x = (v + H) * cos(phi) * cos(lambda);
79    double y = (v + H) * cos(phi) * sin(lambda);
80    double z = ((1 - eSquared) * v + H) * sin(phi);
81
82    double tx = 446.448;
83    double ty = -124.157;
84    double tz = 542.060;
85    double s = -0.0000204894;
86    double rx = toRadians(0.00004172222);
87    double ry = toRadians(0.00006861111);
88    double rz = toRadians(0.00023391666);
89
90    double xB = tx + (x * (1 + s)) + (-rx * y) + (ry * z);
91    double yB = ty + (rz * x) + (y * (1 + s)) + (-rx * z);
92    double zB = tz + (-ry * x) + (rx * y) + (z * (1 + s));
93
94         RefEll WGS84     (6378137.000, 6356752.3141);
95    a = WGS84.getMaj();
96    eSquared = WGS84.getEcc();
97
98    double lambdaB = toDegrees(atan(yB / xB));
99    double p = sqrt((xB * xB) + (yB * yB));
100    double phiN = atan(zB / (p * (1 - eSquared)));
101    for (int i = 1; i < 10; i++) {
102      v = a / (sqrt(1 - eSquared * sinSquared(phiN)));
103      double phiN1 = atan((zB + (eSquared * v * sin(phiN))) / p);
104      phiN = phiN1;
105    }
106
107    double phiB = toDegrees(phiN);
108
109    lat = phiB;
110    lng = lambdaB;
111  }
112
113
114  /**
115   * Convert this->LatLng from the WGS84 datum to the OSGB36 datum using an
116   * approximate Helmert transformation.
117   *
118   * @since 1.0
119   */
120  void toOSGB36() {
121    RefEll wgs84 (6378137.000, 6356752.3141);
122    double a = wgs84.getMaj();
123    double eSquared = wgs84.getEcc();
124    double phi = toRadians(lat);
125    double lambda = toRadians(lng);
126    double v = a / (sqrt(1 - eSquared * sinSquared(phi)));
127    double H = 0; // height
128    double x = (v + H) * cos(phi) * cos(lambda);
129    double y = (v + H) * cos(phi) * sin(lambda);
130    double z = ((1 - eSquared) * v + H) * sin(phi);
131
132    double tx = -446.448;
133    double ty = 124.157;
134    double tz = -542.060;
135    double s = 0.0000204894;
136    double rx = toRadians(-0.00004172222);
137    double ry = toRadians(-0.00006861111);
138    double rz = toRadians(-0.00023391666);
139
140    double xB = tx + (x * (1 + s)) + (-rx * y) + (ry * z);
141    double yB = ty + (rz * x) + (y * (1 + s)) + (-rx * z);
142    double zB = tz + (-ry * x) + (rx * y) + (z * (1 + s));
143
144    RefEll airy1830 (6377563.396, 6356256.909);
145    a = airy1830.getMaj();
146    eSquared = airy1830.getEcc();
147
148    double lambdaB = toDegrees(atan(yB / xB));
149    double p = sqrt((xB * xB) + (yB * yB));
150    double phiN = atan(zB / (p * (1 - eSquared)));
151    for (int i = 1; i < 10; i++) {
152      v = a / (sqrt(1 - eSquared * sinSquared(phiN)));
153      double phiN1 = atan((zB + (eSquared * v * sin(phiN))) / p);
154      phiN = phiN1;
155    }
156
157    double phiB = toDegrees(phiN);
158
159    lat = phiB;
160    lng = lambdaB;
161  }
162
163
164  /**
165   * Calculate the surface distance in kilometres from the this->LatLng to the
166   * given LatLng.
167   *
168   * @param ll
169   * @return the surface distance in km
170   * @since 1.0
171   */
172  double distance(LatLng ll) {
173    double er = 6366.707;
174
175    double latFrom = toRadians(getLat());
176    double latTo = toRadians(ll.getLat());
177    double lngFrom = toRadians(getLng());
178    double lngTo = toRadians(ll.getLng());
179
180    double d =
181        acos(sin(latFrom) * sin(latTo) + cos(latFrom)
182            * cos(latTo) * cos(lngTo - lngFrom))
183            * er;
184
185    return d;
186  }
187
188
189  /**
190   * Return the latitude in degrees.
191   *
192   * @return the latitude in degrees
193   * @since 1.0
194   */
195  double getLat() {
196    return lat;
197  }
198
199
200  /**
201   * Return the longitude in degrees.
202   *
203   * @return the longitude in degrees
204   * @since 1.0
205   */
206  double getLng() {
207    return lng;
208  }
209};
210
211#endif
Note: See TracBrowser for help on using the repository browser.