|
|
@@ -0,0 +1,796 @@
|
|
|
+package com.skyversation.poiaddr.util;
|
|
|
+
|
|
|
+import org.opengis.referencing.operation.Projection;
|
|
|
+import org.osgeo.proj4j.*;
|
|
|
+
|
|
|
+import java.lang.reflect.InvocationTargetException;
|
|
|
+import java.util.ArrayList;
|
|
|
+import java.util.List;
|
|
|
+
|
|
|
+public class CoordTransform2 {
|
|
|
+
|
|
|
+ private final double x_PI = (Math.PI * 3000.0) / 180.0;
|
|
|
+ private final double ee = 0.00669342162296594323;
|
|
|
+ private final double a = 6378245.0;
|
|
|
+
|
|
|
+
|
|
|
+ private final double BD_FACTOR = (3.14159265358979324 * 3000.0) / 180.0;
|
|
|
+ private final double PI = Math.PI;
|
|
|
+ private final double RADIUS = 6378245.0;
|
|
|
+ private final double EE = 0.00669342162296594323;
|
|
|
+
|
|
|
+
|
|
|
+ // 定义常量
|
|
|
+ private final double A = 6378245.0;
|
|
|
+
|
|
|
+ public enum Projs {
|
|
|
+ WGS84, GCJ02, BD09, UTM4, SHCJ, METER, DEGREE
|
|
|
+ }
|
|
|
+
|
|
|
+ public enum EPSG {
|
|
|
+ WGS84("+proj=longlat +datum=WGS84 +no_defs"),
|
|
|
+ MERCATOR("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"),
|
|
|
+ GAUSS_KRUGER("+proj=tmerc +lat_0=0 +lon_0=120 +k=1 +x_0=40500000 +y_0=0 +ellps=GRS80 +units=m +no_defs");
|
|
|
+
|
|
|
+ private final String projString;
|
|
|
+
|
|
|
+ EPSG(String projString) {
|
|
|
+ this.projString = projString;
|
|
|
+ }
|
|
|
+
|
|
|
+ public String getProjString() {
|
|
|
+ return projString;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ public static boolean outOfSh(double lng, double lat){
|
|
|
+ lat += lat;
|
|
|
+ lng +=lng;
|
|
|
+ return !(lng > 115.487 && lng < 123.696 && lat > 28.260 && lat < 33.521);
|
|
|
+ }
|
|
|
+
|
|
|
+ private static CoordTransform2 instance = new CoordTransform2();
|
|
|
+ private CoordTransform2(){}
|
|
|
+ public static CoordTransform2 getInstance(){
|
|
|
+ if(instance == null) {
|
|
|
+ instance = new CoordTransform2();
|
|
|
+ }
|
|
|
+ return instance;
|
|
|
+ }
|
|
|
+
|
|
|
+
|
|
|
+ // 定义投影字典
|
|
|
+ private static final List<Projection> projs = new ArrayList<>();
|
|
|
+
|
|
|
+ // 定义坐标转换方法
|
|
|
+ public double[] wgs84ToShcj(double x, double y) {
|
|
|
+ if (outOfChina(x, y)) {
|
|
|
+ return new double[]{x, y};
|
|
|
+ } else {
|
|
|
+ double[] dlatLng = transformLatLng(x - 105.0, y - 35.0);
|
|
|
+ double radLat = (y / 180.0) * PI;
|
|
|
+ double magic = Math.sin(radLat);
|
|
|
+ magic = 1 - ee * magic * magic;
|
|
|
+ double sqrtMagic = Math.sqrt(magic);
|
|
|
+ dlatLng[0] = (dlatLng[0] * 180.0) / (((a * (1 - ee)) / (magic * sqrtMagic)) * PI);
|
|
|
+ dlatLng[1] = (dlatLng[1] * 180.0) / ((a / sqrtMagic) * Math.cos(radLat) * PI);
|
|
|
+ double mglat = y + dlatLng[0];
|
|
|
+ double mglng = x + dlatLng[1];
|
|
|
+ return new double[]{mglng, mglat};
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义坐标转换方法
|
|
|
+ public double[] gcj02ToWgs84(double x, double y) {
|
|
|
+ if (outOfChina(x, y)) {
|
|
|
+ return new double[]{x, y};
|
|
|
+ } else {
|
|
|
+ double[] dlatLng = transformLatLng(x - 105.0, y - 35.0);
|
|
|
+ double radLat = (y / 180.0) * PI;
|
|
|
+ double magic = Math.sin(radLat);
|
|
|
+ magic = 1 - ee * magic * magic;
|
|
|
+ double sqrtMagic = Math.sqrt(magic);
|
|
|
+ dlatLng[0] = (dlatLng[0] * 180.0) / (((a * (1 - ee)) / (magic * sqrtMagic)) * PI);
|
|
|
+ dlatLng[1] = (dlatLng[1] * 180.0) / ((a / sqrtMagic) * Math.cos(radLat) * PI);
|
|
|
+ double mglat = y + dlatLng[0];
|
|
|
+ double mglng = x + dlatLng[1];
|
|
|
+ return new double[]{x * 2 - mglng, y * 2 - mglat};
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义判断是否在中国境内的方法
|
|
|
+ public boolean outOfChina(double x, double y) {
|
|
|
+ // 纬度 3.86~53.55，经度 73.66~135.05
|
|
|
+ return!(x > 73.66 && x < 135.05 && y > 3.86 && y < 53.55);
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义坐标转换方法
|
|
|
+ public double[] transformLatLng(double x, double y) {
|
|
|
+ double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.sqrt(Math.abs(x));
|
|
|
+ ret += ((20.0 * Math.sin(6.0 * x * PI) + 20.0 * Math.sin(2.0 * x * PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((20.0 * Math.sin(y * PI) + 40.0 * Math.sin((y / 3.0) * PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((160.0 * Math.sin((y / 12.0) * PI) + 320 * Math.sin((y * PI) / 30.0)) * 2.0) / 3.0;
|
|
|
+ return new double[]{ret, ret};
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义经纬度转米的方法
|
|
|
+ public double[] degreeToMeter(double x, double y) {
|
|
|
+ double lon = (x * 20037508.34) / 180;
|
|
|
+ double lat = Math.log(Math.tan(((90 + y) * PI) / 360)) / (PI / 180);
|
|
|
+ lat = (lat * 20037508.34) / 180;
|
|
|
+ return new double[]{lon, lat};
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义米转经纬度的方法
|
|
|
+ public double[] meterToDegree(double x, double y) {
|
|
|
+ double lon = (x / 20037508.34) * 180;
|
|
|
+ double lat = (y / 20037508.34) * 180;
|
|
|
+ lat = (180 / PI) * (2 * Math.atan(Math.exp((lat * PI) / 180)) - PI / 2);
|
|
|
+ return new double[]{lon, lat};
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义 WGS84 坐标转上海城建的方法
|
|
|
+ public double[] wgs84ToShcj2(double x, double y) {
|
|
|
+ double[] xy = shcjGetUTMFromWGS(x, y);
|
|
|
+ return utmToShcj4(xy[0], xy[1]);
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义 UTM 转上海城建的方法
|
|
|
+ public double[] utmToShcj4(double x, double y) {
|
|
|
+ double DX = -500199.29965;
|
|
|
+ double DY = -3457078.805985;
|
|
|
+ double T = 0.0000001755;
|
|
|
+ double K = 1.0004000106;
|
|
|
+ return covertByFourParm(x, y, DX, DY, T, K);
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义上海城建坐标转 UTM 的方法
|
|
|
+ public double[] shcjGetUTMFromWGS(double lon, double lat) {
|
|
|
+ double a = 6378137;
|
|
|
+ double b = 6356752.3142451;
|
|
|
+ double f = (a - b) / a;
|
|
|
+
|
|
|
+ double eSquare = 2 * f - f * f;
|
|
|
+ double k0 = 0.9996;
|
|
|
+ double lonOrigin = 121.46714714;
|
|
|
+ double FN = 0;
|
|
|
+ // # 确保 longtitude 位于-180.00----179.9 之间
|
|
|
+ double lonTemp = lon + 180 - Math.floor((lon + 180) / 360) * 360 - 180;
|
|
|
+ double latRad = (lat * PI) / 180;
|
|
|
+ double lonRad = (lonTemp * PI) / 180;
|
|
|
+ double lonOriginRad = (lonOrigin * PI) / 180;
|
|
|
+ double e2Square = eSquare / (1 - eSquare);
|
|
|
+
|
|
|
+ double V = a / Math.sqrt(1 - eSquare * Math.pow(Math.sin(latRad), 2));
|
|
|
+ double T = Math.pow(Math.tan(latRad), 2);
|
|
|
+ double C = e2Square * Math.pow(Math.cos(latRad), 2);
|
|
|
+ double A = Math.cos(latRad) * (lonRad - lonOriginRad);
|
|
|
+ double M = a * ((1 - eSquare / 4 - (3 * Math.pow(eSquare, 2)) / 64 - (5 * Math.pow(eSquare, 3)) / 256) * latRad - ((3 * eSquare) / 8 + (3 * Math.pow(eSquare, 2)) / 32 + (45 * Math.pow(eSquare, 3)) / 1024) * Math.sin(2 * latRad) + ((15 * Math.pow(eSquare, 2)) / 256 + (45 * Math.pow(eSquare, 3)) / 1024) * Math.sin(4 * latRad) - ((35 * Math.pow(eSquare, 3)) / 3072) * Math.sin(6 * latRad));
|
|
|
+
|
|
|
+ // # x
|
|
|
+ double UTMEasting = k0 * V * (A + ((1 - T + C) * Math.pow(A, 3)) / 6 + ((5 - 18 * T + Math.pow(T, 2) + 72 * C - 58 * e2Square) * Math.pow(A, 5)) / 120) + 500000.0;
|
|
|
+ // # y
|
|
|
+ double UTMNorthing = k0 * (M + V * Math.tan(latRad) * (Math.pow(A, 2) / 2 + ((5 - T + 9 * C + 4 * Math.pow(C, 2)) * Math.pow(A, 4)) / 24 + ((61 - 58 * T + Math.pow(T, 2) + 600 * C - 330 * e2Square) * Math.pow(A, 6)) / 720));
|
|
|
+ //# 南半球纬度起点为 10000000.0m
|
|
|
+ UTMNorthing += FN;
|
|
|
+ double[] xy = new double[2];
|
|
|
+ xy[0] = UTMEasting;
|
|
|
+ xy[1] = UTMNorthing;
|
|
|
+ return xy;
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义上海城建坐标转 WGS84 的方法
|
|
|
+ public double[] shcjToWgs84(double x, double y) {
|
|
|
+ double[] xy = shcjToUtm4(x, y);
|
|
|
+ return shcjGetWGSFromUTM(xy[0], xy[1]);
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义上海城建坐标转 UTM 的方法
|
|
|
+ public double[] shcjToUtm4(double x, double y) {
|
|
|
+ double DX = 499999.90104;
|
|
|
+ double DY = 3455696.403019;
|
|
|
+ double T = -0.0000001755;
|
|
|
+ double K = 0.999600149344;
|
|
|
+ return covertByFourParm(x, y, DX, DY, T, K);
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义四参数公式
|
|
|
+ public double[] covertByFourParm(double x, double y, double dx, double dy, double a, double k) {
|
|
|
+ double px = 0;
|
|
|
+ double py = 0;
|
|
|
+ px = x * k * Math.cos(a) - y * k * Math.sin(a) + dx;
|
|
|
+ py = x * k * Math.sin(a) + y * k * Math.cos(a) + dy;
|
|
|
+ double[] xy = new double[2];
|
|
|
+ xy[0] = px;
|
|
|
+ xy[1] = py;
|
|
|
+ return xy;
|
|
|
+ }
|
|
|
+
|
|
|
+ // 定义上海城建坐标转 WGS84 的方法
|
|
|
+ public double[] shcjGetWGSFromUTM(double x, double y) {
|
|
|
+ // WGS84
|
|
|
+ double a = 6378137; // 椭球体长半轴
|
|
|
+ double b = 6356752.3142451; // 椭球体短半轴
|
|
|
+ x = 500000 - x;
|
|
|
+ double k0 = 0.9996;
|
|
|
+ double e = Math.sqrt(1 - Math.pow(b, 2) / Math.pow(a, 2));
|
|
|
+ // # calculate the meridional arc
|
|
|
+ double M = y / k0;
|
|
|
+ //# calculate footprint latitude
|
|
|
+ double mu = M / (a * (1 - Math.pow(e, 2) / 4 - (3 * Math.pow(e, 4)) / 64 - (5 * Math.pow(e, 6)) / 256));
|
|
|
+ double e1 = (1 - Math.pow(1 - Math.pow(e, 2), 1.0 / 2)) / (1 + Math.pow(1 - Math.pow(e, 2), 1.0 / 2));
|
|
|
+ double J1 = (3 * e1) / 2 - (27 * Math.pow(e1, 3)) / 32;
|
|
|
+ double J2 = (21 * Math.pow(e1, 2)) / 16 - (55 * Math.pow(e1, 4)) / 32;
|
|
|
+ double J3 = (151 * Math.pow(e1, 3)) / 96;
|
|
|
+ double J4 = (1097 * Math.pow(e1, 4)) / 512;
|
|
|
+ double fp = mu + J1 * Math.sin(2 * mu) + J2 * Math.sin(4 * mu) + J3 * Math.sin(6 * mu) + J4 * Math.sin(8 * mu);
|
|
|
+ // # Calculate Latitude and Longitude
|
|
|
+ double e2 = Math.pow(e, 2) / (1 - Math.pow(e, 2));
|
|
|
+ double C1 = e2 * Math.pow(Math.cos(fp), 2);
|
|
|
+ double T1 = Math.pow(Math.tan(fp), 2);
|
|
|
+ double R1 = (a * (1 - Math.pow(e, 2))) / Math.pow(1 - Math.pow(e * Math.sin(fp), 2), 3.0 / 2);
|
|
|
+ //# This is the same as rho in the forward conversion formulas above, but calculated for fp instead of lat.
|
|
|
+ double N1 = a / Math.pow(1 - Math.pow(e * Math.sin(fp), 2), 1.0 / 2);
|
|
|
+ //# This is the same as nu in the forward conversion formulas above, but calculated for fp instead of lat.
|
|
|
+ double D = x / (N1 * k0);
|
|
|
+ double Q1 = (N1 * Math.tan(fp)) / R1;
|
|
|
+ double Q2 = Math.pow(D, 2) / 2;
|
|
|
+ double Q3 = ((5 + 3 * T1 + 10 * C1 - 4 * Math.pow(C1, 2) - 9 * e2) * Math.pow(D, 4)) / 24;
|
|
|
+ double Q4 = ((61 + 90 * T1 + 298 * C1 + 45 * Math.pow(T1, 2) - 3 * Math.pow(C1, 2) - 252 * e2) * Math.pow(D, 6)) / 720;
|
|
|
+ double lat = ((fp - Q1 * (Q2 - Q3 + Q4)) * 180) / PI;
|
|
|
+ // System.out.println("lat===="+Math.toRadians(fp - Q1*(Q2 - Q3 + Q4)));
|
|
|
+ double Q5 = D;
|
|
|
+ double Q6 = ((1 + 2 * T1 + C1) * Math.pow(D, 3)) / 6;
|
|
|
+ double Q7 = ((5 - 2 * C1 + 28 * T1 - 3 * Math.pow(C1, 2) + 8 * e2 + 24 * Math.pow(T1, 2)) * Math.pow(D, 5)) / 120;
|
|
|
+ double lonMid = 121.46714714;
|
|
|
+ double lon = lonMid - (((Q5 - Q6 + Q7) / Math.cos(fp)) * 180) / PI;
|
|
|
+ double[] xy = new double[2];
|
|
|
+ xy[0] = lon;
|
|
|
+ xy[1] = lat;
|
|
|
+ return xy;
|
|
|
+ }
|
|
|
+
|
|
|
+ /***
|
|
|
+ * 计算两个经纬度之间的距离，返回单位米
|
|
|
+ * @param lat1
|
|
|
+ * @param lon1
|
|
|
+ * @param lat2
|
|
|
+ * @param lon2
|
|
|
+ * @return
|
|
|
+ */
|
|
|
+ public double calculateDistance(double lat1, double lon1, double lat2, double lon2) {
|
|
|
+ // 将经纬度转换为弧度
|
|
|
+ double lat1Rad = Math.toRadians(lat1);
|
|
|
+ double lon1Rad = Math.toRadians(lon1);
|
|
|
+ double lat2Rad = Math.toRadians(lat2);
|
|
|
+ double lon2Rad = Math.toRadians(lon2);
|
|
|
+
|
|
|
+// // 地球半径（单位：千米）
|
|
|
+// double earthRadius = 6371.0;
|
|
|
+
|
|
|
+ // Haversine公式计算两点之间的距离
|
|
|
+ double dLat = lat2Rad - lat1Rad;
|
|
|
+ double dLon = lon2Rad - lon1Rad;
|
|
|
+ double a = Math.pow(Math.sin(dLat / 2), 2) + Math.cos(lat1Rad) * Math.cos(lat2Rad) * Math.pow(Math.sin(dLon / 2), 2);
|
|
|
+ double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
|
|
|
+ double distance = RADIUS * c;
|
|
|
+
|
|
|
+ return distance;
|
|
|
+ }
|
|
|
+
|
|
|
+ public boolean checkLatlon(double lat, double lon){
|
|
|
+ if(lat >85 || lat < -85|| lon >180 || lon < -180){
|
|
|
+ return false;
|
|
|
+ }
|
|
|
+ return true;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] bd09_to_gcj02(double bd_lon, double bd_lat) {
|
|
|
+ double x = bd_lon - 0.0065;
|
|
|
+ double y = bd_lat - 0.006;
|
|
|
+ double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * x_PI);
|
|
|
+ double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * x_PI);
|
|
|
+ double gg_lng = z * Math.cos(theta);
|
|
|
+ double gg_lat = z * Math.sin(theta);
|
|
|
+ return new double[]{gg_lng, gg_lat};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] gcj02_to_bd09(double lng, double lat) {
|
|
|
+ double z = Math.sqrt(lng * lng + lat * lat) + 0.00002 * Math.sin(lat * x_PI);
|
|
|
+ double theta = Math.atan2(lat, lng) + 0.000003 * Math.cos(lng * x_PI);
|
|
|
+ double bd_lng = z * Math.cos(theta) + 0.0065;
|
|
|
+ double bd_lat = z * Math.sin(theta) + 0.006;
|
|
|
+ return new double[]{bd_lng, bd_lat};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] wgs84_to_gcj02(double lng, double lat) {
|
|
|
+ if (out_of_china(lng, lat)) {
|
|
|
+ return new double[]{lng, lat};
|
|
|
+ } else {
|
|
|
+ double dlat = transformlat(lng - 105.0, lat - 35.0);
|
|
|
+ double dlng = transformlng(lng - 105.0, lat - 35.0);
|
|
|
+ double radlat = (lat / 180.0) * PI;
|
|
|
+ double magic = Math.sin(radlat);
|
|
|
+ magic = 1 - ee * magic * magic;
|
|
|
+ double sqrtmagic = Math.sqrt(magic);
|
|
|
+ dlat = (dlat * 180.0) / (((a * (1 - ee)) / (magic * sqrtmagic)) * PI);
|
|
|
+ dlng = (dlng * 180.0) / ((a / sqrtmagic) * Math.cos(radlat) * PI);
|
|
|
+ double mglat = lat + dlat;
|
|
|
+ double mglng = lng + dlng;
|
|
|
+ return new double[]{mglng, mglat};
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+// public double[] gcj02_to_wgs84(double lng, double lat) {
|
|
|
+// if (out_of_china(lng, lat)) {
|
|
|
+// return new double[]{lng, lat};
|
|
|
+// } else {
|
|
|
+// double dlat = transformlat(lng - 105.0, lat - 35.0);
|
|
|
+// double dlng = transformlng(lng - 105.0, lat - 35.0);
|
|
|
+// double radlat = (lat / 180.0) * PI;
|
|
|
+// double magic = Math.sin(radlat);
|
|
|
+// magic = 1 - ee * magic * magic;
|
|
|
+// double sqrtmagic = Math.sqrt(magic);
|
|
|
+// dlat = (dlat * 180.0) / (((a * (1 - ee)) / (magic * sqrtmagic)) * PI);
|
|
|
+// dlng = (dlng * 180.0) / ((a / sqrtmagic) * Math.cos(radlat) * PI);
|
|
|
+// double mglat = lat + dlat;
|
|
|
+// double mglng = lng + dlng;
|
|
|
+// return new double[]{lng * 2 - mglng, lat * 2 - mglat};
|
|
|
+// }
|
|
|
+// }
|
|
|
+
|
|
|
+// public boolean out_of_china(double lng, double lat) {
|
|
|
+// return !(lng > 73.66 && lng < 135.05 && lat > 3.86 && lat < 53.55);
|
|
|
+// }
|
|
|
+
|
|
|
+// public double transformlat(double lng, double lat) {
|
|
|
+// double ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * Math.sqrt(Math.abs(lng));
|
|
|
+// ret += ((20.0 * Math.sin(6.0 * lng * PI) + 20.0 * Math.sin(2.0 * lng * PI)) * 2.0) / 3.0;
|
|
|
+// ret += ((20.0 * Math.sin(lat * PI) + 40.0 * Math.sin((lat / 3.0) * PI)) * 2.0) / 3.0;
|
|
|
+// ret += ((160.0 * Math.sin((lat / 12.0) * PI) + 320 * Math.sin((lat * PI) / 30.0)) * 2.0) / 3.0;
|
|
|
+// return ret;
|
|
|
+// }
|
|
|
+//
|
|
|
+// public double transformlng(double lng, double lat) {
|
|
|
+// double ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * Math.sqrt(Math.abs(lng));
|
|
|
+// ret += ((20.0 * Math.sin(6.0 * lng * PI) + 20.0 * Math.sin(2.0 * lng * PI)) * 2.0) / 3.0;
|
|
|
+// ret += ((20.0 * Math.sin(lng * PI) + 40.0 * Math.sin((lng / 3.0) * PI)) * 2.0) / 3.0;
|
|
|
+// ret += ((150.0 * Math.sin((lng / 12.0) * PI) + 300.0 * Math.sin((lng / 30.0) * PI)) * 2.0) / 3.0;
|
|
|
+// return ret;
|
|
|
+// }
|
|
|
+
|
|
|
+ public double[] degree_to_meter(double lng, double lat) {
|
|
|
+ double x = (lng * 20037508.34) / 180;
|
|
|
+ double y = Math.log(Math.tan(((90 + lat) * PI) / 360)) / (PI / 180);
|
|
|
+ y = (y * 20037508.34) / 180;
|
|
|
+ return new double[]{x, y};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] meter_to_degree(double x, double y) {
|
|
|
+ double lon = (x / 20037508.34) * 180;
|
|
|
+ double lat = (y / 20037508.34) * 180;
|
|
|
+ lat = (180 / PI) * (2 * Math.atan(Math.exp((lat * PI) / 180)) - PI / 2);
|
|
|
+ return new double[]{lon, lat};
|
|
|
+ }
|
|
|
+
|
|
|
+// public double[] wgs84_to_shcj(double x, double y) {
|
|
|
+// double[] xy = new double[2];
|
|
|
+// xy = shcj_get_UTM_from_WGS(x, y);
|
|
|
+// return utm_to_shcj4(xy[0], xy[1]);
|
|
|
+// }
|
|
|
+
|
|
|
+ public double[] utm_to_shcj4(double x, double y) {
|
|
|
+ double DX = -500199.29965;
|
|
|
+ double DY = -3457078.805985;
|
|
|
+ double T = 0.0000001755;
|
|
|
+ double K = 1.0004000106;
|
|
|
+ return covert_by_four_parm(x, y, DX, DY, T, K);
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] shcj_get_UTM_from_WGS(double lon, double lat) {
|
|
|
+ double a = 6378137;
|
|
|
+ double b = 6356752.3142451;
|
|
|
+ double f = (a - b) / a;
|
|
|
+ double eSquare = 2 * f - f * f;
|
|
|
+ double k0 = 0.9996;
|
|
|
+ double lonOrigin = 121.46714714;
|
|
|
+ double FN = 0;
|
|
|
+ double lonTemp = lon + 180 - Math.floor((lon + 180) / 360) * 360 - 180;
|
|
|
+ double latRad = (lat * PI) / 180;
|
|
|
+ double lonRad = (lonTemp * PI) / 180;
|
|
|
+ double lonOriginRad = (lonOrigin * PI) / 180;
|
|
|
+ double e2Square = eSquare / (1 - eSquare);
|
|
|
+ double V = a / Math.sqrt(1 - eSquare * Math.pow(Math.sin(latRad), 2));
|
|
|
+ double T = Math.pow(Math.tan(latRad), 2);
|
|
|
+ double C = e2Square * Math.pow(Math.cos(latRad), 2);
|
|
|
+ double A = Math.cos(latRad) * (lonRad - lonOriginRad);
|
|
|
+ double M = a * ((1 - eSquare / 4 - (3 * Math.pow(eSquare, 2)) / 64 - (5 * Math.pow(eSquare, 3)) / 256) * latRad -
|
|
|
+ ((3 * eSquare) / 8 + (3 * Math.pow(eSquare, 2)) / 32 + (45 * Math.pow(eSquare, 3)) / 1024) *
|
|
|
+ Math.sin(2 * latRad) +
|
|
|
+ ((15 * Math.pow(eSquare, 2)) / 256 + (45 * Math.pow(eSquare, 3)) / 1024) *
|
|
|
+ Math.sin(4 * latRad) -
|
|
|
+ ((35 * Math.pow(eSquare, 3)) / 3072) * Math.sin(6 * latRad));
|
|
|
+ double UTMEasting = k0 * V * (A + ((1 - T + C) * Math.pow(A, 3)) / 6 +
|
|
|
+ ((5 - 18 * T + Math.pow(T, 2) + 72 * C - 58 * e2Square) * Math.pow(A, 5)) / 120) + 500000.0;
|
|
|
+ double UTMNorthing = k0 * (M + V * Math.tan(latRad) * (Math.pow(A, 2) / 2 +
|
|
|
+ ((5 - T + 9 * C + 4 * Math.pow(C, 2)) * Math.pow(A, 4)) / 24 +
|
|
|
+ ((61 - 58 * T + Math.pow(T, 2) + 600 * C - 330 * e2Square) * Math.pow(A, 6)) / 720));
|
|
|
+ UTMNorthing += FN;
|
|
|
+ return new double[]{UTMEasting, UTMNorthing};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] shcj_to_wgs84(double x, double y) {
|
|
|
+ double[] xy = new double[2];
|
|
|
+ xy = shcj_to_utm4(x, y);
|
|
|
+ return shcj_get_WGS_from_UTM(xy[0], xy[1]);
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] shcj_to_utm4(double x, double y) {
|
|
|
+ double DX = 499999.90104;
|
|
|
+ double DY = 3455696.403019;
|
|
|
+ double T = -0.0000001755;
|
|
|
+ double K = 0.999600149344;
|
|
|
+ return covert_by_four_parm(x, y, DX, DY, T, K);
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] covert_by_four_parm(double x, double y, double dx, double dy, double a, double k) {
|
|
|
+ double px = x * k * Math.cos(a) - y * k * Math.sin(a) + dx;
|
|
|
+ double py = x * k * Math.sin(a) + y * k * Math.cos(a) + dy;
|
|
|
+ return new double[]{px, py};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] shcj_get_WGS_from_UTM(double x, double y) {
|
|
|
+ double a = 6378137;
|
|
|
+ double b = 6356752.3142451;
|
|
|
+ x = 500000 - x;
|
|
|
+ double k0 = 0.9996;
|
|
|
+ double e = Math.sqrt(1 - Math.pow(b, 2) / Math.pow(a, 2));
|
|
|
+ double M = y / k0;
|
|
|
+ double mu = M / (a * (1 - Math.pow(e, 2) / 4 - (3 * Math.pow(e, 4)) / 64 - (5 * Math.pow(e, 6)) / 256));
|
|
|
+ double e1 = (1 - Math.pow(1 - Math.pow(e, 2), 1.0 / 2)) / (1 + Math.pow(1 - Math.pow(e, 2), 1.0 / 2));
|
|
|
+ double J1 = (3 * e1) / 2 - (27 * Math.pow(e1, 3)) / 32;
|
|
|
+ double J2 = (21 * Math.pow(e1, 2)) / 16 - (55 * Math.pow(e1, 4)) / 32;
|
|
|
+ double J3 = (151 * Math.pow(e1, 3)) / 96;
|
|
|
+ double J4 = (1097 * Math.pow(e1, 4)) / 512;
|
|
|
+ double fp = mu + J1 * Math.sin(2 * mu) + J2 * Math.sin(4 * mu) + J3 * Math.sin(6 * mu) + J4 * Math.sin(8 * mu);
|
|
|
+ double e2 = Math.pow(e, 2) / (1 - Math.pow(e, 2));
|
|
|
+ double C1 = e2 * Math.pow(Math.cos(fp), 2);
|
|
|
+ double T1 = Math.pow(Math.tan(fp), 2);
|
|
|
+ double R1 = (a * (1 - Math.pow(e, 2))) / Math.pow(1 - Math.pow(e * Math.sin(fp), 2), 3.0 / 2);
|
|
|
+ double N1 = a / Math.pow(1 - Math.pow(e * Math.sin(fp), 2), 1.0 / 2);
|
|
|
+ double D = x / (N1 * k0);
|
|
|
+ double Q1 = (N1 * Math.tan(fp)) / R1;
|
|
|
+ double Q2 = Math.pow(D, 2) / 2;
|
|
|
+ double Q3 = ((5 + 3 * T1 + 10 * C1 - 4 * Math.pow(C1, 2) - 9 * e2) * Math.pow(D, 4)) / 24;
|
|
|
+ double Q4 = ((61 + 90 * T1 + 298 * C1 + 45 * Math.pow(T1, 2) - 3 * Math.pow(C1, 2) - 252 * e2) * Math.pow(D, 6)) / 720;
|
|
|
+ double lat = ((fp - Q1 * (Q2 - Q3 + Q4)) * 180) / PI;
|
|
|
+ double Q5 = D;
|
|
|
+ double Q6 = ((1 + 2 * T1 + C1) * Math.pow(D, 3)) / 6;
|
|
|
+ double Q7 = ((5 - 2 * C1 + 28 * T1 - 3 * Math.pow(C1, 2) + 8 * e2 + 24 * Math.pow(T1, 2)) * Math.pow(D, 5)) / 120;
|
|
|
+ double lonmid = 121.46714714;
|
|
|
+ double lon = lonmid - (((Q5 - Q6 + Q7) / Math.cos(fp)) * 180) / PI;
|
|
|
+ return new double[]{lon, lat};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] wgs84_to_bd09(double x, double y) {
|
|
|
+ double[] ll = wgs84_to_gcj02(x, y);
|
|
|
+ ll = gcj02_to_bd09(ll[0], ll[1]);
|
|
|
+ return ll;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] bd09_to_wgs84(double x, double y) {
|
|
|
+ double[] ll = bd09_to_gcj02(x, y);
|
|
|
+ ll = gcj02_to_wgs84(ll[0], ll[1]);
|
|
|
+ return ll;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] gcj02_to_shcj(double x, double y) {
|
|
|
+ double[] ll = gcj02_to_wgs84(x, y);
|
|
|
+ ll = wgs84_to_shcj(ll[0], ll[1]);
|
|
|
+ return ll;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] shcj_to_gcj02(double x, double y) {
|
|
|
+ double[] ll = shcj_to_wgs84(x, y);
|
|
|
+ ll = wgs84_to_gcj02(ll[0], ll[1]);
|
|
|
+ return ll;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] convert_proj_from_to(Projs from, Projs to, double[] xy) {
|
|
|
+ if (from == to) {
|
|
|
+ return xy;
|
|
|
+ } else {
|
|
|
+ String fromProj = from.toString().toLowerCase();
|
|
|
+ String toProj = to.toString().toLowerCase();
|
|
|
+ String targetMethod = fromProj + "_to_" + toProj;
|
|
|
+ try {
|
|
|
+ return (double[]) this.getClass().getMethod(targetMethod, double.class, double.class).invoke(this, xy[0], xy[1]);
|
|
|
+ } catch (IllegalAccessException e) {
|
|
|
+ throw new RuntimeException(e);
|
|
|
+ } catch (InvocationTargetException e) {
|
|
|
+ throw new RuntimeException(e);
|
|
|
+ } catch (NoSuchMethodException e) {
|
|
|
+ throw new RuntimeException(e);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] WGS84ToSH2000(double lng, double lat) {
|
|
|
+ if (out_of_sh(lng, lat)) {
|
|
|
+ return new double[]{lng, lat};
|
|
|
+ } else {
|
|
|
+ double[] templatlng = wgs84_to_shcj(lng, lat);
|
|
|
+ templatlng = meter_to_degree(templatlng[0], templatlng[1]);
|
|
|
+ return new double[]{templatlng[0], templatlng[1]};
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] SH2000ToWGS84(double x, double y) {
|
|
|
+ double[] templatlng = meter_to_degree(x, y);
|
|
|
+ if (out_of_sh2000(templatlng[0], templatlng[1])) {
|
|
|
+ return new double[]{x, y};
|
|
|
+ } else {
|
|
|
+ templatlng = shcj_to_wgs84(x, y);
|
|
|
+ return templatlng;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] SH2000lnglatToWGS84(double lng, double lat) {
|
|
|
+ if (out_of_sh2000(lng, lat)) {
|
|
|
+ return new double[]{lng, lat};
|
|
|
+ } else {
|
|
|
+ double[] templatlng = degree_to_meter(lng, lat);
|
|
|
+ templatlng = shcj_to_wgs84(templatlng[0], templatlng[1]);
|
|
|
+ return templatlng;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ private boolean out_of_sh(double lng, double lat) {
|
|
|
+ // 实现out_of_sh方法
|
|
|
+ return false;
|
|
|
+ }
|
|
|
+
|
|
|
+ private boolean out_of_sh2000(double lng, double lat) {
|
|
|
+ // 实现out_of_sh2000方法
|
|
|
+ return false;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] wgs84_to_shcj(double x, double y) {
|
|
|
+ double[] xy = new double[2];
|
|
|
+ xy = shcj_get_UTM_from_WGS(x, y);
|
|
|
+ return utm_to_shcj4(xy[0], xy[1]);
|
|
|
+ }
|
|
|
+
|
|
|
+
|
|
|
+ public double[] BD09ToGCJ02(double lng, double lat) {
|
|
|
+ double x = lng - 0.0065;
|
|
|
+ double y = lat - 0.006;
|
|
|
+ double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * BD_FACTOR);
|
|
|
+ double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * BD_FACTOR);
|
|
|
+ double gg_lng = z * Math.cos(theta);
|
|
|
+ double gg_lat = z * Math.sin(theta);
|
|
|
+ return new double[]{gg_lng, gg_lat};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] GCJ02ToBD09(double lng, double lat) {
|
|
|
+ lat = +lat;
|
|
|
+ lng = +lng;
|
|
|
+ double z = Math.sqrt(lng * lng + lat * lat) + 0.00002 * Math.sin(lat * BD_FACTOR);
|
|
|
+ double theta = Math.atan2(lat, lng) + 0.000003 * Math.cos(lng * BD_FACTOR);
|
|
|
+ double bd_lng = z * Math.cos(theta) + 0.0065;
|
|
|
+ double bd_lat = z * Math.sin(theta) + 0.006;
|
|
|
+ return new double[]{bd_lng, bd_lat};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] WGS84ToBD09(double lng, double lat) {
|
|
|
+ double[] GCJ02latlng = WGS84ToGCJ02(lng, lat);
|
|
|
+ double[] BD09latlng = GCJ02ToBD09(GCJ02latlng[0], GCJ02latlng[1]);
|
|
|
+ return BD09latlng;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] BD09ToWGS84(double lng, double lat) {
|
|
|
+ double[] GCJ02latlng = BD09ToGCJ02(lng, lat);
|
|
|
+ double[] wgslatlng = GCJ02ToWGS84(GCJ02latlng[0], GCJ02latlng[1]);
|
|
|
+ return wgslatlng;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] WGS84ToGCJ02(double lng, double lat) {
|
|
|
+ lat = +lat;
|
|
|
+ lng = +lng;
|
|
|
+ if (out_of_china(lng, lat)) {
|
|
|
+ return new double[]{lng, lat};
|
|
|
+ } else {
|
|
|
+ double[] d = delta(lng, lat);
|
|
|
+ return new double[]{lng + d[0], lat + d[1]};
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] GCJ02ToWGS84(double lng, double lat) {
|
|
|
+ lat = +lat;
|
|
|
+ lng = +lng;
|
|
|
+ if (out_of_china(lng, lat)) {
|
|
|
+ return new double[]{lng, lat};
|
|
|
+ } else {
|
|
|
+ double[] d = delta(lng, lat);
|
|
|
+ double mgLng = lng + d[0];
|
|
|
+ double mgLat = lat + d[1];
|
|
|
+ return new double[]{lng * 2 - mgLng, lat * 2 - mgLat};
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ private double[] delta(double lng, double lat) {
|
|
|
+ double dLng = transformLng(lng - 105, lat - 35);
|
|
|
+ double dLat = transformLat(lng - 105, lat - 35);
|
|
|
+ double radLat = (lat / 180) * PI;
|
|
|
+ double magic = Math.sin(radLat);
|
|
|
+ magic = 1 - EE * magic * magic;
|
|
|
+ double sqrtMagic = Math.sqrt(magic);
|
|
|
+ dLng = (dLng * 180) / ((RADIUS / sqrtMagic) * Math.cos(radLat) * PI);
|
|
|
+ dLat = (dLat * 180) / (((RADIUS * (1 - EE)) / (magic * sqrtMagic)) * PI);
|
|
|
+ return new double[]{dLng, dLat};
|
|
|
+ }
|
|
|
+
|
|
|
+ private double transformLng(double lng, double lat) {
|
|
|
+ lat = +lat;
|
|
|
+ lng = +lng;
|
|
|
+ double ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * Math.sqrt(Math.abs(lng));
|
|
|
+ ret += ((20.0 * Math.sin(6.0 * lng * PI) + 20.0 * Math.sin(2.0 * lng * PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((20.0 * Math.sin(lng * PI) + 40.0 * Math.sin((lng / 3.0) * PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((150.0 * Math.sin((lng / 12.0) * PI) + 300.0 * Math.sin((lng / 30.0) * PI)) * 2.0) / 3.0;
|
|
|
+ return ret;
|
|
|
+ }
|
|
|
+
|
|
|
+ private double transformLat(double lng, double lat) {
|
|
|
+ lat = +lat;
|
|
|
+ lng = +lng;
|
|
|
+ double ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * Math.sqrt(Math.abs(lng));
|
|
|
+ ret += ((20.0 * Math.sin(6.0 * lng * PI) + 20.0 * Math.sin(2.0 * lng * PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((20.0 * Math.sin(lat * PI) + 40.0 * Math.sin((lat / 3.0) * PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((160.0 * Math.sin((lat / 12.0) * PI) + 320 * Math.sin((lat * PI) / 30.0)) * 2.0) / 3.0;
|
|
|
+ return ret;
|
|
|
+ }
|
|
|
+
|
|
|
+ private boolean out_of_china(double lng, double lat) {
|
|
|
+ lat = +lat;
|
|
|
+ lng = +lng;
|
|
|
+ return !(lng > 73.66 && lng < 135.05 && lat > 3.86 && lat < 53.55);
|
|
|
+ }
|
|
|
+
|
|
|
+ private CRSFactory crsFactory = new CRSFactory();
|
|
|
+
|
|
|
+ // 定义84坐标系
|
|
|
+ private String wgs84Str = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
|
|
|
+// private CoordinateReferenceSystem wgs84 = crsFactory.createFromParameters(
|
|
|
+// "WGS84", wgs84Str);
|
|
|
+ private CoordinateReferenceSystem wgs84 =
|
|
|
+ crsFactory.createFromName("EPSG:4326");
|
|
|
+
|
|
|
+ // 定义CGCS2000坐标系
|
|
|
+ private String cgcs2000Str = "+proj=utm +zone=50 +ellps=GRS80 +units=m +no_defs";
|
|
|
+// private CoordinateReferenceSystem cgcs2000 = crsFactory.createFromParameters(
|
|
|
+// "CGCS2000", cgcs2000Str);
|
|
|
+ private CoordinateReferenceSystem cgcs2000 =
|
|
|
+ crsFactory.createFromName("EPSG:4490");
|
|
|
+
|
|
|
+ public double[] gcj02_to_wgs84(double lng, double lat) {
|
|
|
+ if (out_of_china(lng, lat)) {
|
|
|
+ return new double[]{lng, lat};
|
|
|
+ } else {
|
|
|
+ double dlat = transformlat(lng - 105.0, lat - 35.0);
|
|
|
+ double dlng = transformlng(lng - 105.0, lat - 35.0);
|
|
|
+ double radlat = (lat / 180.0) * Math.PI;
|
|
|
+ double magic = Math.sin(radlat);
|
|
|
+ magic = 1 - EE * magic * magic;
|
|
|
+ double sqrtmagic = Math.sqrt(magic);
|
|
|
+ dlat = (dlat * 180.0) / (((A * (1 - EE)) / (magic * sqrtmagic)) * Math.PI);
|
|
|
+ dlng = (dlng * 180.0) / ((A / sqrtmagic) * Math.cos(radlat) * Math.PI);
|
|
|
+ double mglat = lat + dlat;
|
|
|
+ double mglng = lng + dlng;
|
|
|
+ return new double[]{lng * 2 - mglng, lat * 2 - mglat};
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ public double transformlng(double lng, double lat) {
|
|
|
+ double ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * Math.sqrt(Math.abs(lng));
|
|
|
+ ret += ((20.0 * Math.sin(6.0 * lng * Math.PI) + 20.0 * Math.sin(2.0 * lng * Math.PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((20.0 * Math.sin(lng * Math.PI) + 40.0 * Math.sin((lng / 3.0) * Math.PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((150.0 * Math.sin((lng / 12.0) * Math.PI) + 300.0 * Math.sin((lng / 30.0) * Math.PI)) * 2.0) / 3.0;
|
|
|
+ return ret;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double transformlat(double lng, double lat) {
|
|
|
+ double ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * Math.sqrt(Math.abs(lng));
|
|
|
+ ret += ((20.0 * Math.sin(6.0 * lng * Math.PI) + 20.0 * Math.sin(2.0 * lng * Math.PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((20.0 * Math.sin(lat * Math.PI) + 40.0 * Math.sin((lat / 3.0) * Math.PI)) * 2.0) / 3.0;
|
|
|
+ ret += ((160.0 * Math.sin((lat / 12.0) * Math.PI) + 320 * Math.sin((lat * Math.PI) / 30.0)) * 2.0) / 3.0;
|
|
|
+ return ret;
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] wgs84ToCGCS2000(double lng, double lat) {
|
|
|
+
|
|
|
+ // 定义WGS84坐标点
|
|
|
+ ProjCoordinate wgs84Point = new ProjCoordinate(lng, lat);
|
|
|
+
|
|
|
+ // 将WGS84坐标点转换为CGCS2000坐标点
|
|
|
+ ProjCoordinate cgcs2000Point = new ProjCoordinate();
|
|
|
+
|
|
|
+ // 创建坐标转换器
|
|
|
+ CoordinateTransformFactory ctf = new CoordinateTransformFactory();
|
|
|
+ CoordinateTransform transform = ctf.createTransform(wgs84, cgcs2000);
|
|
|
+
|
|
|
+ transform.transform(wgs84Point, cgcs2000Point);
|
|
|
+ System.out.println(cgcs2000Point.x + "," + cgcs2000Point.y);
|
|
|
+ return new double[]{cgcs2000Point.x, cgcs2000Point.y};
|
|
|
+ }
|
|
|
+
|
|
|
+ public double[] cgcs2000ToWGS84(double lng, double lat) {
|
|
|
+ // 定义CGCS2000坐标点
|
|
|
+ ProjCoordinate cgcs2000Point = new ProjCoordinate(lng, lat);
|
|
|
+
|
|
|
+ // 将CGCS2000坐标点转换为WGS84坐标点
|
|
|
+ ProjCoordinate wgs84Point = new ProjCoordinate();
|
|
|
+
|
|
|
+ // 创建坐标转换器
|
|
|
+ CoordinateTransformFactory ctf = new CoordinateTransformFactory();
|
|
|
+ CoordinateTransform transform = ctf.createTransform(cgcs2000, wgs84);
|
|
|
+
|
|
|
+ transform.transform(cgcs2000Point, wgs84Point);
|
|
|
+ System.out.println(wgs84Point.x + "," + wgs84Point.y);
|
|
|
+ return new double[]{wgs84Point.x, wgs84Point.y};
|
|
|
+ }
|
|
|
+
|
|
|
+
|
|
|
+ public static double[] convertCGCS2000ToWGS84(double X, double Y) {
|
|
|
+ double x,y,L0 = 0,B,L;
|
|
|
+ x = X;
|
|
|
+ y = Y;
|
|
|
+ double p=206264.80625;
|
|
|
+ for (int i = 1; Y/i >=10; i = i * 10)//对Y坐标处理并求出中央子午线经度
|
|
|
+ {
|
|
|
+ y = Y - (int)(Y / i) * i-500000;
|
|
|
+ L0 =120;//中央经线，请完善代码去计算，这里处理浙江东阳的数据，偷懒直接指定了
|
|
|
+ }
|
|
|
+ //按6°带克氏椭球反算
|
|
|
+// double bt = x / 6367558.4969*p;
|
|
|
+// double BT = x / 6367558.4969;
|
|
|
+// double c3=Math.cos(BT)*Math.cos(BT);
|
|
|
+// double c4=Math.sin(BT)*Math.cos(BT);
|
|
|
+// double Bf=(bt+(50221746+(293622+(2350+22*c3)*c3)*c3)*c4*Math.pow(10,-10)*p)/p;
|
|
|
+// double c5=Math.pow(Math.cos(Bf),2);
|
|
|
+// double c6=Math.sin(Bf)*Math.cos(Bf);
|
|
|
+// double Nf=6399698.902-(21562.267-(108.973-0.612*c5)*c5)*c5;
|
|
|
+// double Z=y/(Nf*Math.cos(Bf));
|
|
|
+// double b2 = (0.5 + 0.003369 * c5) * c6;
|
|
|
+// double b3 = 0.333333 - (0.166667 - 0.001123 * c5) * c5;
|
|
|
+// double b4 = 0.25 + (0.16161 + 0.00562 * c5) * c5;
|
|
|
+// double b5=0.2-(0.1667-0.0088*c5)*c5;
|
|
|
+// double z2=Math.pow(Z,2);
|
|
|
+// B = (Bf*p - (1 - (b4 - 0.12 *z2) * z2) * z2 * b2 * p)/3600.0;
|
|
|
+// L = L0+((1 - (b3 - b5 * z2) * z2) * Z * p)/3600.0;
|
|
|
+// System.out.println("纬度："+B+" 经度:"+L+" 中央子午线："+L0);
|
|
|
+
|
|
|
+ // 按 3°带克氏椭球反算
|
|
|
+ double bt = x / 6367558.4969 * p;
|
|
|
+ double BT = x / 6367558.4969;
|
|
|
+ double c3 = Math.cos(BT) * Math.cos(BT);
|
|
|
+ double c4 = Math.sin(BT) * Math.cos(BT);
|
|
|
+ double Bf = (bt + (50221746 + (293622 + (2350 + 22 * c3) * c3) * c3) * c4 * Math.pow(10, -10) * p) / p;
|
|
|
+ double c5 = Math.pow(Math.cos(Bf), 2);
|
|
|
+ double c6 = Math.sin(Bf) * Math.cos(Bf);
|
|
|
+ double Nf = 6399698.902 - (21562.267 - (108.973 - 0.612 * c5) * c5) * c5;
|
|
|
+ double Z = y / (Nf * Math.cos(Bf));
|
|
|
+ double b2 = (0.5 + 0.003369 * c5) * c6;
|
|
|
+ double b3 = 0.333333 - (0.166667 - 0.001123 * c5) * c5;
|
|
|
+ double b4 = 0.25 + (0.16161 + 0.00562 * c5) * c5;
|
|
|
+ double b5 = 0.2 - (0.1667 - 0.0088 * c5) * c5;
|
|
|
+ double z2 = Math.pow(Z, 2);
|
|
|
+ // 3 度带中央子午线计算公式与 6 度带不同
|
|
|
+ double L0ThreeDegree = ((int)(L0 / 3) + 1) * 3 - 3;
|
|
|
+ B = (Bf * p - (1 - (b4 - 0.12 * z2) * z2) * z2 * b2 * p) / 3600.0;
|
|
|
+ L = L0ThreeDegree + ((1 - (b3 - b5 * z2) * z2) * Z * p) / 3600.0;
|
|
|
+// System.out.println("纬度：" + B + " 经度:" + L + " 中央子午线：" + L0ThreeDegree);
|
|
|
+ return new double[]{L, B};
|
|
|
+ }
|
|
|
+
|
|
|
+
|
|
|
+}
|
zhanghan