Implemented KKJ<->WGS84 transformations in C.

This is a direct translation from the Python prototype. There is some
cleaning up to be done next. Includes unit test.

diff --git a/.gitignore b/.gitignore
index 8705fec..5a90e3d 100644 (file)
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,5 @@
 .deps
 .libs
-Makefile
 Makefile.in
 aclocal.m4
 autom4te.cache
@@ -21,3 +20,4 @@ ltmain.sh
 missing
 stamp-h1
 *.desktop
+*~

diff --git a/coordinate-system.c b/coordinate-system.c
index fefb1c4..ed84894 100644 (file)
--- a/coordinate-system.c
+++ b/coordinate-system.c
@@ -1,21 +1,156 @@
 #include "coordinate-system.h"
 
-void wgs84_to_kkj2(double inLatitude, double inLongitude, double *outLatitude, double *outLongitude)
-{
-    if (outLatitude) {
-        *outLatitude = 100;
-    }
-    if (outLongitude) {
-        *outLongitude = 100;
+#include <math.h>
+
+
+// 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;
+}
+
+
+// Constants
+// Longitude0 and Center meridian of KKJ bands
+const double KKJ_ZONE_INFO[6][2] = { {18.0,  500000.0},
+                                     {21.0, 1500000.0},
+                                     {24.0, 2500000.0},
+                                     {27.0, 3500000.0},
+                                     {30.0, 4500000.0},
+                                     {33.0, 5500000.0} };
+
+
+// Function:  KKJ_Zone_I
+int KKJ_Zone_I(KKJ easting) {
+    int zoneNumber = floor(easting / 1000000.0);
+    if (zoneNumber < 0 || zoneNumber > 5) {
+        zoneNumber = -1;
     }
+
+    return zoneNumber;
 }
 
-void kkj2_to_wgs84(double inLatitude, double inLongitude, double *outLatitude, double *outLongitude)
-{
-    if (outLatitude) {
-        *outLatitude = 100;
+
+// Function:  KKJ_Zone_Lo
+int KKJ_Zone_Lo(double kkjlo) {
+    // 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 = 5;
+    while (zoneNumber >= 0) {
+        if (fabs(kkjlo - KKJ_ZONE_INFO[zoneNumber][0]) <= 1.5) {
+            break;
+        }
+        zoneNumber--;
     }
-    if (outLongitude) {
-        *outLongitude = 100;
+
+    return zoneNumber;
+}
+
+
+// Function:  KKJlalo_to_WGS84lalo
+void KKJlola_to_WGS84lola(double kkjlo, double kkjla, double *outLongitude, double *outLatitude) {
+    double dLa = radians(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 = radians(-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;
+
+    *outLatitude = degrees(radians(kkjla) + dLa);
+    *outLongitude = degrees(radians(kkjlo) + dLo);
+}
+
+
+// Function:  WGS84lalo_to_KKJlalo
+void WGS84lola_to_KKJlola(double longitude, double latitude, double *outLongitude, double *outLatitude) {
+    double dLa = radians(-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 = radians(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;
+
+    *outLatitude = degrees(radians(latitude) + dLa);
+    *outLongitude = degrees(radians(longitude) + dLo);
+}
+
+
+// Function:  KKJlalo_to_KKJxy
+void KKJlola_to_KKJxy(double lon, double lat, int zoneNumber, KKJ *outX, KKJ *outY) {
+    // Hayford ellipsoid
+    double a = 6378388.0;
+    double f  = 1.0 / 297.0;
+    double b  = (1.0 - f) * a;
+    double bb = b * b;
+    double c  = (a / b) * a;
+    double ee = (a * a - bb) / bb;
+    double n = (a - b) / (a + b);
+    double nn = n * n;
+
+    double Lo = radians(lon) - radians(KKJ_ZONE_INFO[zoneNumber][0]);
+    double cosLa = cos(radians(lat));
+    double NN = ee * cosLa * cosLa;
+    double LaF = atan(tan(radians(lat)) / cos(Lo * sqrt(1.0 + NN)));
+    double cosLaF = cos(LaF);
+    double t = (tan(Lo) * cosLaF) / sqrt(1.0 + ee * cosLaF * cosLaF);
+    double A = a / (1.0 + n);
+    double A1 = A * (1.0 + nn / 4.0 + nn * nn / 64.0);
+    double A2 = A * 1.5 * n * (1.0 - nn / 8.0);
+    double A3 = A * 0.9375 * nn * (1.0 - nn / 4.0);
+    double A4 = A * 35.0 / 48.0 * nn * n;
+
+    *outY = A1 * LaF - A2 * sin(2.0 * LaF) + A3 * sin(4.0 * LaF) - A4 * sin(6.0 * LaF);
+    *outX = c * log(t + sqrt(1.0 + t * t)) + 500000.0 + zoneNumber * 1000000.0;
+}
+
+// Function:  KKJxy_to_KKJlalo
+void KKJxy_to_KKJlola(KKJ x, KKJ y, double *outLongitude, double *outLatitude) {
+    // Scan iteratively the target area, until find matching
+    // KKJ coordinate value.  Area is defined with Hayford Ellipsoid.
+    int zoneNumber = KKJ_Zone_I(x);
+    double minLo = radians(18.5);
+    double maxLo = radians(32.0);
+    double minLa = radians(59.0);
+    double maxLa = radians(70.5);
+
+    int i = 1;
+    KKJ tmpX, tmpY;
+
+    while (i < 35) {
+        double deltaLo = maxLo - minLo;
+        double deltaLa = maxLa - minLa;
+        *outLongitude = degrees(minLo + 0.5 * deltaLo);
+        *outLatitude = degrees(minLa + 0.5 * deltaLa);
+        KKJlola_to_KKJxy(*outLongitude, *outLatitude, zoneNumber, &tmpX, &tmpY);
+        if (tmpY < y) {
+            minLa = minLa + 0.45 * deltaLa;
+        } else {
+            maxLa = minLa + 0.55 * deltaLa;
+        }
+
+        if (tmpX < x) {
+            minLo = minLo + 0.45 * deltaLo;
+        } else {
+            maxLo = minLo + 0.55 * deltaLo;
+        }
+
+        i++;
     }
 }
+
+  ////////////////////
+ // PUBLIC METHODS //
+////////////////////
+
+void WGS84lola_to_KKJxy(double longitude, double latitude, KKJ *outX, KKJ *outY) {
+    double kkjlo, kkjla;
+
+    WGS84lola_to_KKJlola(longitude, latitude, &kkjlo, &kkjla);
+    int zoneNumber = KKJ_Zone_Lo(kkjlo);
+    KKJlola_to_KKJxy(kkjlo, kkjla, zoneNumber, outX, outY);
+ }
+
+void KKJxy_to_WGS84lola(KKJ x, KKJ y, double *outLongitude, double *outLatitude) {
+    double kkjlo, kkjla;
+
+    KKJxy_to_KKJlola(x, y, &kkjlo, &kkjla);
+    KKJlola_to_WGS84lola(kkjlo, kkjla, outLongitude, outLatitude);
+}

diff --git a/coordinate-system.h b/coordinate-system.h
index 78a4620..7c3d100 100644 (file)
--- a/coordinate-system.h
+++ b/coordinate-system.h
@@ -1,8 +1,25 @@
 #ifndef COORDINATE_SYSTEM_H
 #define COORDINATE_SYSTEM_H
 
-void wgs84_to_kkj2(double inLatitude, double inLongitude, double *outLatitude, double *outLongitude);
+// Type for KKJ x/y coordinates
+typedef unsigned int KKJ;
 
-void kkj2_to_wgs84(double inLatitude, double inLongitude, double *outLatitude, double *outLongitude);
+/**
+ * Transformes WGS84 longitude/latitude coordinates to KKJ x/y coordinates.
+ * @param longitude the input longitude in degrees
+ * @param latitude the input latitude in degrees
+ * @param outX the result x (easting)
+ * @param outY the result y (northing)
+ */
+void WGS84lola_to_KKJxy(double longitude, double latitude, KKJ *outX, KKJ *outY);
+
+/**
+ * Transformes KKJ x/y coordinates to WGS84 longitude/latitude coordinates.
+ * @param x the input x (easting)
+ * @param y the input y (northing)
+ * @param outLongitude the result longitude in degrees
+ * @param outLatitude the result latitude in degrees
+ */
+void KKJxy_to_WGS84lola(KKJ x, KKJ y, double *outLongitude, double *outLatitude);
 
 #endif

diff --git a/tests/ut_coord_trans/.gitignore b/tests/ut_coord_trans/.gitignore
new file mode 100644 (file)
index 0000000..4ec7554
--- /dev/null
+++ b/tests/ut_coord_trans/.gitignore
@@ -0,0 +1 @@
+ut_coord_trans

diff --git a/tests/ut_coord_trans/Makefile b/tests/ut_coord_trans/Makefile
new file mode 100644 (file)
index 0000000..095515a
--- /dev/null
+++ b/tests/ut_coord_trans/Makefile
@@ -0,0 +1,29 @@
+CC=gcc
+RM=rm -f
+INCDIR=../../
+CFLAGS=-I$(INCDIR)
+LDFLAGS=-lm
+DEPS=$(INCDIR)coordinate-system.h
+
+UNIT_TEST_SOURCES=ut_coord_trans.c
+SOURCES=$(UNIT_TEST_SOURCES) \
+  $(INCDIR)coordinate-system.c
+TARGET=ut_coord_trans
+
+OBJS=$(SOURCES:%.c=%.o)
+
+all: $(TARGET)
+
+%.o: %.c $(DEPS)
+       $(CC) -c -o $@ $(CFLAGS) $<
+
+$(TARGET): $(OBJS)
+       $(CC) -o $(TARGET) $(LDFLAGS) $(OBJS)
+
+.PHONY: clean
+
+clean:
+       $(RM) *.o *~ core $(TARGET) 
+
+check: $(TARGET)
+       @./$(TARGET)

diff --git a/tests/ut_coord_trans/ut_coord_trans.c b/tests/ut_coord_trans/ut_coord_trans.c
new file mode 100644 (file)
index 0000000..425d629
--- /dev/null
+++ b/tests/ut_coord_trans/ut_coord_trans.c
@@ -0,0 +1,96 @@
+#include "coordinate-system.h"
+
+#include <math.h>
+#include <stdio.h>
+
+struct TestCoordinate {
+    KKJ x, y;
+    double lon, lat;
+};
+
+// Test data extracted from example on page
+// http://developer.reittiopas.fi/pages/fi/http-get-interface.php
+struct TestCoordinate testData[] = {
+    { 2556686, 6682815, 25.02051, 60.2528 },
+    { 2546340, 6675352, 24.832, 60.18713 },
+    { 2557985, 6685213, 25.04465, 60.27414 },
+    { 2556532, 6682578, 25.01767, 60.2507 },
+    { 2524959, 6686629, 24.44804, 60.2902 },
+    { 2559094, 6693721, 25.06718, 60.35033 },
+    { 2556861, 6683030, 25.02373, 60.25471 },
+    { 2556888, 6682971, 25.0242, 60.25417 },
+    { 2560257, 6698983, 25.08981, 60.39737 },
+    { 2562518, 6686969, 25.12709, 60.28923 },
+    { 2536615, 6673635, 24.65643, 60.1727 },
+    { 2559118, 6693833, 25.06764, 60.35133 },
+    { 2559182, 6693629, 25.06874, 60.34949 },
+    { 2556947, 6682640, 25.02518, 60.25119 },
+    { 2556822, 6682723, 25.02294, 60.25196 },
+    { 2559089, 6693605, 25.06705, 60.34929 },
+    { 2546445, 6675512, 24.83393, 60.18855 },
+    { 2556964, 6682609, 25.02547, 60.25091 },
+    { 2556740, 6682861, 25.0215, 60.25321 },
+    { 2559002, 6694007, 25.06559, 60.35291 }};
+
+
+int testKKJxytoWGS84lola() {
+    double lon, lat;
+    int result = 0;
+    int i;
+
+    for (i = 0; i < sizeof(testData) / sizeof(struct TestCoordinate); ++i) {
+        KKJxy_to_WGS84lola(testData[i].x, testData[i].y, &lon, &lat);
+
+        if (fabs(testData[i].lon - lon) < 0.001 && fabs(testData[i].lat - lat) < 0.001) {
+            ;
+        } else {
+            printf("Got: (%f, %f), expected: (%f, %f)\n", lon, lat, testData[i].lon, testData[i].lat);
+            result = -1;
+        }
+    }
+
+    return result;
+}
+
+
+int testWGS84lolatoKKJxy() {
+    KKJ x, y;
+    int result = 0;
+    int i;
+
+    for (i = 0; i < sizeof(testData) / sizeof(struct TestCoordinate); ++i) {
+        WGS84lola_to_KKJxy(testData[i].lon, testData[i].lat, &x, &y);
+
+        if (abs(testData[i].x - x) < 2 && abs(testData[i].y - y) < 2) {
+            ;
+        } else {
+            printf("Got: (%u, %u), expected: (%u, %u)\n", x, y, testData[i].x, testData[i].y);
+            result = -1;
+        }
+    }
+
+    return result;
+}
+
+
+
+int main(int argc, char* argv[]) {
+    int testResult = 0;
+
+    // Test transforming from KKJxy to WGS84lonlat
+    testResult = testKKJxytoWGS84lola();
+
+    if (testResult == 0) {
+        printf("All tests in testKKJxytoWGS84lola passed\n");
+    }
+
+    // Test transforming from WGS84lonlat to KKJxy
+    testResult = testWGS84lolatoKKJxy();
+
+    if (testResult == 0) {
+        printf("All tests in testWGS84lolatoKKJxy passed\n");
+    }
+
+
+    return 0;
+}