添加坐标系转换算法

coordinate_convert/c++/studio_geo_const.h

@@ -0,0 +1,60 @@
+/*****************************************************************//**
+ * \file   studio_geo_const.h
+ * \brief  
+ * 
+ * \author wangyingjie
+ * \date   April 2025
+ *********************************************************************/
+
+#ifndef STUDIO_GEO_CONST_H
+#define STUDIO_GEO_CONST_H
+
+namespace EarthRadius
+{
+static constexpr double MEAN = 6371000.0;        // 平均半径 米
+static constexpr double EQUATORIAL = 6378137.0;  // 赤道半径 米
+static constexpr double POLAR = 6356752.3142;    // 极地半径 米
+}  // namespace EarthRadius
+namespace WGS84
+{
+static constexpr double A = 6378137.0;         // 长半轴 米
+static constexpr double INVF = 298.257223563;  // 扁率的倒数
+static constexpr double F = (1.0 / INVF);      // 扁率
+static constexpr double B = (A * (1.0 - F));   // 短半轴
+static constexpr double E2 = (2 * F - F * F);  // 第一偏心扁率的平方
+static constexpr double GM = 3.986004418e14;   // 地心引力常数 立方米/秒的平方
+static constexpr double W = 7.2921151467e-5;   // 自转角速度 弧度每秒
+static constexpr double J2 = 1.08262983226e-3;
+}  // namespace WGS84
+namespace CGCS2000
+{
+// CGCS2000 大多数参数跟WGS84 一致
+static constexpr double A = 6378137.0;
+static constexpr double INVF = 298.257222101;
+static constexpr double F = (1.0 / INVF);
+static constexpr double B = (A * (1.0 - F));
+static constexpr double E2 = (2 * F - F * F);
+static constexpr double GM = 3.9860044181e14;   // 地心引力常数 立方米/秒的平方
+static constexpr double W = 7.2921151467e-5;    // 自转角速度 弧度每秒
+static constexpr double J2 = 1.08262983226e-3;  // 重力场谐系数
+}  // namespace CGCS2000
+namespace XIAN80
+{
+
+static constexpr double A = 6378140.0;     // 长半轴
+static constexpr double INVF = 298.257;    // 扁率的倒数
+static constexpr double F = (1.0 / INVF);  // 扁率
+static constexpr double B = (A * (1.0 - F));
+static constexpr double E2 = (2 * F - F * F);
+}  // namespace XIAN80
+namespace BEIJING54
+{
+
+static constexpr double A = 6378245.0;     // 长半轴
+static constexpr double INVF = 298.3;      // 扁率的倒数
+static constexpr double F = (1.0 / INVF);  // 扁率
+static constexpr double B = (A * (1.0 - F));
+static constexpr double E2 = (2 * F - F * F);
+}  // namespace BEIJING54
+
+#endif  // STUDIO_GEO_CONST_H

coordinate_convert/c++/studio_proj.cpp

@@ -0,0 +1,208 @@
+
+#pragma once
+#include "studio_proj.h"
+
+void proj::gauss_to_lonlat(const double& gx, const double& gy, double& lon, double& lat, const proj::param& p)
+{
+    double centralMeridian = p.central;
+
+    // 辅助项
+    double n = (p.major_axis - p.major_axis * (1 - 1 / p.flatten)) / (p.major_axis + p.major_axis * (1 - 1 / p.flatten));
+    double semiMinorAxis = p.major_axis * (1 - 1 / p.flatten);
+    double firstEccentricitySquared = (p.major_axis * p.major_axis - semiMinorAxis * semiMinorAxis) / (p.major_axis * p.major_axis);
+    double secondEccentricitySquared = (p.major_axis * p.major_axis - semiMinorAxis * semiMinorAxis) / (semiMinorAxis * semiMinorAxis);
+
+    // 计算Y坐标对应的纬度
+    double mu = (gy - 10000000.0) / (p.scale * semiMinorAxis);
+    double phi = mu + (3 * n / 2 - 27 * n * n * n / 32) * sin(2 * mu) + (21 * n * n / 16 - 55 * n * n * n / 32) * sin(4 * mu) + (151 * n * n * n / 96) * sin(6 * mu) + (1097 * n * n * n * n / 512) * sin(8 * mu);
+
+    // 计算X坐标对应的经度
+    double lambda = gx - p.easting;
+    double lambdaPrime = lambda / (p.major_axis * cos(phi));
+
+    // 迭代计算纬度
+    double phiPrime = phi;
+    double phiNew = phi + (phiPrime - phi - lambdaPrime * cos(phi)) / (1 - firstEccentricitySquared * sin(phi) * sin(phi) - lambdaPrime * lambdaPrime * cos(phi) * cos(phi));
+
+    while (fabs(phiNew - phiPrime) > 1e-10)
+    {
+        phiPrime = phiNew;
+        phiNew = phi + (phiPrime - phi - lambdaPrime * cos(phiPrime)) / (1 - firstEccentricitySquared * sin(phiPrime) * sin(phiPrime) - lambdaPrime * lambdaPrime * cos(phiPrime) * cos(phiPrime));
+    }
+
+    phi = phiNew;
+
+    // 转换为经纬度
+    lat = RAD2DEG(phi);
+    lon = RAD2DEG((centralMeridian + lambdaPrime * 180.0 / (AO_PI * p.major_axis * cos(phi))));
+}
+
+void proj::gauss_to_lonlat(const double& central, const double& gx, const double& gy, double& lon, double& lat)
+{
+    double Y = gx, X = gy;
+
+    Y -= 500000;
+    double e1 = (1 - sqrt(1 - WGS84::E2)) / (1 + sqrt(1 - WGS84::E2));
+    double M = X;
+    double mu = M / (WGS84::A * (1 - WGS84::E2 / 4.0 - 3 * WGS84::E2 * WGS84::E2 / 64.0 - 5 * WGS84::E2 * WGS84::E2 * WGS84::E2 / 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 * cos(phi1) * cos(phi1) / (1 - WGS84::E2);
+    double T1 = tan(phi1) * tan(phi1);
+    double N1 = WGS84::A / sqrt(1 - WGS84::E2 * sin(phi1) * sin(phi1));
+    double R1 = WGS84::A * (1 - WGS84::E2) / pow(1 - WGS84::E2 * sin(phi1) * sin(phi1), 1.5);
+    double D = Y / N1;
+
+    // 经纬度计算
+    double phi = phi1 - (N1 * tan(phi1) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * e1) * D * D * D * D / 24 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * e1 - 3 * C1 * C1) * D * D * D * D * D * D / 720);
+
+    double lambda = DEG2RAD(central) + (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * e1 + 24 * T1 * T1) * D * D * D * D * D / 120) / cos(phi1);
+
+    // 转换为度
+    lat = RAD2DEG(phi);
+    lon = RAD2DEG(lambda);
+}
+
+void proj::lonlat_to_gauss(const double& lon, const double& lat, double& gx, double& gy, const proj::param& p)
+{
+    // 中央经线
+    double L0 = p.central;
+
+    // 将经纬度转换为弧度
+    double lon_rad = DEG2RAD(lon);
+    double lat_rad = DEG2RAD(lat);
+
+    // 计算经差
+    double lambda = (lon_rad - DEG2RAD(L0));
+
+    // 辅助项
+    double n = (p.major_axis - p.major_axis * (1 - 1 / p.flatten)) / (p.major_axis + p.major_axis * (1 - 1 / p.flatten));
+    double semiMinorAxis = p.major_axis * (1 - 1 / p.flatten);
+    double firstEccentricitySquared = (p.major_axis * p.major_axis - semiMinorAxis * semiMinorAxis) / (p.major_axis * p.major_axis);
+    double secondEccentricitySquared = (p.major_axis * p.major_axis - semiMinorAxis * semiMinorAxis) / (semiMinorAxis * semiMinorAxis);
+
+    // 计算N和T
+    double radiusOfCurvature = p.major_axis / sqrt(1 - firstEccentricitySquared * sin(lat_rad) * sin(lat_rad));
+    double tangentSquared = tan(lat_rad) * tan(lat_rad);
+    double meridianCurvature = secondEccentricitySquared * cos(lat_rad) * cos(lat_rad);
+
+    // 计算A和B
+    double A = lambda * cos(lat_rad);
+    double B = lambda * lambda * cos(lat_rad) * cos(lat_rad) / 2;
+
+    // X坐标
+    gx = radiusOfCurvature * (lambda + (1 - tangentSquared + meridianCurvature) * pow(A, 3) / 6 +
+                              (5 - 18 * tangentSquared + tangentSquared * tangentSquared + 72 * meridianCurvature - 58 * (1 - firstEccentricitySquared) / (1 + firstEccentricitySquared)) * pow(A, 5) / 120) +
+         p.easting;
+
+    // Y坐标
+    gy = p.scale * (semiMinorAxis * (lat_rad - (1 - n + (5 * n * n) / 4 - (5 * n * n * n) / 4) * sin(2 * lat_rad) / 2 + (1 - (3 * n) / 2 + (21 * n * n) / 16 - (55 * n * n * n) / 32) * sin(4 * lat_rad) / 24 -
+                                     (1 - (11 * n) / 4 + (413 * n * n) / 96 - (1231 * n * n * n) / 256) * sin(6 * lat_rad) / 720 + (1 - (15 * n) / 4 + (517 * n * n) / 64 - (5147 * n * n * n) / 256) * sin(8 * lat_rad) / 40320)) +
+         10000000.0;
+}
+
+void proj::lonlat_to_gauss(const double& central, const double& lon, const double& lat, double& gx, double& gy)
+{
+    double lamb = DEG2RAD(lat);
+    double phi = DEG2RAD(lon);
+
+    // 将中央子午线转换为弧度
+    double centralMeridianRad = DEG2RAD(central);
+
+    // 计算高斯-克吕格投影公式中的参数
+    double N = WGS84::A / sqrt(1 - WGS84::E2 * sin(lamb) * sin(lamb));
+    double T = tan(lamb) * tan(lamb);
+    double C = WGS84::E2 * cos(lamb) * cos(lamb) / (1 - WGS84::E2);
+    double A = (phi - centralMeridianRad) * cos(lamb);
+
+    double M = WGS84::A * ((1 - WGS84::E2 / 4.0 - 3.0 * WGS84::E2 * WGS84::E2 / 64.0 - 5.0 * WGS84::E2 * WGS84::E2 * WGS84::E2 / 256.0) * lamb -
+                           (3.0 * WGS84::E2 / 8.0 + 3.0 * WGS84::E2 * WGS84::E2 / 32.0 + 45.0 * WGS84::E2 * WGS84::E2 * WGS84::E2 / 1024.0) * sin(2.0 * lamb) +
+                           (15.0 * WGS84::E2 * WGS84::E2 / 256.0 + 45.0 * WGS84::E2 * WGS84::E2 * WGS84::E2 / 1024.0) * sin(4.0 * lamb) - (35.0 * WGS84::E2 * WGS84::E2 * WGS84::E2 / 3072.0) * sin(6.0 * lamb));
+
+    // 计算 X, Y 坐标
+    gy = M + N * tan(lamb) * (A * A / 2.0 + (5.0 - T + 9.0 * C + 4.0 * C * C) * A * A * A * A / 24.0 + (61.0 - 58.0 * T + T * T + 600.0 * C - 330.0 * WGS84::E2) * A * A * A * A * A * A / 720.0);
+    gx = N * (A + (1.0 - T + C) * A * A * A / 6.0 + (5.0 - 18.0 * T + T * T + 72.0 * C - 58.0 * WGS84::E2) * A * A * A * A * A / 120.0) + 500000.0;  // 中央子午线偏移+500000.0
+}
+
+void proj::mercator_to_lonlat(const double& mctx, const double& mcty, double& lon, double& lat)
+{
+    lon = RAD2DEG(mctx / WGS84::A);
+    lat = RAD2DEG(2 * atan(exp(mcty / WGS84::A)) - AO_PI / 2);
+}
+
+void proj::lonlat_to_mercator(const double& lon, const double& lat, double& mctx, double& mcty)
+{
+    double lamb = DEG2RAD(lat);
+    double phi = DEG2RAD(lon);
+
+    // 墨卡托投影公式
+    mctx = WGS84::A * phi;
+    mcty = WGS84::A * log(tan(AO_PI / 4 + lamb / 2));
+}
+
+void proj::mercator_to_gauss(const double& mctx, const double& mcty, double& gx, double& gy, const proj::param& p)
+{
+    double lon, lat;
+    mercator_to_lonlat(mctx, mcty, lon, lat);
+    lonlat_to_gauss(lon, lat, gx, gy, p);
+}
+
+void proj::gauss_to_mercator(const double& gx, const double& gy, double& mctx, double& mcty, const proj::param& p)
+{
+    double lon, lat;
+    gauss_to_lonlat(gx, gy, lon, lat, p);
+    lonlat_to_mercator(lon, lat, mctx, mcty);
+}
+
+void proj::mercator_to_gauss(const double& central, const double& mctx, const double& mcty, double& gx, double& gy)
+{
+    double lon, lat;
+    mercator_to_lonlat(mctx, mcty, lon, lat);
+    lonlat_to_gauss(central, lon, lat, gx, gy);
+}
+
+void proj::gauss_to_mercator(const double& central, const double& gx, const double& gy, double& mctx, double& mcty)
+{
+    double lon, lat;
+    gauss_to_lonlat(central, gx, gy, lon, lat);
+    lonlat_to_mercator(lon, lat, mctx, mcty);
+}
+
+void proj::ecef_to_lonlat(const double& x, const double& y, const double& z, double& lon, double& lat, double& height)
+{
+    // 计算经度
+    lon = std::atan2(y, x);
+
+    // 计算初始纬度估计
+    double p = std::sqrt(x * x + y * y);
+    double theta = std::atan2(z, p * (1 - WGS84::F));
+    lat = std::atan2(z + WGS84::E2 * WGS84::B * std::pow(std::sin(theta), 3), p - WGS84::E2 * WGS84::B * std::pow(std::cos(theta), 3));
+
+    // 迭代计算纬度，直到收敛
+    double previousLatitude;
+    do
+    {
+        previousLatitude = lat;
+        double N = WGS84::A / std::sqrt(1 - WGS84::E2 * std::sin(lat) * std::sin(lat));
+        height = p / std::cos(lat) - N;
+        lat = std::atan2(z + WGS84::E2 * N * sin(lat), p);
+    } while (std::fabs(lat - previousLatitude) > 1e-12);  // 收敛条件
+
+    // 将纬度和经度转换为度
+    lat = RAD2DEG(lat);
+    lon = RAD2DEG(lon);
+}
+
+void proj::lonlat_to_ecef(const double& lon, const double& lat, const double& height, double& x, double& y, double& z)
+{
+    double lamb = DEG2RAD(lat);
+    double phi = DEG2RAD(lon);
+
+    // 计算N（曲率半径）
+    double N = WGS84::A / std::sqrt(1 - WGS84::E2 * std::sin(lamb) * std::sin(lamb));
+
+    // 计算XYZ坐标
+    x = (N + height) * std::cos(lamb) * std::cos(phi);
+    y = (N + height) * std::cos(lamb) * std::sin(phi);
+    z = (WGS84::B * WGS84::B / (WGS84::A * WGS84::A) * N + height) * std::sin(lamb);
+}

coordinate_convert/c++/studio_proj.h

@@ -0,0 +1,160 @@
+/**
+******************************************************************************
+* @file           : studio_proj.h
+* @author         : wangyingjie
+* @brief          : 经纬度 转 高斯投影，墨卡托投影，地心地固直角坐标系 C 语言版
+ *                : 中央经线计算公式： double central = static_cast<int>(lon / 3) * 3;
+* @attention      : None
+* @date           : 2025/5/12
+******************************************************************************
+*/
+
+#ifndef STUDIO_PROJ_H
+#define STUDIO_PROJ_H
+
+#include "studio_macros.h"
+#include "studio_geo_const.h"
+
+#define AO_GAUSS3_NO(lon) static_cast<int>(lon / 3.0 + 1)
+#define AO_GAUSS6_NO(lon) static_cast<int>(lon / 6.0 + 1)
+#define AO_GAUSS3_L0(no) static_cast<double>(3 * no - 1.5)
+#define AO_GAUSS6_L0(no) static_cast<double>(6 * no - 3.0)
+
+class proj
+{
+  public:
+    struct param
+    {
+        // 中央经度线
+        double central{117.0};
+        double scale{0.9996};
+        double easting{500000.0};  // 东偏移
+        // 长半轴
+        double major_axis{WGS84::A};
+        // 扁率
+        double flatten{WGS84::F};
+    };
+
+  public:
+    /// <summary>
+    /// 高斯投影反算经纬度
+    /// </summary>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    /// <param name="p"></param>
+    static void gauss_to_lonlat(const double& gx, const double& gy, double& lon, double& lat, const param& p);
+
+    /// <summary>
+    /// 高斯投影反算经纬度,使用静态值速度会更快
+    /// </summary>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    static void gauss_to_lonlat(const double& central, const double& gx, const double& gy, double& lon, double& lat);
+
+    /// <summary>
+    /// 经纬度转高斯投影
+    /// </summary>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    /// <param name="p"></param>
+    static void lonlat_to_gauss(const double& lon, const double& lat, double& gx, double& gy, const param& p);
+
+    /// <summary>
+    /// 经纬度转高斯投影,使用静态值速度会更快
+    /// </summary>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    static void lonlat_to_gauss(const double& central, const double& lon, const double& lat, double& gx, double& gy);
+
+    /// <summary>
+    /// 墨卡托转经纬度
+    /// </summary>
+    /// <param name="mctx"></param>
+    /// <param name="mcty"></param>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    static void mercator_to_lonlat(const double& mctx, const double& mcty, double& lon, double& lat);
+
+    /// <summary>
+    /// 经纬度转墨卡托
+    /// </summary>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    /// <param name="mctx"></param>
+    /// <param name="mcty"></param>
+    static void lonlat_to_mercator(const double& lon, const double& lat, double& mctx, double& mcty);
+
+    /// <summary>
+    /// 墨卡托转高斯投影
+    /// </summary>
+    /// <param name="mctx"></param>
+    /// <param name="mcty"></param>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    /// <param name="p"></param>
+    static void mercator_to_gauss(const double& mctx, const double& mcty, double& gx, double& gy, const param& p);
+
+    /// <summary>
+    /// 高斯投影转墨卡托
+    /// </summary>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    /// <param name="mctx"></param>
+    /// <param name="mcty"></param>
+    /// <param name="p"></param>
+    static void gauss_to_mercator(const double& gx, const double& gy, double& mctx, double& mcty, const param& p);
+
+    /// <summary>
+    /// 墨卡托转高斯投影,使用静态值速度会更快
+    /// </summary>
+    /// <param name="central"></param>
+    /// <param name="mctx"></param>
+    /// <param name="mcty"></param>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    static void mercator_to_gauss(const double& central, const double& mctx, const double& mcty, double& gx, double& gy);
+
+    /// <summary>
+    /// 高斯投影转墨卡托,使用静态值速度会更快
+    /// </summary>
+    /// <param name="central"></param>
+    /// <param name="gx"></param>
+    /// <param name="gy"></param>
+    /// <param name="mctx"></param>
+    /// <param name="mcty"></param>
+    static void gauss_to_mercator(const double& central, const double& gx, const double& gy, double& mctx, double& mcty);
+
+    /// <summary>
+    /// 地心地固直角坐标系 转 经纬坐标系
+    /// </summary>
+    /// <param name="x"></param>
+    /// <param name="y"></param>
+    /// <param name="z"></param>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    /// <param name="height"></param>
+    static void ecef_to_lonlat(const double& x, const double& y, const double& z, double& lon, double& lat, double& height);
+
+    /// <summary>
+    /// 经纬坐标系 转 地心地固直角坐标系
+    /// </summary>
+    /// <param name="lon"></param>
+    /// <param name="lat"></param>
+    /// <param name="height"></param>
+    /// <param name="x"></param>
+    /// <param name="y"></param>
+    /// <param name="z"></param>
+    static void lonlat_to_ecef(const double& lon, const double& lat, const double& height, double& x, double& y, double& z);
+
+  private:
+};
+
+#endif  //

coordinate_convert/c/studio_proj_c.c

@@ -0,0 +1,140 @@
+/**
+******************************************************************************
+* @file           : studio_proj_c.c
+* @author         : wyj
+* @brief          : 高斯投影 C 语言版
+* @attention      : None
+* @date           : 2025/5/9
+******************************************************************************
+*/
+
+#include "studio_proj_c.h"
+
+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);
+
+    // 将中央子午线转换为弧度
+    double centralMeridianRad = DEG2RAD_C(central);
+
+    // 计算高斯-克吕格投影公式中的参数
+    double N = WGS84_A_C / sqrt(1 - WGS84_E2_C * sin(lamb) * sin(lamb));
+    double T = tan(lamb) * tan(lamb);
+    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));
+
+    // 计算 X, Y 坐标
+    *gy = M + N * tan(lamb) * (A * A / 2.0 + (5.0 - T + 9.0 * C + 4.0 * C * C) * A * A * A * A / 24.0 + (61.0 - 58.0 * T + T * T + 600.0 * C - 330.0 * WGS84_E2_C) * A * A * A * A * A * A / 720.0);
+    *gx = N * (A + (1.0 - T + C) * A * A * A / 6.0 + (5.0 - 18.0 * T + T * T + 72.0 * C - 58.0 * WGS84_E2_C) * A * A * A * A * A / 120.0) + 500000.0;  // 中央子午线偏移+500000.0
+}
+
+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;
+
+    // 经纬度计算
+    double phi = phi1 - (N1 * tan(phi1) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * e1) * D * D * D * D / 24 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * e1 - 3 * C1 * C1) * D * D * D * D * D * D / 720);
+
+    double lambda = DEG2RAD_C(central) + (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * e1 + 24 * T1 * T1) * D * D * D * D * D / 120) / cos(phi1);
+
+    // 转换为度
+    *lat = RAD2DEG_C(phi);
+    *lon = RAD2DEG_C(lambda);
+}
+
+void lonlat_to_mercator(const double lon, const double lat, double *mctx, double *mcty)
+{
+    // 经度转换：phi = lon * (M_PI / 180.0)
+    double phi = DEG2RAD_C(lon);
+
+    // 纬度转换：lamb = lat * (M_PI / 180.0)
+    double lamb = DEG2RAD_C(lat);
+
+    // 墨卡托投影公式：mctx = WGS84_A_C * phi
+    *mctx = WGS84_A_C * phi;
+
+    // 墨卡托投影公式：mcty = WGS84_A_C * log(tan(M_PI / 4 + lamb / 2))
+    *mcty = WGS84_A_C * log(tan(M_PI / 4 + lamb / 2));
+}
+
+void mercator_to_lonlat(const double mctx, const double mcty, double *lon, double *lat)
+{
+    // 经度转换：mctx / WGS84_A_C
+    *lon = RAD2DEG_C(mctx / WGS84_A_C);
+
+    // 纬度转换：使用逆墨卡托公式
+    *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);
+
+    // 计算N (曲率半径)
+    double N = WGS84_A_C / sqrt(1 - WGS84_E2_C * sin(lamb) * sin(lamb));
+
+    // 计算XYZ坐标
+    *x = (N + height) * cos(lamb) * cos(phi);
+    *y = (N + height) * cos(lamb) * sin(phi);
+    *z = (WGS84_B_C * WGS84_B_C / (WGS84_A_C * WGS84_A_C) * N + height) * sin(lamb);
+}
+
+void ecef_to_lonlat(const double x, const double y, const double z, double* lon, double* lat, double* height)
+{
+    // 计算经度
+    *lon = atan2(y, x);
+
+    // 计算初始纬度估计
+    double p = sqrt(x * x + y * y);
+    double theta = atan2(z, p * (1 - WGS84_F_C));
+    *lat = atan2(z + WGS84_E2_C * WGS84_B_C * pow(sin(theta), 3), p - WGS84_E2_C * WGS84_B_C * pow(cos(theta), 3));
+
+    // 迭代计算纬度，直到收敛
+    double previousLatitude;
+    do
+    {
+        previousLatitude = *lat;
+        double N = WGS84_A_C / sqrt(1 - WGS84_E2_C * sin(*lat) * sin(*lat));
+        *height = p / cos(*lat) - N;
+        *lat = atan2(z + WGS84_E2_C * N * sin(*lat), p);
+    } while (fabs(*lat - previousLatitude) > 1e-12);  // 收敛条件
+
+    // 将纬度和经度转换为度
+    *lat = RAD2DEG_C(*lat);
+    *lon = RAD2DEG_C(*lon);
+}

coordinate_convert/c/studio_proj_c.h

@@ -0,0 +1,91 @@
+/**
+******************************************************************************
+* @file           : studio_proj_c.h
+* @author         : wangyingjie
+* @brief          : 经纬度 转 高斯投影，墨卡托投影，地心地固直角坐标系 C 语言版
+ *                : 中央经线计算公式： double central = (int)(lon / 3) * 3;
+* @attention      : None
+* @date           : 2025/5/9
+******************************************************************************
+*/
+
+#ifndef STUDIO_PROJ_C_H
+#define STUDIO_PROJ_C_H
+
+#include <math.h>
+
+#define DEG2RAD_C(deg) ((deg) * M_PI / 180.0)
+#define RAD2DEG_C(rad) ((rad) * 180.0 / M_PI)
+
+#define WGS84_A_C 6378137.0                                     // 长半轴
+#define WGS84_E2_C 0.0066943799901413165                        // 第一偏心扁率的平方
+#define INVF_C 298.257223563                                    // 扁率的倒数
+#define WGS84_B_C 6356752.3142                                  // 短半轴
+#define WGS84_F_C 1 / INVF_C                                    // 扁率
+#define WGS84_E2_C ((2 * WGS84_F_C) - (WGS84_F_C * WGS84_F_C))  // 第一偏心扁率的平方
+
+/// 经纬度 转 高斯克吕格投影
+/// \param central 中央经线 (度)
+/// \param lon 经度 (度)
+/// \param lat 纬度 (度)
+/// \param gx 高斯坐标  x(m)
+/// \param gy 高斯坐标 y(m)
+void lonlat_to_gauss(const double central, const double lon, const double lat, double *gx, double *gy);
+
+/// 高斯克吕格投影 转 经纬度
+/// \param central 中央经线 (度)
+/// \param gx 高斯坐标  x(m)
+/// \param gy 高斯坐标 y(m)
+/// \param lon 经度 (度)
+/// \param lat 纬度 (度)
+void gauss_to_lonlat(const double central, const double gx, const double gy, double *lon, double *lat);
+
+/// 经纬度 转 墨卡托投影
+/// \param lon 经度 (度)
+/// \param lat 纬度 (度)
+/// \param mctx 墨卡托投影 x(m)
+/// \param mcty 墨卡托投影 y(m)
+void lonlat_to_mercator(const double lon, const double lat, double *mctx, double *mcty);
+
+/// 墨卡托投影 转 经纬度
+/// \param mctx 墨卡托投影 x(m)
+/// \param mcty 墨卡托投影 y(m)
+/// \param lon 经度 (度)
+/// \param lat 纬度 (度)
+void mercator_to_lonlat(const double mctx, const double mcty, double *lon, double *lat);
+
+/// 墨卡托投影 转 高斯投影
+/// \param central 中央经线 (度)
+/// \param gx 墨卡托投影 x(m)
+/// \param gy 墨卡托投影 y(m)
+/// \param mctx 高斯投影 x(m)
+/// \param mcty 高斯投影 y(m)
+void gauss_to_mercator(const double central, const double gx, const double gy, double *mctx, double *mcty);
+
+/// 高斯投影 转 墨卡托投影
+/// \param central 中央经线 (度)
+/// \param mctx 墨卡托投影 x(m)
+/// \param mcty 墨卡托投影 y(m)
+/// \param gx 高斯投影 x(m)
+/// \param gy 高斯投影 y(m)
+void mercator_to_gauss(const double central, const double mctx, const double mcty, double *gx, double *gy);
+
+/// 经纬坐标系 转 地心地固直角坐标系
+/// \param lon 经度 (度)
+/// \param lat 纬度 (度)
+/// \param height 高程 (m)
+/// \param x 地心地固直角坐标系 x (m)
+/// \param y 地心地固直角坐标系 y (m)
+/// \param z 地心地固直角坐标系 z (m)
+void lonlat_to_ecef(const double lon, const double lat, const double height, double *x, double *y, double *z);
+
+/// 经纬坐标系 转 地心地固直角坐标系
+/// \param x 地心地固直角坐标系 x (m)
+/// \param y 地心地固直角坐标系 y (m)
+/// \param z 地心地固直角坐标系 z (m)
+/// \param lon 经度 (度)
+/// \param lat 纬度 (度)
+/// \param height 高程 (m)
+void ecef_to_lonlat(const double x, const double y, const double z, double *lon, double *lat, double *height);
+
+#endif  // STUDIO_PROJ_C_H

coordinate_convert/c/task_alog_c.cpp

@@ -0,0 +1,88 @@
+/**
+******************************************************************************
+* @file           : task_alog_c.cpp
+* @author         : wyj
+* @brief          : C语言语法测试
+* @attention      : None
+* @date           : 2025/5/9
+******************************************************************************
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "geography/studio_proj_c.h"
+
+#ifdef __cplusplus
+}
+#endif
+
+#include <stdio.h>
+#include <stdio.h>
+#include "geometry/studio_geo_utils.h"
+#include "geometry/studio_geo_algo.h"
+
+int main() {
+    printf("\n\n===================== %s =====================\n\n", __FILE__);
+    silly::geo::utils::init_gdal_env();
+    std::string path;
+#ifdef IS_WINDOWS
+    path = "D:/5_file/2_readfile/geojson/multi_point/fitting_examples.geojson";
+#else
+    path = "/home/wyj/myself/2_data/2_geojson/multi_point/fitting_examples.geojson";
+#endif
+
+    std::cout << "path: " << path << std::endl;
+    std::vector<studio_geo_coll> res_collections;
+    std::vector<studio_geo_coll> collections;
+    silly::geo::utils::read_geo_coll(path, collections);
+    for (auto &coll: collections) {
+        // ------------- 转换为高斯投影 -------------
+        studio_line gauss_line;
+        double central = static_cast<int>(coll.m_line[0].x / 3) * 3;
+        for (auto &point: coll.m_line) {
+            double gx = 0.0;
+            double gy = 0.0;
+            lonlat_to_gauss(central, point.x, point.y, &gx, &gy);
+            gauss_line.push_back(studio_point(gx, gy));
+        }
+
+        // 简化线段，目标点数为28个
+        int max_points = 28;
+        double epsilon = 0.001;
+
+        // studio_line gs_simplified_line = simplify_line_2(gauss_line, max_points);
+        studio_line gs_simplified_line;
+        bool res = geo_vacuate::vacuate(gauss_line, max_points, epsilon, gs_simplified_line);
+        if (!res) {
+            std::cout << "Failed to simplify line." << std::endl;
+            return 1;
+        }
+        studio_line simplified_line;
+        // 高斯投影在转回经纬度
+        for (auto &point: gs_simplified_line) {
+            double lon = 0.0;
+            double lat = 0.0;
+            gauss_to_lonlat(central, point.x, point.y, &lon, &lat);
+            simplified_line.push_back(studio_point(lon, lat));
+        }
+
+        studio_geo_coll temp;
+        temp.m_type = enum_geometry_type::egtLineString;
+        temp.m_line = simplified_line;
+        res_collections.push_back(temp);
+        break;
+    }
+    std::string output_path;
+#ifdef IS_WINDOWS
+    output_path = "D:/5_file/2_readfile/geojson/multi_point/fitting_examples_res_1.geojson";
+#else
+    output_path = "/home/wyj/myself/2_data/2_geojson/multi_point/fitting_examples_res_1_c.geojson";
+#endif
+    silly::geo::utils::write_geo_coll(output_path, res_collections);
+
+    silly::geo::utils::destroy_gdal_env();
+
+    return 0;
+}