Coordinate system transformer

Makes transformations between WGS84 and KKJ rectangular grid coordinate systems

diff --git a/src/coordinatesystemtransformer.cpp b/src/coordinatesystemtransformer.cpp
new file mode 100644 (file)
index 0000000..1632133
--- /dev/null
+++ b/src/coordinatesystemtransformer.cpp
@@ -0,0 +1,348 @@
+#include "coordinatesystemtransformer.h"
+
+#include <math.h>
+
+
+using namespace QTM_NAMESPACE;
+
+/**
+ * A class representing geographical coordinates from the Finnish KKJ coordinate system.
+ */
+class KKJGeoCoordinate {
+    // The latitude of the coordinate
+    double mLatitude;
+    // The longitude of the coordinate
+    double mLongitude;
+
+public:
+    /**
+     * An empty constructor.
+     * Constructs a new KKJ geographical coordinate with default coordinate values.
+     * Note that the constructed coordinate doesn't necessarily represent any real location
+     * on earth.
+     */
+    KKJGeoCoordinate();
+
+    /**
+     * Constructs a new KKJ geographical coordinate.
+     * @param latitude the latitude of the coordinate
+     * @param longitude the longitude of the coordinate
+     */
+    KKJGeoCoordinate(double latitude, double longitude);
+
+    /**
+     * Gets the latitude of this coordinate.
+     * @return the latitude of this coordinate
+     */
+    double latitude() const;
+
+    /**
+     * Gets the longitude of this coordinate.
+     * @return the longitude of this coordinate
+     */
+    double longitude() const;
+
+    /**
+     * Sets the latitude of this coordinate.
+     * @param latitude the new latitude
+     */
+    void setLatitude(double latitude);
+
+    /**
+     * Sets the longitude of this coordinate.
+     * @param longitude the new latitude
+     */
+    void setLongitude(double longitude);
+
+};
+
+
+KKJGeoCoordinate::KKJGeoCoordinate() :
+        mLatitude(0.0),
+        mLongitude(0.0)
+{
+}
+
+KKJGeoCoordinate::KKJGeoCoordinate(double latitude, double longitude) :
+        mLatitude(latitude),
+        mLongitude(longitude)
+{
+}
+
+double KKJGeoCoordinate::latitude() const
+{
+    return mLatitude;
+}
+
+double KKJGeoCoordinate::longitude() const
+{
+    return mLongitude;
+}
+
+void KKJGeoCoordinate::setLatitude(double latitude)
+{
+    mLatitude = latitude;
+}
+
+void KKJGeoCoordinate::setLongitude(double longitude)
+{
+    mLongitude = longitude;
+}
+
+
+/**
+ * A Hayford reference ellipsoid
+ */
+class HayfordEllipsoid {
+public:
+    // Equatorial radius
+    static const double a;
+    // Flattening
+    static const double f;
+    // Polar radius
+    static const double b;
+    // Polar radius squared
+    static const double bb;
+    // Polar radius of curvature
+    static const double c;
+    // First eccentricity squared
+    static const double ee;
+    // Second flattening
+    static const double n;
+    // Second flattening squared
+    static const double nn;
+};
+
+const double HayfordEllipsoid::a = 6378388.0;
+const double HayfordEllipsoid::f  = 1.0 / 297.0;
+const double HayfordEllipsoid::b  = (1.0 - f) * a;
+const double HayfordEllipsoid::bb = b * b;
+const double HayfordEllipsoid::c  = (a / b) * a;
+const double HayfordEllipsoid::ee = (a * a - bb) / bb; // should probably be (a * a - bb) / (a * a)
+const double HayfordEllipsoid::n = (a - b) / (a + b);
+const double HayfordEllipsoid::nn = n * n;
+
+
+
+// Degrees to radians
+double radians(double deg) {
+    return deg * M_PI / 180.0;
+}
+
+// Radians to degrees
+double degrees(double rad) {
+    return rad * 180.0 / M_PI;
+}
+
+
+/**
+ * A class providing KKJ zone information.
+ */
+class KKJZone {
+private:
+    // Storage for the central meridians and false eastings of the zones
+    static const double KKJ_ZONE_INFO[6][2];
+
+    // Minimum zone number
+    static const int MIN_ZONE = 0;
+    // Maximum zone number
+    static const int MAX_ZONE = 5;
+
+public:
+    /**
+     * Determines a zone number from the KKJ easting value. If the easting is not within any of
+     * the zones, -1 is returned.
+     * @param kkj the KKJ coordinate
+     * @return the zone number or -1 on error
+     */
+    static int getZoneNumberFromEasting(const KKJ &kkj);
+
+    /**
+     * Determines a zone number from the KKJ longitude value. If the longitude is not within any of
+     * the zones, -1 is returned.
+     * @param kkj the KKJ coordinate
+     * @return the zone number or -1 on error
+     */
+    static int getZoneNumberFromLongitude(const KKJGeoCoordinate &kkj);
+
+    /**
+     * Gets the central meridian in degrees of the given zone. The zone number must be
+     * on the interval [0, 5]. If an invalid zone number is given, 0.0 is returned.
+     * @param zoneNumber the zone number
+     * @return the central meridian of the zone or 0.0 on error
+     */
+    static double getCentralMeridianOfZone(int zoneNumber);
+
+    /**
+     * Gets the false easting in metres of the given zone. The zone number must be
+     * on the interval [0, 5]. If an invalid zone number is given, 0.0 is returned.
+     * @param zoneNumber the zone number
+     * @return the false easting of the zone or 0.0 on error
+     */
+    static double getFalseEastingOfZone(int zoneNumber);
+};
+
+                                        // central meridian
+                                                  // false easting
+const double KKJZone::KKJ_ZONE_INFO[6][2] = { {18.0,  500000.0},  // zone 0
+                                              {21.0, 1500000.0},
+                                              {24.0, 2500000.0},
+                                              {27.0, 3500000.0},
+                                              {30.0, 4500000.0},
+                                              {33.0, 5500000.0} };// zone 5
+
+
+int KKJZone::getZoneNumberFromEasting(const KKJ &kkj) {
+    int zoneNumber = floor(kkj.easting() / 1000000.0);
+    if (zoneNumber < MIN_ZONE || zoneNumber > MAX_ZONE) {
+        zoneNumber = -1;
+    }
+
+    return zoneNumber;
+}
+
+int KKJZone::getZoneNumberFromLongitude(const KKJGeoCoordinate &kkj) {
+    // determine the zonenumber from KKJ easting
+    // takes KKJ zone which has center meridian
+    // longitude nearest (in math value) to
+    // the given KKJ longitude
+    int zoneNumber = MAX_ZONE;
+    while (zoneNumber >= MIN_ZONE) {
+        if (fabs(kkj.longitude() - KKJ_ZONE_INFO[zoneNumber][0]) <= 1.5) {
+            break;
+        }
+        zoneNumber--;
+    }
+
+    return zoneNumber;
+}
+
+double KKJZone::getCentralMeridianOfZone(int zoneNumber)
+{
+    if (zoneNumber >= MIN_ZONE && zoneNumber <= MAX_ZONE) {
+        return KKJ_ZONE_INFO[zoneNumber][0];
+    }
+
+    return 0.0;
+}
+
+double KKJZone::getFalseEastingOfZone(int zoneNumber)
+{
+    if (zoneNumber >= MIN_ZONE && zoneNumber <= MAX_ZONE) {
+        return KKJ_ZONE_INFO[zoneNumber][1];
+    }
+
+    return 0.0;
+}
+
+
+/**
+ * Transforms a KKJ geographical coordinate to a WGS84 geographical coordinate.
+ * @param fromCoordinate the input coordinate
+ * @return the transformed coordinate
+ */
+QGeoCoordinate transformToWGS84GeoCoordinate(const KKJGeoCoordinate &fromCoordinate) {
+    double kkjla = fromCoordinate.latitude();
+    double kkjlo = fromCoordinate.longitude();
+
+    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;
+    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;
+
+    return QGeoCoordinate(kkjla + dLa, kkjlo + dLo);
+}
+
+
+/**
+ * Transforms a WGS84 geographical coordinate to a KKJ geographical coordinate.
+ * @param fromCoordinate the input coordinate
+ * @return the transformed coordinate
+ */
+KKJGeoCoordinate transformToKKJGeoCoordinate(const QGeoCoordinate &fromCoordinate) {
+    double longitude = fromCoordinate.longitude();
+    double latitude = fromCoordinate.latitude();
+
+    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;
+    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;
+
+    return KKJGeoCoordinate(latitude + dLa, longitude + dLo);
+}
+
+
+/**
+ * Transforms a KKJ geographical coordinate to a KKJ rectangular grid coordinate.
+ * @param fromCoordinate the input coordinate
+ * @param zoneNumber the zone number in which the input coordinate resides
+ * @return the transformed coordinate
+ */
+KKJ transformToKKJGridCoordinate(const KKJGeoCoordinate &fromCoordinate) {
+    int zoneNumber = KKJZone::getZoneNumberFromLongitude(fromCoordinate);
+    double Lo = radians(fromCoordinate.longitude()) - radians(KKJZone::getCentralMeridianOfZone(zoneNumber));
+    double cosLa = cos(radians(fromCoordinate.latitude()));
+    double NN = HayfordEllipsoid::ee * cosLa * cosLa;
+    double LaF = atan(tan(radians(fromCoordinate.latitude())) / cos(Lo * sqrt(1.0 + NN)));
+    double cosLaF = cos(LaF);
+    double t = (tan(Lo) * cosLaF) / sqrt(1.0 + HayfordEllipsoid::ee * cosLaF * cosLaF);
+    double A = HayfordEllipsoid::a / (1.0 + HayfordEllipsoid::n);
+    double A1 = A * (1.0 + HayfordEllipsoid::nn / 4.0 + HayfordEllipsoid::nn * HayfordEllipsoid::nn / 64.0);
+    double A2 = A * 1.5 * HayfordEllipsoid::n * (1.0 - HayfordEllipsoid::nn / 8.0);
+    double A3 = A * 0.9375 * HayfordEllipsoid::nn * (1.0 - HayfordEllipsoid::nn / 4.0);
+    double A4 = A * 35.0 / 48.0 * HayfordEllipsoid::nn * HayfordEllipsoid::n;
+
+    unsigned int outY = A1 * LaF - A2 * sin(2.0 * LaF) + A3 * sin(4.0 * LaF) - A4 * sin(6.0 * LaF);
+    unsigned int outX = HayfordEllipsoid::c * log(t + sqrt(1.0 + t * t)) + 500000.0 + zoneNumber * 1000000.0;
+    return KKJ(outY, outX);
+}
+
+/**
+ * Transforms a KKJ rectangular grid coordinate to a KKJ geographical coordinate.
+ * @param fromCoordinate the input coordinate
+ * @return the transformed coordinate
+ */
+KKJGeoCoordinate transformToKKJGeoCoordinate(const KKJ &fromCoordinate) {
+    // Scan iteratively the target area, until find matching
+    // KKJ coordinate value. Area is defined with Hayford Ellipsoid.
+    double minLo = 18.5;
+    double maxLo = 32.0;
+    double minLa = 59.0;
+    double maxLa = 70.5;
+
+    int i = 1;
+
+    KKJGeoCoordinate ret;
+
+    while (i < 35) {
+        double deltaLo = maxLo - minLo;
+        double deltaLa = maxLa - minLa;
+        ret.setLongitude(minLo + 0.5 * deltaLo);
+        ret.setLatitude(minLa + 0.5 * deltaLa);
+        KKJ kkj = transformToKKJGridCoordinate(ret);
+        if (kkj.northing() < fromCoordinate.northing()) {
+            minLa = minLa + 0.45 * deltaLa;
+        } else {
+            maxLa = minLa + 0.55 * deltaLa;
+        }
+
+        if (kkj.easting() < fromCoordinate.easting()) {
+            minLo = minLo + 0.45 * deltaLo;
+        } else {
+            maxLo = minLo + 0.55 * deltaLo;
+        }
+
+        i++;
+    }
+
+    return ret;
+}
+
+
+KKJ CoordinateSystemTransformer::transformToKKJ(const QGeoCoordinate &fromCoordinate)
+{
+    KKJGeoCoordinate tmpKKJ = transformToKKJGeoCoordinate(fromCoordinate);
+    return transformToKKJGridCoordinate(tmpKKJ);
+}
+
+QGeoCoordinate CoordinateSystemTransformer::transformToWGS84(const KKJ &fromCoordinate)
+{
+    KKJGeoCoordinate tmpKKJ = transformToKKJGeoCoordinate(fromCoordinate);
+    return transformToWGS84GeoCoordinate(tmpKKJ);
+}

diff --git a/src/coordinatesystemtransformer.h b/src/coordinatesystemtransformer.h
new file mode 100644 (file)
index 0000000..b8f9dfe
--- /dev/null
+++ b/src/coordinatesystemtransformer.h
@@ -0,0 +1,30 @@
+#ifndef COORDINATESYSTEMTRANSFORMER_H
+#define COORDINATESYSTEMTRANSFORMER_H
+
+#include "kkj.h"
+#include <QGeoCoordinate>
+
+/**
+ * A utility class for transforming coordinates from one coordinate system to another.
+ */
+class CoordinateSystemTransformer
+{
+public:
+    /**
+     * Makes a coordinate transformation from WGS84 coordinate system to KKJ rectangular
+     * grid coordinates.
+     * @param fromCoordinate the WGS84 coordinate that will be transformed.
+     * @return the transformed coordinate in KKJ coordinate system.
+     */
+    static KKJ transformToKKJ(const QTM_NAMESPACE::QGeoCoordinate &fromCoordinate);
+
+    /**
+     * Makes a coordinate transformation from KKJ rectangular grid coordinate system
+     * to WGS84 coordinate system.
+     * @param fromCoordinate the KKJ coordinate that will be transformed.
+     * @return the transformed coordinate in WGS84 coordinate system.
+     */
+    static QTM_NAMESPACE::QGeoCoordinate transformToWGS84(const KKJ &fromCoordinate);
+};
+
+#endif // COORDINATESYSTEMTRANSFORMER_H

diff --git a/tests/tests.pro b/tests/tests.pro
index 2e22c4a..573f500 100644 (file)
--- a/tests/tests.pro
+++ b/tests/tests.pro
@@ -1,5 +1,6 @@
 TEMPLATE = subdirs
 SUBDIRS = \
+    ut_coordinatesystemtransformer \
     ut_kkj
 
 check.target = check

diff --git a/tests/ut_common.pri b/tests/ut_common.pri
index 3937725..390dd94 100644 (file)
--- a/tests/ut_common.pri
+++ b/tests/ut_common.pri
@@ -1,3 +1,5 @@
 SRCDIR = ../../src
+MOCKSDIR = ../mocks
 INCLUDEPATH += ../util \
-    $$SRCDIR
+    $$SRCDIR \
+    $$MOCKSDIR

diff --git a/tests/ut_coordinatesystemtransformer/.gitignore b/tests/ut_coordinatesystemtransformer/.gitignore
new file mode 100644 (file)
index 0000000..dd22e19
--- /dev/null
+++ b/tests/ut_coordinatesystemtransformer/.gitignore
@@ -0,0 +1,2 @@
+Makefile
+ut_coordinatesystemtransformer

diff --git a/tests/ut_coordinatesystemtransformer/ut_coordinatesystemtransformer.cpp b/tests/ut_coordinatesystemtransformer/ut_coordinatesystemtransformer.cpp
new file mode 100644 (file)
index 0000000..4658d60
--- /dev/null
+++ b/tests/ut_coordinatesystemtransformer/ut_coordinatesystemtransformer.cpp
@@ -0,0 +1,73 @@
+#include "coordinatesystemtransformer.h"
+
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+#include "stlhelpers.h"
+
+#include <QPair>
+
+class CoordinateSystemTransformerTest : public ::testing::Test
+{
+public:
+    QList<QPair<QTM_NAMESPACE::QGeoCoordinate, KKJ> > testData;
+
+    CoordinateSystemTransformerTest() {
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.2528, 25.02051), KKJ(6682815, 2556686));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.18713, 24.832), KKJ(6675352, 2546340));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.27414, 25.04465), KKJ(6685213, 2557985));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.2507, 25.01767), KKJ(6682578, 2556532));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.2902, 24.44804), KKJ(6686629, 2524959));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.35033, 25.06718), KKJ(6693721, 2559094));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.25471, 25.02373), KKJ(6683030, 2556861));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.25417, 25.0242), KKJ(6682971, 2556888));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.39737, 25.08981), KKJ(6698983, 2560257));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.28923, 25.12709), KKJ(6686969, 2562518));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.1727, 24.65643), KKJ(6673635, 2536615));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.35133, 25.06764), KKJ(6693833, 2559118));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.34949, 25.06874), KKJ(6693629, 2559182));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.25119, 25.02518), KKJ(6682640, 2556947));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.25196, 25.02294), KKJ(6682723, 2556822));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.34929, 25.06705), KKJ(6693605, 2559089));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.18855, 24.83393), KKJ(6675512, 2546445));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.25091, 25.02547), KKJ(6682609, 2556964));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.25321, 25.0215), KKJ(6682861, 2556740));
+        testData << qMakePair(QTM_NAMESPACE::QGeoCoordinate(60.35291, 25.06559), KKJ(6694007, 2559002));
+    }
+};
+
+TEST_F(CoordinateSystemTransformerTest, WGS84CoordinatesToKKJCoordinates)
+{
+    QListIterator<QPair<QTM_NAMESPACE::QGeoCoordinate, KKJ> > it(testData);
+    while (it.hasNext()) {
+        const QPair<QTM_NAMESPACE::QGeoCoordinate, KKJ> &datum = it.next();
+        KKJ result = CoordinateSystemTransformer::transformToKKJ(datum.first);
+        KKJ expected = datum.second;
+        // Allow one unit difference from the expected
+        int northDiff = abs((long)expected.northing() - (long)result.northing());
+        int eastDiff = abs((long)expected.easting() - (long)result.easting());
+        EXPECT_LE(northDiff, 1);
+        EXPECT_LE(eastDiff, 1);
+    }
+}
+
+TEST_F(CoordinateSystemTransformerTest, KKJCoordinatesToWGS84Coordinates)
+{
+    QListIterator<QPair<QTM_NAMESPACE::QGeoCoordinate, KKJ> > it(testData);
+    while (it.hasNext()) {
+        const QPair<QTM_NAMESPACE::QGeoCoordinate, KKJ> &datum = it.next();
+        QTM_NAMESPACE::QGeoCoordinate result = CoordinateSystemTransformer::transformToWGS84(datum.second);
+        QTM_NAMESPACE::QGeoCoordinate expected = datum.first;
+        // Allow small difference from the expected
+        double latitudeDiff = fabs(expected.latitude() - result.latitude());
+        double longitudeDiff = fabs(expected.longitude() - result.longitude());
+        EXPECT_LE(latitudeDiff, 0.00001);
+        EXPECT_LE(longitudeDiff, 0.00001);
+    }
+}
+
+int main(int argc, char *argv[])
+{
+    ::testing::InitGoogleMock(&argc, argv);
+    return RUN_ALL_TESTS();
+}

diff --git a/tests/ut_coordinatesystemtransformer/ut_coordinatesystemtransformer.pro b/tests/ut_coordinatesystemtransformer/ut_coordinatesystemtransformer.pro
new file mode 100644 (file)
index 0000000..213f01d
--- /dev/null
+++ b/tests/ut_coordinatesystemtransformer/ut_coordinatesystemtransformer.pro
@@ -0,0 +1,22 @@
+include(../ut_common.pri)
+
+TARGET = ut_coordinatesystemtransformer
+QT += testlib
+QT -= gui
+CONFIG += console \
+    mobility
+CONFIG -= app_bundle
+MOBILITY = location
+TEMPLATE = app
+OBJECTS_DIR = .obj
+MOC_DIR = .moc
+SOURCES += \
+    ut_coordinatesystemtransformer.cpp \
+    $$SRCDIR/coordinatesystemtransformer.cpp \
+    $$SRCDIR/kkj.cpp
+HEADERS += \
+    $$SRCDIR/coordinatesystemtransformer.h \
+    $$SRCDIR/kkj.h
+
+include(../gmock.pri)
+include(../check.pri)

diff --git a/tests/util/stlhelpers.h b/tests/util/stlhelpers.h
new file mode 100644 (file)
index 0000000..d51798b
--- /dev/null
+++ b/tests/util/stlhelpers.h
@@ -0,0 +1,20 @@
+#ifndef STLHELPERS_H
+#define STLHELPERS_H
+
+#include <ostream>
+#include "kkj.h"
+#include <QGeoCoordinate>
+
+std::ostream& operator<<(std::ostream& stream, const KKJ &val)
+{
+    stream << val.northing() << ", " << val.easting();
+    return stream;
+}
+
+std::ostream& operator<<(std::ostream& stream, const QTM_NAMESPACE::QGeoCoordinate &val)
+{
+    stream << val.latitude() << ", " << val.longitude();
+    return stream;
+}
+
+#endif // STLHELPERS_H