package com.ximple.eofms.util; import com.vividsolutions.jts.geom.Coordinate; /** * TWDDatumConverter * User: Ulysses * Date: 2007/10/8 * Time: 下午 01:35:03 * 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)); * } */ /** * 由TM2座標轉換至TWD97座標 * @param pt TM2位置 * @return 新的TWD97座標 */ public static Coordinate fromTM2ToTWD97(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; return pt97TM2; } }