From: Aki Koskinen Date: Sun, 21 Mar 2010 15:45:37 +0000 (+0200) Subject: Coordinate system transformer X-Git-Url: http://vcs.maemo.org/git/?p=ptas;a=commitdiff_plain;h=ba2b89b302c4e9d32edddc9b1d1845dcd967dce4 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 index 0000000..1632133 --- /dev/null +++ b/src/coordinatesystemtransformer.cpp @@ -0,0 +1,348 @@ +#include "coordinatesystemtransformer.h" + +#include + + +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 index 0000000..b8f9dfe --- /dev/null +++ b/src/coordinatesystemtransformer.h @@ -0,0 +1,30 @@ +#ifndef COORDINATESYSTEMTRANSFORMER_H +#define COORDINATESYSTEMTRANSFORMER_H + +#include "kkj.h" +#include + +/** + * 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 --- 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 --- 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 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 index 0000000..4658d60 --- /dev/null +++ b/tests/ut_coordinatesystemtransformer/ut_coordinatesystemtransformer.cpp @@ -0,0 +1,73 @@ +#include "coordinatesystemtransformer.h" + +#include +#include + +#include "stlhelpers.h" + +#include + +class CoordinateSystemTransformerTest : public ::testing::Test +{ +public: + QList > 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 > it(testData); + while (it.hasNext()) { + const QPair &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 > it(testData); + while (it.hasNext()) { + const QPair &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 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 index 0000000..d51798b --- /dev/null +++ b/tests/util/stlhelpers.h @@ -0,0 +1,20 @@ +#ifndef STLHELPERS_H +#define STLHELPERS_H + +#include +#include "kkj.h" +#include + +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