From 4c16ea9637aef4ee32b83af36eb86e5e85a889ea Mon Sep 17 00:00:00 2001 From: Aki Koskinen Date: Sun, 7 Feb 2010 12:21:07 +0200 Subject: [PATCH] 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. --- .gitignore | 2 +- coordinate-system.c | 161 ++++++++++++++++++++++++++++++--- coordinate-system.h | 21 ++++- tests/ut_coord_trans/.gitignore | 1 + tests/ut_coord_trans/Makefile | 29 ++++++ tests/ut_coord_trans/ut_coord_trans.c | 96 ++++++++++++++++++++ 6 files changed, 294 insertions(+), 16 deletions(-) create mode 100644 tests/ut_coord_trans/.gitignore create mode 100644 tests/ut_coord_trans/Makefile create mode 100644 tests/ut_coord_trans/ut_coord_trans.c diff --git a/.gitignore b/.gitignore index 8705fec..5a90e3d 100644 --- 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 --- 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 + + +// 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 --- 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 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 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 index 0000000..425d629 --- /dev/null +++ b/tests/ut_coord_trans/ut_coord_trans.c @@ -0,0 +1,96 @@ +#include "coordinate-system.h" + +#include +#include + +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; +} -- 1.7.9.5