1 #include "coordinatesystemtransformer.h"
6 using namespace QTM_NAMESPACE;
9 * A class representing geographical coordinates from the Finnish KKJ coordinate system.
11 class KKJGeoCoordinate {
12 // The latitude of the coordinate
14 // The longitude of the coordinate
19 * An empty constructor.
20 * Constructs a new KKJ geographical coordinate with default coordinate values.
21 * Note that the constructed coordinate doesn't necessarily represent any real location
27 * Constructs a new KKJ geographical coordinate.
28 * @param latitude the latitude of the coordinate
29 * @param longitude the longitude of the coordinate
31 KKJGeoCoordinate(double latitude, double longitude);
34 * Gets the latitude of this coordinate.
35 * @return the latitude of this coordinate
37 double latitude() const;
40 * Gets the longitude of this coordinate.
41 * @return the longitude of this coordinate
43 double longitude() const;
46 * Sets the latitude of this coordinate.
47 * @param latitude the new latitude
49 void setLatitude(double latitude);
52 * Sets the longitude of this coordinate.
53 * @param longitude the new latitude
55 void setLongitude(double longitude);
60 KKJGeoCoordinate::KKJGeoCoordinate() :
66 KKJGeoCoordinate::KKJGeoCoordinate(double latitude, double longitude) :
72 double KKJGeoCoordinate::latitude() const
77 double KKJGeoCoordinate::longitude() const
82 void KKJGeoCoordinate::setLatitude(double latitude)
87 void KKJGeoCoordinate::setLongitude(double longitude)
89 mLongitude = longitude;
94 * A Hayford reference ellipsoid
96 class HayfordEllipsoid {
99 static const double a;
101 static const double f;
103 static const double b;
104 // Polar radius squared
105 static const double bb;
106 // Polar radius of curvature
107 static const double c;
108 // First eccentricity squared
109 static const double ee;
111 static const double n;
112 // Second flattening squared
113 static const double nn;
116 const double HayfordEllipsoid::a = 6378388.0;
117 const double HayfordEllipsoid::f = 1.0 / 297.0;
118 const double HayfordEllipsoid::b = (1.0 - f) * a;
119 const double HayfordEllipsoid::bb = b * b;
120 const double HayfordEllipsoid::c = (a / b) * a;
121 const double HayfordEllipsoid::ee = (a * a - bb) / bb; // should probably be (a * a - bb) / (a * a)
122 const double HayfordEllipsoid::n = (a - b) / (a + b);
123 const double HayfordEllipsoid::nn = n * n;
127 // Degrees to radians
128 double radians(double deg) {
129 return deg * M_PI / 180.0;
132 // Radians to degrees
133 double degrees(double rad) {
134 return rad * 180.0 / M_PI;
139 * A class providing KKJ zone information.
143 // Storage for the central meridians and false eastings of the zones
144 static const double KKJ_ZONE_INFO[6][2];
146 // Minimum zone number
147 static const int MIN_ZONE = 0;
148 // Maximum zone number
149 static const int MAX_ZONE = 5;
153 * Determines a zone number from the KKJ easting value. If the easting is not within any of
154 * the zones, -1 is returned.
155 * @param kkj the KKJ coordinate
156 * @return the zone number or -1 on error
158 static int getZoneNumberFromEasting(const KKJ &kkj);
161 * Determines a zone number from the KKJ longitude value. If the longitude is not within any of
162 * the zones, -1 is returned.
163 * @param kkj the KKJ coordinate
164 * @return the zone number or -1 on error
166 static int getZoneNumberFromLongitude(const KKJGeoCoordinate &kkj);
169 * Gets the central meridian in degrees of the given zone. The zone number must be
170 * on the interval [0, 5]. If an invalid zone number is given, 0.0 is returned.
171 * @param zoneNumber the zone number
172 * @return the central meridian of the zone or 0.0 on error
174 static double getCentralMeridianOfZone(int zoneNumber);
177 * Gets the false easting in metres of the given zone. The zone number must be
178 * on the interval [0, 5]. If an invalid zone number is given, 0.0 is returned.
179 * @param zoneNumber the zone number
180 * @return the false easting of the zone or 0.0 on error
182 static double getFalseEastingOfZone(int zoneNumber);
187 const double KKJZone::KKJ_ZONE_INFO[6][2] = { {18.0, 500000.0}, // zone 0
192 {33.0, 5500000.0} };// zone 5
195 int KKJZone::getZoneNumberFromEasting(const KKJ &kkj) {
196 int zoneNumber = floor(kkj.easting() / 1000000.0);
197 if (zoneNumber < MIN_ZONE || zoneNumber > MAX_ZONE) {
204 int KKJZone::getZoneNumberFromLongitude(const KKJGeoCoordinate &kkj) {
205 // determine the zonenumber from KKJ easting
206 // takes KKJ zone which has center meridian
207 // longitude nearest (in math value) to
208 // the given KKJ longitude
209 int zoneNumber = MAX_ZONE;
210 while (zoneNumber >= MIN_ZONE) {
211 if (fabs(kkj.longitude() - KKJ_ZONE_INFO[zoneNumber][0]) <= 1.5) {
220 double KKJZone::getCentralMeridianOfZone(int zoneNumber)
222 if (zoneNumber >= MIN_ZONE && zoneNumber <= MAX_ZONE) {
223 return KKJ_ZONE_INFO[zoneNumber][0];
229 double KKJZone::getFalseEastingOfZone(int zoneNumber)
231 if (zoneNumber >= MIN_ZONE && zoneNumber <= MAX_ZONE) {
232 return KKJ_ZONE_INFO[zoneNumber][1];
240 * Transforms a KKJ geographical coordinate to a WGS84 geographical coordinate.
241 * @param fromCoordinate the input coordinate
242 * @return the transformed coordinate
244 QGeoCoordinate transformToWGS84GeoCoordinate(const KKJGeoCoordinate &fromCoordinate) {
245 double kkjla = fromCoordinate.latitude();
246 double kkjlo = fromCoordinate.longitude();
248 double dLa = (0.124867E+01 + -0.269982E+00 * kkjla + 0.191330E+00 * kkjlo + 0.356119E-02 * kkjla * kkjla + -0.122312E-02 * kkjla * kkjlo + -0.335514E-03 * kkjlo * kkjlo) / 3600.0;
249 double dLo = (-0.286111E+02 + 0.114183E+01 * kkjla + -0.581428E+00 * kkjlo + -0.152421E-01 * kkjla * kkjla + 0.118177E-01 * kkjla * kkjlo + 0.826646E-03 * kkjlo * kkjlo) / 3600.0;
251 return QGeoCoordinate(kkjla + dLa, kkjlo + dLo);
256 * Transforms a WGS84 geographical coordinate to a KKJ geographical coordinate.
257 * @param fromCoordinate the input coordinate
258 * @return the transformed coordinate
260 KKJGeoCoordinate transformToKKJGeoCoordinate(const QGeoCoordinate &fromCoordinate) {
261 double longitude = fromCoordinate.longitude();
262 double latitude = fromCoordinate.latitude();
264 double dLa = (-0.124766E+01 + 0.269941E+00 * latitude + -0.191342E+00 * longitude + -0.356086E-02 * latitude * latitude + 0.122353E-02 * latitude * longitude + 0.335456E-03 * longitude * longitude) / 3600.0;
265 double dLo = (0.286008E+02 + -0.114139E+01 * latitude + 0.581329E+00 * longitude + 0.152376E-01 * latitude * latitude + -0.118166E-01 * latitude * longitude + -0.826201E-03 * longitude * longitude) / 3600.0;
267 return KKJGeoCoordinate(latitude + dLa, longitude + dLo);
272 * Transforms a KKJ geographical coordinate to a KKJ rectangular grid coordinate.
273 * @param fromCoordinate the input coordinate
274 * @param zoneNumber the zone number in which the input coordinate resides
275 * @return the transformed coordinate
277 KKJ transformToKKJGridCoordinate(const KKJGeoCoordinate &fromCoordinate) {
278 int zoneNumber = KKJZone::getZoneNumberFromLongitude(fromCoordinate);
279 double Lo = radians(fromCoordinate.longitude()) - radians(KKJZone::getCentralMeridianOfZone(zoneNumber));
280 double cosLa = cos(radians(fromCoordinate.latitude()));
281 double NN = HayfordEllipsoid::ee * cosLa * cosLa;
282 double LaF = atan(tan(radians(fromCoordinate.latitude())) / cos(Lo * sqrt(1.0 + NN)));
283 double cosLaF = cos(LaF);
284 double t = (tan(Lo) * cosLaF) / sqrt(1.0 + HayfordEllipsoid::ee * cosLaF * cosLaF);
285 double A = HayfordEllipsoid::a / (1.0 + HayfordEllipsoid::n);
286 double A1 = A * (1.0 + HayfordEllipsoid::nn / 4.0 + HayfordEllipsoid::nn * HayfordEllipsoid::nn / 64.0);
287 double A2 = A * 1.5 * HayfordEllipsoid::n * (1.0 - HayfordEllipsoid::nn / 8.0);
288 double A3 = A * 0.9375 * HayfordEllipsoid::nn * (1.0 - HayfordEllipsoid::nn / 4.0);
289 double A4 = A * 35.0 / 48.0 * HayfordEllipsoid::nn * HayfordEllipsoid::n;
291 unsigned int outY = A1 * LaF - A2 * sin(2.0 * LaF) + A3 * sin(4.0 * LaF) - A4 * sin(6.0 * LaF);
292 unsigned int outX = HayfordEllipsoid::c * log(t + sqrt(1.0 + t * t)) + 500000.0 + zoneNumber * 1000000.0;
293 return KKJ(outY, outX);
297 * Transforms a KKJ rectangular grid coordinate to a KKJ geographical coordinate.
298 * @param fromCoordinate the input coordinate
299 * @return the transformed coordinate
301 KKJGeoCoordinate transformToKKJGeoCoordinate(const KKJ &fromCoordinate) {
302 // Scan iteratively the target area, until find matching
303 // KKJ coordinate value. Area is defined with Hayford Ellipsoid.
311 KKJGeoCoordinate ret;
314 double deltaLo = maxLo - minLo;
315 double deltaLa = maxLa - minLa;
316 ret.setLongitude(minLo + 0.5 * deltaLo);
317 ret.setLatitude(minLa + 0.5 * deltaLa);
318 KKJ kkj = transformToKKJGridCoordinate(ret);
319 if (kkj.northing() < fromCoordinate.northing()) {
320 minLa = minLa + 0.45 * deltaLa;
322 maxLa = minLa + 0.55 * deltaLa;
325 if (kkj.easting() < fromCoordinate.easting()) {
326 minLo = minLo + 0.45 * deltaLo;
328 maxLo = minLo + 0.55 * deltaLo;
338 KKJ CoordinateSystemTransformer::transformToKKJ(const QGeoCoordinate &fromCoordinate)
340 KKJGeoCoordinate tmpKKJ = transformToKKJGeoCoordinate(fromCoordinate);
341 return transformToKKJGridCoordinate(tmpKKJ);
344 QGeoCoordinate CoordinateSystemTransformer::transformToWGS84(const KKJ &fromCoordinate)
346 KKJGeoCoordinate tmpKKJ = transformToKKJGeoCoordinate(fromCoordinate);
347 return transformToWGS84GeoCoordinate(tmpKKJ);