package com.ximple.eofms.util;
|
|
import com.vividsolutions.jts.geom.Coordinate;
|
|
/**
|
* TWDDatumConverter
|
* User: Ulysses
|
* Date: 2007/10/8
|
* To change this template use File | Settings | File Templates.
|
*/
|
public abstract class TWDDatumConverter {
|
/*
|
* Definition of math related value
|
*/
|
private static final double COS67_5 = 0.3826834323650897717284599840304e0;
|
private static final double PI = 3.14159265358979323e0;
|
private static final double HALF_PI = 1.570796326794896615e0;
|
private static final double DEG_RAD = 0.01745329251994329572e0;
|
private static final double RAD_DEG = 57.295779513082321031e0;
|
|
/*
|
* Definition of datum related value
|
*/
|
private static final double AD_C = 1.0026000e0;
|
private static final double TWD67_A = 6378160.0e0;
|
private static final double TWD67_B = 6356774.7192e0;
|
private static final double TWD67_ECC = 0.00669454185458e0;
|
private static final double TWD67_ECC2 = 0.00673966079586e0;
|
|
// different from garmin and already knowned value, but those all value only
|
private static final double TWD67_DX = -752.32e0;
|
|
// got 5-15m accuracy. the real offical value is holded by somebody and not
|
private static final double TWD67_DY = -361.32e0;
|
|
// release to public. if can got more enough twd67/twd97 control point coordinare,
|
private static final double TWD67_DZ = -180.51e0;
|
private static final double TWD67_RX = -0.00000117e0; // then we can calculate a better value than now.
|
private static final double TWD67_RY = 0.00000184e0; //
|
private static final double TWD67_RZ = 0.00000098e0; // and, also lack twd67/twd97 altitude convertion value...
|
private static final double TWD67_S = 0.00002329e0; //
|
private static final double TWD97_A = 6378137.0e0;
|
private static final double TWD97_B = 6356752.3141e0;
|
private static final double TWD97_ECC = 0.00669438002290e0;
|
private static final double TWD97_ECC2 = 0.00673949677556e0;
|
private static final double TWD67_TM2 = 0.9999e0; // TWD67->TM2 scale
|
private static final double TWD97_TM2 = 0.9999e0; // TWD97->TM2 scale
|
|
/*
|
* datum convert function
|
*/
|
public static Coordinate toTWD97(Coordinate pt) {
|
double newX, newY, newZ;
|
double r, pole, sin_lat, cos_lat;
|
double lat, lon, height;
|
double x1, y1, z1, x2, y2, z2;
|
double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;
|
|
lon = pt.x * DEG_RAD;
|
lat = pt.y * DEG_RAD;
|
height = pt.z * DEG_RAD;
|
|
if ((lat < -HALF_PI) && (lat > -1.001 * HALF_PI)) {
|
lat = -HALF_PI;
|
} else if ((lat > HALF_PI) && (lat < 1.001 * HALF_PI)) {
|
lat = HALF_PI;
|
} else if ((lat < -HALF_PI) || (lat > HALF_PI)) {
|
return null;
|
}
|
|
if (lon > PI) {
|
lon -= (2 * PI);
|
}
|
|
sin_lat = Math.sin(lat);
|
cos_lat = Math.cos(lat);
|
r = TWD67_A / (Math.sqrt(1.0 - TWD67_ECC * sin_lat * sin_lat));
|
x1 = (r + height) * cos_lat * Math.cos(lon);
|
y1 = (r + height) * cos_lat * Math.sin(lon);
|
z1 = ((r * (1 - TWD67_ECC)) + height) * sin_lat;
|
x2 = x1 + TWD67_DX + TWD67_S * (lon + TWD67_RZ * lat - TWD67_RY * height);
|
y2 = y1 + TWD67_DY + TWD67_S * (-TWD67_RZ * lon + lat + TWD67_RX * height);
|
z2 = z1 + TWD67_DZ + TWD67_S * (TWD67_RY * lon - TWD67_RX * lat + height);
|
pole = 0.0;
|
|
if (x2 != 0.0) {
|
lon = Math.atan2(y2, x2);
|
} else {
|
if (y2 > 0) {
|
lon = HALF_PI;
|
} else if (y2 < 0) {
|
lon = -HALF_PI;
|
} else {
|
pole = 1;
|
lon = 0;
|
|
if (z2 > 0) {
|
lat = HALF_PI;
|
} else if (z2 < 0) {
|
lat = -HALF_PI;
|
} else {
|
lat = HALF_PI;
|
newX = lon * RAD_DEG;
|
newY = lat * RAD_DEG;
|
newZ = -TWD97_B;
|
|
return new Coordinate(newX, newY, newZ);
|
}
|
}
|
}
|
|
q2 = x2 * x2 + y2 * y2;
|
q = Math.sqrt(q2);
|
t = z2 * AD_C;
|
s = Math.sqrt(t * t + q2);
|
sin_b = t / s;
|
cos_b = q / s;
|
t1 = z2 + TWD97_B * TWD97_ECC2 * sin_b * sin_b * sin_b;
|
sum = q - TWD97_A * TWD97_ECC * cos_b * cos_b * cos_b;
|
s1 = Math.sqrt(t1 * t1 + sum * sum);
|
sin_p = t1 / s1;
|
cos_p = sum / s1;
|
r = TWD97_A / Math.sqrt(1.0 - TWD97_ECC * sin_p * sin_p);
|
|
if (cos_p >= COS67_5) {
|
height = q / cos_p - r;
|
} else if (cos_p <= -COS67_5) {
|
height = q / -cos_p - r;
|
} else {
|
height = z2 / sin_p + r * (TWD97_ECC - 1.0);
|
}
|
|
if (pole != 0.0) {
|
lat = Math.atan(sin_p / cos_p);
|
}
|
|
newX = lon * RAD_DEG;
|
newY = lat * RAD_DEG;
|
newZ = height;
|
return new Coordinate(newX, newY, newZ);
|
}
|
|
public static Coordinate toTWD67(Coordinate pt) {
|
double newX, newY, newZ;
|
double r, pole, sin_lat, cos_lat;
|
double lat, lon, height;
|
double x1, y1, z1, x2, y2, z2;
|
double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;
|
|
lon = pt.x * DEG_RAD;
|
lat = pt.y * DEG_RAD;
|
height = pt.z * DEG_RAD;
|
|
if ((lat < -HALF_PI) && (lat > -1.001 * HALF_PI)) {
|
lat = -HALF_PI;
|
} else if ((lat > HALF_PI) && (lat < 1.001 * HALF_PI)) {
|
lat = HALF_PI;
|
} else if ((lat < -HALF_PI) || (lat > HALF_PI)) {
|
return null;
|
}
|
|
if (lon > PI) {
|
lon -= (2 * PI);
|
}
|
|
sin_lat = Math.sin(lat);
|
cos_lat = Math.cos(lat);
|
r = TWD97_A / (Math.sqrt(1.0 - TWD97_ECC * sin_lat * sin_lat));
|
x1 = (r + height) * cos_lat * Math.cos(lon);
|
y1 = (r + height) * cos_lat * Math.sin(lon);
|
z1 = ((r * (1 - TWD97_ECC)) + height) * sin_lat;
|
x2 = x1 - TWD67_DX - TWD67_S * (lon + TWD67_RZ * lat - TWD67_RY * height);
|
y2 = y1 - TWD67_DY - TWD67_S * (-TWD67_RZ * lon + lat + TWD67_RX * height);
|
z2 = z1 - TWD67_DZ - TWD67_S * (TWD67_RY * lon - TWD67_RX * lat + height);
|
pole = 0;
|
|
if (x2 != 0.0) {
|
lon = Math.atan2(y2, x2);
|
} else {
|
if (y2 > 0) {
|
lon = HALF_PI;
|
} else if (y2 < 0) {
|
lon = -HALF_PI;
|
} else {
|
pole = 1;
|
lon = 0;
|
|
if (z2 > 0) {
|
lat = HALF_PI;
|
} else if (z2 < 0) {
|
lat = -HALF_PI;
|
} else {
|
lat = HALF_PI;
|
newX = lon * RAD_DEG;
|
newY = lat * RAD_DEG;
|
newZ = -TWD67_B;
|
|
return new Coordinate(newX, newY, newZ);
|
}
|
}
|
}
|
|
q2 = x2 * x2 + y2 * y2;
|
q = Math.sqrt(q2);
|
t = z2 * AD_C;
|
s = Math.sqrt(t * t + q2);
|
sin_b = t / s;
|
cos_b = q / s;
|
t1 = z2 + TWD67_B * TWD67_ECC2 * sin_b * sin_b * sin_b;
|
sum = q - TWD67_A * TWD67_ECC * cos_b * cos_b * cos_b;
|
s1 = Math.sqrt(t1 * t1 + sum * sum);
|
sin_p = t1 / s1;
|
cos_p = sum / s1;
|
r = TWD67_A / Math.sqrt(1.0 - TWD67_ECC * sin_p * sin_p);
|
|
if (cos_p >= COS67_5) {
|
height = q / cos_p - r;
|
} else if (cos_p <= -COS67_5) {
|
height = q / -cos_p - r;
|
} else {
|
height = z2 / sin_p + r * (TWD67_ECC - 1.0);
|
}
|
|
if (pole != 0.0) {
|
lat = Math.atan(sin_p / cos_p);
|
}
|
|
newX = lon * RAD_DEG;
|
newY = lat * RAD_DEG;
|
newZ = height;
|
return new Coordinate(newX, newY, newZ);
|
}
|
|
public static Coordinate toTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double x, double y) {
|
double x0, y0, x1, y1, m0, m1;
|
double n, t, c, A;
|
double newX, newY;
|
|
x0 = x * DEG_RAD;
|
y0 = y * DEG_RAD;
|
x1 = lon * DEG_RAD;
|
y1 = lat * DEG_RAD;
|
m0 = mercator(y1, a, ecc);
|
m1 = mercator(y0, a, ecc);
|
n = a / Math.sqrt(1 - ecc * Math.pow(Math.sin(y0), 2.0));
|
t = Math.pow(Math.tan(y0), 2.0);
|
c = ecc2 * Math.pow(Math.cos(y0), 2.0);
|
A = (x0 - x1) * Math.cos(y0);
|
newX = scale * n
|
* (A + (1.0 - t + c) * A * A * A / 6.0
|
+ (5.0 - 18.0 * t + t * t + 72.0 * c - 58.0 * ecc2) * Math.pow(A, 5.0) / 120.0);
|
newY = scale
|
* (m1 - m0
|
+ n * Math.tan(y0)
|
* (A * A / 2.0 + (5.0 - t + 9.0 * c + 4 * c * c) * Math.pow(A, 4.0) / 24.0
|
+ (61.0 - 58.0 * t + t * t + 600.0 * c - 330.0 * ecc2) * Math.pow(A, 6.0) / 720.0));
|
return new Coordinate(newX, newY);
|
}
|
|
public static Coordinate fromTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double x, double y) {
|
double newX, newY;
|
double x0, y0, x1, y1, phi, m, m0, mu, e1;
|
double c1, t1, n1, r1, d;
|
|
x0 = x;
|
y0 = y;
|
x1 = lon * DEG_RAD;
|
y1 = lat * DEG_RAD;
|
m0 = mercator(y1, a, ecc);
|
m = m0 + y0 / scale;
|
e1 = (1.0 - Math.sqrt(1.0 - ecc)) / (1.0 + Math.sqrt(1.0 - ecc));
|
mu = m / (a * (1.0 - ecc / 4.0 - 3.0 * ecc * ecc / 64.0 - 5.0 * ecc * ecc * ecc / 256.0));
|
phi = mu + (3.0 * e1 / 2.0 - 27.0 * Math.pow(e1, 3.0) / 32.0) * Math.sin(2.0 * mu)
|
+ (21.0 * e1 * e1 / 16.0 - 55.0 * Math.pow(e1, 4.0) / 32.0) * Math.sin(4.0 * mu)
|
+ 151.0 * Math.pow(e1, 3.0) / 96.0 * Math.sin(6.0 * mu) + 1097.0 * Math.pow(e1, 4.0) / 512.0 * Math.sin(8.0 * mu);
|
c1 = ecc2 * Math.pow(Math.cos(phi), 2.0);
|
t1 = Math.pow(Math.tan(phi), 2.0);
|
n1 = a / Math.sqrt(1 - ecc * Math.pow(Math.sin(phi), 2.0));
|
r1 = a * (1.0 - ecc) / Math.pow(1.0 - ecc * Math.pow(Math.sin(phi), 2.0), 1.5);
|
d = x0 / (n1 * scale);
|
newX = (x1 + (d - (1.0 + 2.0 * t1 + c1) * Math.pow(d, 3.0) / 6.0
|
+ (5.0 - 2.0 * c1 + 28.0 * t1 - 3.0 * c1 * c1 + 8.0 * ecc2 + 24.0 * t1 * t1) * Math.pow(d, 5.0)
|
/ 120.0) / Math.cos(phi)) * RAD_DEG;
|
newY = (phi
|
- n1 * Math.tan(phi) / r1
|
* (d * d / 2.0 - (5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1 * c1 - 9.0 * ecc2) * Math.pow(d, 4.0) / 24.0
|
+ (61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1 * t1 - 252.0 * ecc2 - 3.0 * c1 * c1) * Math.pow(d, 6.0)
|
/ 72.0)) * RAD_DEG;
|
return new Coordinate(newX, newY);
|
}
|
|
private static double mercator(double y, double a, double ecc) {
|
if (y == 0.0) {
|
return 0.0;
|
} else {
|
return a * ((1.0 - ecc / 4.0 - 3.0 * ecc * ecc / 64.0 - 5.0 * ecc * ecc * ecc / 256.0) * y
|
- (3.0 * ecc / 8.0 + 3.0 * ecc * ecc / 32.0 + 45.0 * ecc * ecc * ecc / 1024.0) * Math.sin(2.0 * y)
|
+ (15.0 * ecc * ecc / 256.0 + 45.0 * ecc * ecc * ecc / 1024.0) * Math.sin(4.0 * y)
|
- (35.0 * ecc * ecc * ecc / 3072.0) * Math.sin(6.0 * y));
|
}
|
}
|
|
/**
|
*
|
* Sample code below, using coordinate in Dan Jacob's website.
|
*
|
*
|
* int main()
|
* {
|
* double x1, y1, z1, x2, y2, z2;
|
* double tx1, ty1, tx2, ty2;
|
* double dx, dy, dz, dx1, dy1;
|
*
|
* x1 = 120.85788004; // TWD67
|
* y1 = 24.18347242;
|
* z1 = 777;
|
*
|
* x2 = 120.86603958; // TWD97
|
* y2 = 24.18170479;
|
* z2 = 777;
|
*
|
* tx1 = 235561; // TWD67->TM2
|
* ty1 = 2675359;
|
*
|
* tx2 = 236389.849; // TWD97->TM2
|
* ty2 = 2675153.168;
|
*
|
* ////////////////////////////////////////////
|
* /
|
* /
|
* // convert TWD67->TM2
|
* /
|
* /
|
* ////////////////////////////////////////////
|
*
|
* dx = x1;
|
* dy = y1;
|
*
|
* toTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx, &dy);
|
* // center longitude of taiwan is 121, for penghu is 119
|
*
|
* dx += 250000; // TM2 in Taiwan should add 250000
|
*
|
* printf("TWD67->TM2nTWD67 (%f, %f)nConvert (%.3f, %.3f)nOrigin (%.3f, %.3f)n", x1, y1, dx, dy, tx1, ty1);
|
* printf("Acuuracy (%.3f, X:%.3f, Y:%.3f)nn", sqrt((dx-tx1)*(dx-tx1)+(dy-ty1)*(dy-ty1)), (dx-tx1), (dy-ty1));
|
*
|
* ////////////////////////////////////////////
|
* /
|
* /
|
* // convert TWD97->TM2
|
* /
|
* /
|
* ////////////////////////////////////////////
|
*
|
* dx = x2;
|
* dy = y2;
|
*
|
* toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx, &dy);
|
* // center longitude of taiwan is 121, for penghu is 119
|
*
|
* dx += 250000; // TM2 in Taiwan should add 250000
|
*
|
* printf("TWD97->TM2nTWD97 (%f, %f)nConvert (%.3f, %.3f)nOrigin (%.3f, %.3f)n", x2, y2, dx, dy, tx2, ty2);
|
* printf("Acuuracy (%.3f, X:%.3f, Y:%.3f)nn", sqrt((dx-tx2)*(dx-tx2)+(dy-ty2)*(dy-ty2)), (dx-tx2), (dy-ty2));
|
*
|
* ////////////////////////////////////////////
|
* /
|
* /
|
* // convert TM2->TWD67
|
* /
|
* /
|
* ////////////////////////////////////////////
|
*
|
* dx = tx1-250000; // should minus 250000 first in Taiwan
|
* dy = ty1;
|
*
|
* fromTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx, &dy);
|
*
|
* printf("TM2->TWD67nTM2 (%f, %f)nConvert (%.9f, %.9f)nOrigin (%.9f, %.9f)n", tx1, ty1, dx, dy, x1, y1);
|
* printf("Acuuracy (%.9f, X:%.9f, Y:%.9f)nn", sqrt((dx-x1)*(dx-x1)+(dy-y1)*(dy-y1)), (dx-x1), (dy-y1));
|
*
|
* ////////////////////////////////////////////
|
* /
|
* /
|
* // convert TM2->TWD97
|
* /
|
* /
|
* ////////////////////////////////////////////
|
*
|
* dx = tx2-250000; // should minus 250000 first in Taiwan
|
* dy = ty2;
|
*
|
* fromTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx, &dy);
|
*
|
* printf("TM2->TWD97nTM2 (%f, %f)nConvert (%.9f, %.9f)\nOrigin (%.9f, %.9f)\n", tx2, ty2, dx, dy, x2, y2);
|
* printf("Acuuracy (%.9f, X:%.9f, Y:%.9f)nn", sqrt((dx-x2)*(dx-x2)+(dy-y2)*(dy-y2)), (dx-x2), (dy-y2));
|
*
|
* ////////////////////////////////////////////
|
* /
|
* /
|
* // convert TWD67->TWD97
|
* /
|
* /
|
* ////////////////////////////////////////////
|
*
|
* dx = x1;
|
* dy = y1;
|
* dz = z1;
|
*
|
* toTWD97(&dx, &dy, &dz);
|
*
|
* dx1 = dx;
|
* dy1 = dy;
|
*
|
* toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx1, &dy1);
|
*
|
* dx1 += 250000; // TM2 in Taiwan should add 250000
|
*
|
* printf("TWD67->TWD97\nTWD67 (%.9f, %.9f, %6.2f) (%.3f, %.3f)n", x1, y1, z1, tx1, ty1);
|
* printf("Convert (%.9f, %.9f, %6.2f) (%.3f, %.3f)n", dx, dy, dz, dx1, dy1);
|
* printf("Origin (%.9f, %.9f, %6.2f) (%.3f, %.3f)n", x2, y2, z2, tx2, ty2);
|
* printf("Acuuracy (%.4f, X:%.4f, Y:%.4f)nn", sqrt((dx1-tx2)*(dx1-tx2)+(dy1-ty2)*(dy1-ty2)), (dx1-tx2), (dy1-ty2));
|
*
|
* ////////////////////////////////////////////
|
* /
|
* /
|
* // convert TWD97->TWD67
|
* /
|
* /
|
* ////////////////////////////////////////////
|
*
|
* dx = x2;
|
* dy = y2;
|
* dz = z2;
|
*
|
* toTWD67(&dx, &dy, &dz);
|
*
|
* dx1 = dx;
|
* dy1 = dy;
|
*
|
* toTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx1, &dy1);
|
*
|
* dx1 += 250000; // TM2 in Taiwan should add 250000
|
*
|
* printf("TWD97->TWD67nTWD97 (%.9f, %.9f, %6.2f) (%.3f, %.3f)n", x2, y2, z2, tx2, ty2);
|
* printf("Convert (%.9f, %.9f, %6.2f) (%.3f, %.3f)n", dx, dy, dz, dx1, dy1);
|
* printf("Origin (%.9f, %.9f, %6.2f) (%.3f, %.3f)n", x1, y1, z1, tx1, ty1);
|
* printf("Acuuracy (%.4f, X:%.4f, Y:%.4f)nn", sqrt((dx1-tx1)*(dx1-tx1)+(dy1-ty1)*(dy1-ty1)), (dx1-tx1), (dy1-ty1));
|
* }
|
*/
|
|
public static Coordinate fromTM2ToEPSG3826(Coordinate pt) {
|
Coordinate ptTWD67 = fromTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, pt.x - 250000.0, pt.y);
|
ptTWD67.z = 0;
|
Coordinate ptTWD97 = toTWD97(ptTWD67);
|
Coordinate pt97TM2 = toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, ptTWD97.x, ptTWD97.y);
|
pt97TM2.x += 250000.0;
|
pt97TM2.y -= 200.0;
|
return pt97TM2;
|
}
|
|
public static Coordinate fromTM2ToEPSG3825(Coordinate pt) {
|
Coordinate ptTWD67 = fromTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 191, TWD67_TM2, pt.x - 250000.0, pt.y);
|
ptTWD67.z = 0;
|
Coordinate ptTWD97 = toTWD97(ptTWD67);
|
Coordinate pt97TM2 = toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 191, TWD97_TM2, ptTWD97.x, ptTWD97.y);
|
pt97TM2.x += 250000.0;
|
pt97TM2.y -= 200.0;
|
return pt97TM2;
|
}
|
}
|
