|
@@ -12,6 +12,7 @@
|
|
|
|
|
|
void lonlat_to_gauss(const double central, const double lon, const double lat, double *gx, double *gy)
|
|
|
{
|
|
|
+ // 将纬度和经度转换为弧度
|
|
|
double lamb = DEG2RAD_C(lat);
|
|
|
double phi = DEG2RAD_C(lon);
|
|
|
|
|
@@ -24,6 +25,7 @@ void lonlat_to_gauss(const double central, const double lon, const double lat, d
|
|
|
double C = WGS84_E2_C * cos(lamb) * cos(lamb) / (1 - WGS84_E2_C);
|
|
|
double A = (phi - centralMeridianRad) * cos(lamb);
|
|
|
|
|
|
+ // 计算子午线弧长
|
|
|
double M = WGS84_A_C * ((1 - WGS84_E2_C / 4.0 - 3.0 * WGS84_E2_C * WGS84_E2_C / 64.0 - 5.0 * WGS84_E2_C * WGS84_E2_C * WGS84_E2_C / 256.0) * lamb -
|
|
|
(3.0 * WGS84_E2_C / 8.0 + 3.0 * WGS84_E2_C * WGS84_E2_C / 32.0 + 45.0 * WGS84_E2_C * WGS84_E2_C * WGS84_E2_C / 1024.0) * sin(2.0 * lamb) +
|
|
|
(15.0 * WGS84_E2_C * WGS84_E2_C / 256.0 + 45.0 * WGS84_E2_C * WGS84_E2_C * WGS84_E2_C / 1024.0) * sin(4.0 * lamb) - (35.0 * WGS84_E2_C * WGS84_E2_C * WGS84_E2_C / 3072.0) * sin(6.0 * lamb));
|
|
@@ -35,19 +37,30 @@ void lonlat_to_gauss(const double central, const double lon, const double lat, d
|
|
|
|
|
|
void gauss_to_lonlat(const double central, const double gx, const double gy, double *lon, double *lat)
|
|
|
{
|
|
|
+ // 将高斯平面坐标转换为经纬度坐标
|
|
|
double Y = gx;
|
|
|
double X = gy;
|
|
|
|
|
|
+ // 去除中央经线的偏移量
|
|
|
Y -= 500000;
|
|
|
+ // 计算第一偏心率的平方根倒数
|
|
|
double e1 = (1 - sqrt(1 - WGS84_E2_C)) / (1 + sqrt(1 - WGS84_E2_C));
|
|
|
+ // 计算子午线弧长
|
|
|
double M = X;
|
|
|
+ // 计算子午线曲率半径
|
|
|
double mu = M / (WGS84_A_C * (1 - WGS84_E2_C / 4.0 - 3 * WGS84_E2_C * WGS84_E2_C / 64.0 - 5 * WGS84_E2_C * WGS84_E2_C * WGS84_E2_C / 256.0));
|
|
|
+ // 计算纬度
|
|
|
double phi1 = mu + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * sin(2 * mu) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * sin(4 * mu) + (151 * e1 * e1 * e1 / 96) * sin(6 * mu) + (1097 * e1 * e1 * e1 * e1 / 512) * sin(8 * mu);
|
|
|
|
|
|
+ // 计算子午线曲率半径
|
|
|
double C1 = WGS84_E2_C * cos(phi1) * cos(phi1) / (1 - WGS84_E2_C);
|
|
|
+ // 计算纬度曲率半径
|
|
|
double T1 = tan(phi1) * tan(phi1);
|
|
|
+ // 计算卯酉圈曲率半径
|
|
|
double N1 = WGS84_A_C / sqrt(1 - WGS84_E2_C * sin(phi1) * sin(phi1));
|
|
|
+ // 计算子午线弧长
|
|
|
double R1 = WGS84_A_C * (1 - WGS84_E2_C) / pow(1 - WGS84_E2_C * sin(phi1) * sin(phi1), 1.5);
|
|
|
+ // 计算子午线弧长
|
|
|
double D = Y / N1;
|
|
|
|
|
|
// 经纬度计算
|
|
@@ -84,24 +97,33 @@ void mercator_to_lonlat(const double mctx, const double mcty, double *lon, doubl
|
|
|
*lat = RAD2DEG_C(2 * atan(exp(mcty / WGS84_A_C)) - M_PI / 2);
|
|
|
}
|
|
|
|
|
|
+// 将高斯投影坐标转换为墨卡托投影坐标
|
|
|
void gauss_to_mercator(const double central, const double gx, const double gy, double *mctx, double *mcty)
|
|
|
{
|
|
|
+ // 定义经度和纬度变量
|
|
|
double lon;
|
|
|
double lat;
|
|
|
+ // 将高斯投影坐标转换为经纬度坐标
|
|
|
gauss_to_lonlat(central, gx, gy, &lon, &lat);
|
|
|
+ // 将经纬度坐标转换为墨卡托投影坐标
|
|
|
lonlat_to_mercator(lon, lat, mctx, mcty);
|
|
|
}
|
|
|
|
|
|
+// 将墨卡托坐标转换为高斯坐标
|
|
|
void mercator_to_gauss(const double central, const double mctx, const double mcty, double *gx, double *gy)
|
|
|
{
|
|
|
+ // 定义经纬度变量
|
|
|
double lon;
|
|
|
double lat;
|
|
|
+ // 将墨卡托坐标转换为经纬度
|
|
|
mercator_to_lonlat(mctx, mcty, &lon, &lat);
|
|
|
+ // 将经纬度转换为高斯坐标
|
|
|
lonlat_to_gauss(central, lon, lat, gx, gy);
|
|
|
}
|
|
|
|
|
|
void lonlat_to_ecef(const double lon, const double lat, const double height, double* x, double* y, double* z)
|
|
|
{
|
|
|
+ // 将经纬度转换为弧度
|
|
|
double lamb = DEG2RAD_C(lat);
|
|
|
double phi = DEG2RAD_C(lon);
|
|
|
|