【未完成】GPS坐标转换为平面坐标

最近基于位置服务比较火，为了赶个时髦，我也尝试研究了一下位置数据。

想到的第一件事就是如果把GPS所代表的球面坐标转换为我们熟悉平面坐标。但是就这第一步的想法可就把我难坏了，因为理论上并不能把球面展开为平面，最后翻阅各种论文发现需要借助“高斯-克吕格坐标系转换算法”将球面局部近似展开为平面。

不得不提及一些老一辈的学者，本文主要参考杨启和教授的论文，他治学的严谨远远超过我的想象。当时计算机并不普及，公式的推导完全依靠人工，后面程序里将会展示公式有多么复杂。

代码如下：

package location.maputil;import java.text.DecimalFormat;import org.junit.Test;import org.slf4j.Logger;import org.slf4j.LoggerFactory;// 参考《GPS经纬度坐标转平面坐标的简化计算方法及精度分析》// 参考《WGS84、北京54、高斯-克吕格坐标系转换算法》// 参考《高斯-克吕格投影常系数和变系数正反解公式的讨论和应用》// 球面坐标理论上不能展开为平面坐标，高斯-克吕格坐标系是将球面坐标近似展开为多平面坐标系，本工具将每个城市作为一个独立的坐标系// 本工具目前只针对北京市、上海市// 北京：B~[39.4, 41.0] L~[115.4, 117.5] L0=116// 上海：B~[30.6,32] L~[120.7,122.5] L0=121.5// 重庆：B~[28.1,32.25] L~[105.1,110.2] L0=107.5// 需要其他城市的坐标系请与程序开发者联系public class Transformer {private static final Logger LOGGER = LoggerFactory.getLogger(Transformer.class);private static int Datum = 84; // 坐标系，默认为WGS84DecimalFormat dGPS = new DecimalFormat("0.000000");DecimalFormat dXY = new DecimalFormat("0");/* 参考椭球体基本量 */private double a; // 长半轴private double f; // 扁率 f=(a-b)/aprivate double b; // 短半轴 b=(1-f)aprivate double e; // 第一偏心率private double e1; // 第二偏心率private double L0; // 中央经度private double W0; // 原点维度private double k0; // 比例因子private double FE; // 东偏移private double FN; // 北偏移double PI = 3.141592653589794;double iPI = 0.0174532925199433; // PI/180public Transformer(String city) {LOGGER.info(city + " map is loading");// Datum 投影基准面类型：北京54基准面为54，西安80基准面为80，WGS84基准面为84if (Datum == 84) {a = 6378137;f = 1 / 298.257223563;} else if (Datum == 54) {a = 6378245;f = 1 / 298.3;} else if (Datum == 80) {a = 6378140;f = 1 / 298.257;} else {// 其他参数按照WGS84基准面为84处理a = 6378137;f = 1 / 298.257223563;}b = (1 - f) * a;e = Math.sqrt(2 * f - f * f);e1 = e / Math.sqrt(1 - e * e);L0 = 0; // 暂时处理W0 = 0; //k0 = 1; //FE = 0; //FN = 0; //if ("beijing".equals(city)) {L0 = 116;} else if ("shanghai".equals(city)) {L0 = 121.5;} else if ("chongqing".equals(city)) {L0 = 107.5;} else {LOGGER.info(city + " is not found, precision may be very low");}LOGGER.info(city + " map is loaded");}// 正解public String getXY(String lat, String lng) {double latitude = Double.valueOf(lat); // 纬度 B Wdouble longitude = Double.valueOf(lng); // 经度 L J/* 必要变量准备 */double B = (latitude - W0) * iPI; // 纬差弧度double L = (longitude - L0) * iPI; // 经差弧度double sinB = Math.sin(B);double cosB = Math.cos(B);double tanB = Math.tan(B);double N = a / Math.sqrt(1 - Math.pow(e * sinB, 2)); // 卯酉圈曲率半径double g = e1 * cosB;/* 求解参数s */double s; // 赤道至纬度latitude的经线弧长double B0;double B2;double B4;double B6;double B8;double C = Math.pow(a, 2) / b;B0 = 1 - 3.0 / 4.0 * Math.pow(e1, 2) + 45.0 / 64.0 * Math.pow(e1, 4)- 175.0 / 256.0 * Math.pow(e1, 6) + 11025.0 / 16384.0* Math.pow(e1, 8);B2 = B0 - 1;B4 = 15.0 / 32.0 * Math.pow(e1, 4) - 175.0 / 384.0 * Math.pow(e1, 6)+ 3675.0 / 8192.0 * Math.pow(e1, 8);B6 = 0 - 35.0 / 96.0 * Math.pow(e1, 6) + 735.0 / 2048.0* Math.pow(e1, 8);B8 = 315.0 / 1024.0 * Math.pow(e1, 8);s = C* (B0 * B + sinB* (B2 * cosB + B4 * Math.pow(cosB, 3) + B6* Math.pow(cosB, 5) + B8 * Math.pow(cosB, 7)));/* 求解平面直角坐标系坐标 */double xTemp = s+ Math.pow(L, 2)* N* sinB* cosB/ 2.0+ Math.pow(L, 4)* N* sinB* Math.pow(cosB, 3)* (5 - Math.pow(tanB, 2) + 9 * Math.pow(g, 2) + 4 * Math.pow(g,4)) / 24.0 + Math.pow(L, 6) * N * sinB* Math.pow(cosB, 5)* (61 - 58 * Math.pow(tanB, 2) + Math.pow(tanB, 4)) / 720.0;double yTemp = L* N* cosB+ Math.pow(L, 3)* N* Math.pow(cosB, 3)* (1 - Math.pow(tanB, 2) + Math.pow(g, 2))/ 6.0+ Math.pow(L, 5)* N* Math.pow(cosB, 5)* (5 - 18 * Math.pow(tanB, 2) + Math.pow(tanB, 4) + 14* Math.pow(g, 2) - 58 * Math.pow(g, 2)* Math.pow(tanB, 2)) / 120.0;double x = xTemp + FN;double y = yTemp + FE;return dXY.format(x) + "," + dXY.format(y);}// 反解public String getGPS(String xStr, String yStr) {double x = Double.valueOf(xStr);double y = Double.valueOf(yStr);double El1 = (1 - Math.sqrt(1 - Math.pow(e, 2)))/ (1 + Math.sqrt(1 - Math.pow(e, 2)));double Mf = (x - FN) / k0; // 真实坐标值double Q = Mf/ (a * (1 - Math.pow(e, 2) / 4.0 - 3 * Math.pow(e, 4) / 64.0 - 5 * Math.pow(e, 6) / 256.0));double Bf = Q + (3 * El1 / 2.0 - 27 * Math.pow(El1, 3) / 32.0)* Math.sin(2 * Q)+ (21 * Math.pow(El1, 2) / 16.0 - 55 * Math.pow(El1, 4) / 32.0)* Math.sin(4 * Q) + (151 * Math.pow(El1, 3) / 96.0)* Math.sin(6 * Q) + 1097 / 512.0 * Math.pow(El1, 4)* Math.sin(8 * Q);double sinBf = Math.sin(Bf);double tanBf = Math.tan(Bf);double cosBf = Math.cos(Bf);double Rf = a * (1 - Math.pow(e, 2))/ Math.sqrt(Math.pow(1 - Math.pow(e * sinBf, 2), 3));double Nf = a / Math.sqrt(1 - Math.pow(e * sinBf, 2)); // 卯酉圈曲率半径double Tf = Math.pow(tanBf, 2);double D = (y - FE) / (k0 * Nf);double Cf = Math.pow(e1, 2) * Math.pow(cosBf, 2);double B = Bf- Nf* tanBf/ Rf* (Math.pow(D, 2)/ 2.0- (5 + 3 * Tf + 10 * Cf - 9 * Tf * Cf - 4* Math.pow(Cf, 2) - 9 * Math.pow(e1, 2))* Math.pow(D, 4) / 24.0 + (61 + 90 * Tf + 45* Math.pow(Tf, 2) - 256 * Math.pow(e1, 2) - 3 * Math.pow(Cf, 2)) * Math.pow(D, 6) / 720.0);double L = L0* iPI+ 1/ cosBf* (D - (1 + 2 * Tf + Cf) * Math.pow(D, 3) / 6.0 + (5 - 2 * Cf+ 28 * Tf - 3 * Math.pow(Cf, 2) + 8 * Math.pow(e1, 2) + 24 * Math.pow(Tf, 2)) * Math.pow(D, 5) / 120.0);double Bangle = B / iPI;double Langle = L / iPI;double latitude = Bangle + W0; // 纬度 B W xdouble longitude = Langle; // 经度 L J yreturn dGPS.format(latitude) + "," + dGPS.format(longitude);}public static double distacne(String str1, String str2) {String[] strArr1 = str1.split(",");String[] strArr2 = str2.split(",");double x1 = Double.valueOf(strArr1[0]);double y1 = Double.valueOf(strArr1[1]);double x2 = Double.valueOf(strArr2[0]);double y2 = Double.valueOf(strArr2[1]);return Math.sqrt(Math.pow(x1 - x2, 2) + Math.pow(y1 - y2, 2));}@Testpublic void test() {Transformer transfer = new Transformer("beijing");// String str1 = transfer.getGPS("39.927938","116.338967");// String str2 = transfer.getGPS("39.926573","116.338581");// System.out.println(distacne(str1, str2));}public static void main(String[] args) {Transformer transfer = new Transformer("beijing");System.out.println(transfer.getXY("41", "120"));System.out.println(transfer.getGPS("4548288", "336579"));Transformer transfer1 = new Transformer("shanghai");System.out.println(transfer1.getXY("30.7", "122"));System.out.println(transfer1.getGPS("3397821", "47901"));Transformer transfer2 = new Transformer("chongqing");System.out.println(transfer2.getXY("29.490295", "106.486654"));System.out.println(transfer2.getXY("29.615467", "106.581515"));System.out.println(distacne("3359469,-1464555", "3372354,-1453338"));Transformer transfer3 = new Transformer("c");System.out.println(transfer3.getXY("10", "127"));Transformer transfer4 = new Transformer("beijing");String str1 = transfer4.getXY("39.927938", "116.338967");String str2 = transfer4.getXY("41", "116.338581");System.out.println(distacne(str1, str2));}}

之后会讲讲原理和推导，敬请期待

参考《GPS经纬度坐标转平面坐标的简化计算方法及精度分析》

参考《WGS84、北京54、高斯-克吕格坐标系转换算法》

参考《高斯-克吕格投影常系数和变系数正反解公式的讨论和应用》