git.maemo.org Git - ptas/blob - src/coordinatesystemtransformer.cpp

   1 #include "coordinatesystemtransformer.h"
   2
   3 #include <math.h>
   4
   5
   6 using namespace QTM_NAMESPACE;
   7
   8 /**
   9  * A class representing geographical coordinates from the Finnish KKJ coordinate system.
  10  */
  11 class KKJGeoCoordinate {
  12     // The latitude of the coordinate
  13     double mLatitude;
  14     // The longitude of the coordinate
  15     double mLongitude;
  16
  17 public:
  18     /**
  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
  22      * on earth.
  23      */
  24     KKJGeoCoordinate();
  25
  26     /**
  27      * Constructs a new KKJ geographical coordinate.
  28      * @param latitude the latitude of the coordinate
  29      * @param longitude the longitude of the coordinate
  30      */
  31     KKJGeoCoordinate(double latitude, double longitude);
  32
  33     /**
  34      * Gets the latitude of this coordinate.
  35      * @return the latitude of this coordinate
  36      */
  37     double latitude() const;
  38
  39     /**
  40      * Gets the longitude of this coordinate.
  41      * @return the longitude of this coordinate
  42      */
  43     double longitude() const;
  44
  45     /**
  46      * Sets the latitude of this coordinate.
  47      * @param latitude the new latitude
  48      */
  49     void setLatitude(double latitude);
  50
  51     /**
  52      * Sets the longitude of this coordinate.
  53      * @param longitude the new latitude
  54      */
  55     void setLongitude(double longitude);
  56
  57 };
  58
  59
  60 KKJGeoCoordinate::KKJGeoCoordinate() :
  61         mLatitude(0.0),
  62         mLongitude(0.0)
  63 {
  64 }
  65
  66 KKJGeoCoordinate::KKJGeoCoordinate(double latitude, double longitude) :
  67         mLatitude(latitude),
  68         mLongitude(longitude)
  69 {
  70 }
  71
  72 double KKJGeoCoordinate::latitude() const
  73 {
  74     return mLatitude;
  75 }
  76
  77 double KKJGeoCoordinate::longitude() const
  78 {
  79     return mLongitude;
  80 }
  81
  82 void KKJGeoCoordinate::setLatitude(double latitude)
  83 {
  84     mLatitude = latitude;
  85 }
  86
  87 void KKJGeoCoordinate::setLongitude(double longitude)
  88 {
  89     mLongitude = longitude;
  90 }
  91
  92
  93 /**
  94  * A Hayford reference ellipsoid
  95  */
  96 class HayfordEllipsoid {
  97 public:
  98     // Equatorial radius
  99     static const double a;
 100     // Flattening
 101     static const double f;
 102     // Polar radius
 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;
 110     // Second flattening
 111     static const double n;
 112     // Second flattening squared
 113     static const double nn;
 114 };
 115
 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;
 124
 125
 126
 127 // Degrees to radians
 128 double radians(double deg) {
 129     return deg * M_PI / 180.0;
 130 }
 131
 132 // Radians to degrees
 133 double degrees(double rad) {
 134     return rad * 180.0 / M_PI;
 135 }
 136
 137
 138 /**
 139  * A class providing KKJ zone information.
 140  */
 141 class KKJZone {
 142 private:
 143     // Storage for the central meridians and false eastings of the zones
 144     static const double KKJ_ZONE_INFO[6][2];
 145
 146     // Minimum zone number
 147     static const int MIN_ZONE = 0;
 148     // Maximum zone number
 149     static const int MAX_ZONE = 5;
 150
 151 public:
 152     /**
 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
 157      */
 158     static int getZoneNumberFromEasting(const KKJGridCoordinate &kkj);
 159
 160     /**
 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
 165      */
 166     static int getZoneNumberFromLongitude(const KKJGeoCoordinate &kkj);
 167
 168     /**
 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
 173      */
 174     static double getCentralMeridianOfZone(int zoneNumber);
 175
 176     /**
 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
 181      */
 182     static double getFalseEastingOfZone(int zoneNumber);
 183 };
 184
 185                                         // central meridian
 186                                                   // false easting
 187 const double KKJZone::KKJ_ZONE_INFO[6][2] = { {18.0,  500000.0},  // zone 0
 188                                               {21.0, 1500000.0},
 189                                               {24.0, 2500000.0},
 190                                               {27.0, 3500000.0},
 191                                               {30.0, 4500000.0},
 192                                               {33.0, 5500000.0} };// zone 5
 193
 194
 195 int KKJZone::getZoneNumberFromEasting(const KKJGridCoordinate &kkj) {
 196     int zoneNumber = floor(kkj.easting() / 1000000.0);
 197     if (zoneNumber < MIN_ZONE || zoneNumber > MAX_ZONE) {
 198         zoneNumber = -1;
 199     }
 200
 201     return zoneNumber;
 202 }
 203
 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) {
 212             break;
 213         }
 214         zoneNumber--;
 215     }
 216
 217     return zoneNumber;
 218 }
 219
 220 double KKJZone::getCentralMeridianOfZone(int zoneNumber)
 221 {
 222     if (zoneNumber >= MIN_ZONE && zoneNumber <= MAX_ZONE) {
 223         return KKJ_ZONE_INFO[zoneNumber][0];
 224     }
 225
 226     return 0.0;
 227 }
 228
 229 double KKJZone::getFalseEastingOfZone(int zoneNumber)
 230 {
 231     if (zoneNumber >= MIN_ZONE && zoneNumber <= MAX_ZONE) {
 232         return KKJ_ZONE_INFO[zoneNumber][1];
 233     }
 234
 235     return 0.0;
 236 }
 237
 238
 239 /**
 240  * Transforms a KKJ geographical coordinate to a WGS84 geographical coordinate.
 241  * @param fromCoordinate the input coordinate
 242  * @return the transformed coordinate
 243  */
 244 QGeoCoordinate transformToWGS84GeoCoordinate(const KKJGeoCoordinate &fromCoordinate) {
 245     double kkjla = fromCoordinate.latitude();
 246     double kkjlo = fromCoordinate.longitude();
 247
 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;
 250
 251     return QGeoCoordinate(kkjla + dLa, kkjlo + dLo);
 252 }
 253
 254
 255 /**
 256  * Transforms a WGS84 geographical coordinate to a KKJ geographical coordinate.
 257  * @param fromCoordinate the input coordinate
 258  * @return the transformed coordinate
 259  */
 260 KKJGeoCoordinate transformToKKJGeoCoordinate(const QGeoCoordinate &fromCoordinate) {
 261     double longitude = fromCoordinate.longitude();
 262     double latitude = fromCoordinate.latitude();
 263
 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;
 266
 267     return KKJGeoCoordinate(latitude + dLa, longitude + dLo);
 268 }
 269
 270
 271 /**
 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
 276  */
 277 KKJGridCoordinate 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;
 290
 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 KKJGridCoordinate(outY, outX);
 294 }
 295
 296 /**
 297  * Transforms a KKJ rectangular grid coordinate to a KKJ geographical coordinate.
 298  * @param fromCoordinate the input coordinate
 299  * @return the transformed coordinate
 300  */
 301 KKJGeoCoordinate transformToKKJGeoCoordinate(const KKJGridCoordinate &fromCoordinate) {
 302     // Scan iteratively the target area, until find matching
 303     // KKJ coordinate value. Area is defined with Hayford Ellipsoid.
 304     double minLo = 18.5;
 305     double maxLo = 32.0;
 306     double minLa = 59.0;
 307     double maxLa = 70.5;
 308
 309     int i = 1;
 310
 311     KKJGeoCoordinate ret;
 312
 313     while (i < 35) {
 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         KKJGridCoordinate kkj = transformToKKJGridCoordinate(ret);
 319         if (kkj.northing() < fromCoordinate.northing()) {
 320             minLa = minLa + 0.45 * deltaLa;
 321         } else {
 322             maxLa = minLa + 0.55 * deltaLa;
 323         }
 324
 325         if (kkj.easting() < fromCoordinate.easting()) {
 326             minLo = minLo + 0.45 * deltaLo;
 327         } else {
 328             maxLo = minLo + 0.55 * deltaLo;
 329         }
 330
 331         i++;
 332     }
 333
 334     return ret;
 335 }
 336
 337
 338 KKJGridCoordinate CoordinateSystemTransformer::transformToKKJ(const QGeoCoordinate &fromCoordinate)
 339 {
 340     KKJGeoCoordinate tmpKKJ = transformToKKJGeoCoordinate(fromCoordinate);
 341     return transformToKKJGridCoordinate(tmpKKJ);
 342 }
 343
 344 QGeoCoordinate CoordinateSystemTransformer::transformToWGS84(const KKJGridCoordinate &fromCoordinate)
 345 {
 346     KKJGeoCoordinate tmpKKJ = transformToKKJGeoCoordinate(fromCoordinate);
 347     return transformToWGS84GeoCoordinate(tmpKKJ);
 348 }