--- /dev/null
+#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);
+}
projects / ptas / blobdiff