WGS84、BD09、GCJ02座標轉換

汗牛充栋發表於2024-09-04
GC

名詞解釋

WGS84

此座標系解釋參考筆者另一篇部落格GIS座標系、投影與轉換

GCJ02

GCJ-02是由中國國家測繪局（G表示Guojia國家，C表示Cehui測繪，J表示Ju局）制訂的地理資訊系統的座標系統。 是中國大陸地區的地圖資料使用的座標系。基於 WGS84 進行加密後形成。

BD09

BD09 是中國的百度地圖使用的一種座標系，全稱是百度座標系（BD-09）。它是在 GCJ-02（國測局座標系）基礎上進行的進一步加密處理。百度地圖使用 BD-09 座標系，這種座標系主要用於在百度地圖上的定位和導航。

座標轉換

WGS84轉GCJ02

需要經過一個偏移演算法

double offset_gcj02_latitude(double x, double y) {
    double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * std::sqrt(abs(x));
    ret += (20.0 * std::sin(6.0 * x * pi) + 20.0 * std::sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * std::sin(y * pi) + 40.0 * std::sin(y / 3.0 * pi)) * 2.0 / 3.0;
    ret += (160.0 * std::sin(y / 12.0 * pi) + 320 * std::sin(y * pi / 30.0)) * 2.0 / 3.0;
    return ret;
}

double offset_gcj02_longitude(double x, double y) {
    double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * std::sqrt(abs(x));
    ret += (20.0 * std::sin(6.0 * x * pi) + 20.0 * std::sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * std::sin(x * pi) + 40.0 * std::sin(x / 3.0 * pi)) * 2.0 / 3.0;
    ret += (150.0 * std::sin(x / 12.0 * pi) + 300.0 * std::sin(x / 30.0 * pi)) * 2.0 / 3.0;
    return ret;
}

void gcj02_delta(double lat, double lon, double& dLat, double& dLon) {
    const double a = 6378245.0;                 // 長半軸
    const double ee = 0.00669342162296594323;   // 扁率

    double radLat = lat / 180.0 * pi;
    double magic = std::sin(radLat);
    magic = 1 - ee * magic * magic;
    double sqrtMagic = std::sqrt(magic);
    dLat = (offset_gcj02_latitude(lon - 105.0, lat - 35.0) * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
    dLon = (offset_gcj02_longitude(lon - 105.0, lat - 35.0) * 180.0) / (a / sqrtMagic * std::cos(radLat) * pi);
}

void wgs84_to_gcj02(double wgsLat, double wgsLon, double& gcjLat, double& gcjLon) {
    if (outOfChina(wgsLat, wgsLon)) {
        gcjLat = wgsLat;
        gcjLon = wgsLon;
        return;
    }

    double lat_offset=0, lon_offset=0;
    gcj02_delta(wgsLat, wgsLon, lat_offset, lon_offset);
    gcjLat = wgsLat + lat_offset;
    gcjLon = wgsLon + lon_offset;
}

如上述程式碼所示， 其轉換僅需在WGS84座標上加上對應偏移即可。 因為是中國國測局制定的， 所以在轉換開始判斷座標是否位於中國大陸境內。

逆轉的話， 減去對應的偏移即可

WGS84轉BD09

WGS84轉為BD09需首先將其轉換為GCJ02, 然後再在GCJ02基礎上作加密偏移。以下程式碼僅僅示例瞭如何從GCJ02轉換到BD09

void gcj02_to_bd09(double gcjLat, double gcjLon, double& bdLat, double& bdLon) {
    double z = std::sqrt(gcjLon * gcjLon + gcjLat * gcjLat) + 0.00002 * std::sin(gcjLat * pi);
    double theta = std::atan2(gcjLat, gcjLon) + 0.000003 * std::cos(gcjLon * pi);
    bdLon = z * std::cos(theta) + 0.0065;
    bdLat = z * std::sin(theta) + 0.006;
}

其逆轉也是一個反向過程

void bd09_to_gcj02(double bdLat, double bdLon, double& gcjLat, double& gcjLon) {
    double x = bdLon - 0.0065;
    double y = bdLat - 0.006;
    double z = std::sqrt(x * x + y * y) - 0.00002 * std::sin(y * pi);
    double theta = std::atan2(y, x) - 0.000003 * std::cos(x * pi);
    gcjLon = z * std::cos(theta);
    gcjLat = z * std::sin(theta);
}

轉換為笛卡爾座標系

由於GCJ02與BD09都是基於WGS84的偏移加密， 所以其大地座標系的基準都是與WGS84一致。所以轉換為笛卡爾座標系需要透過WGS84進行轉換。

WGS84 大地座標系轉換為 ECEF 空間直角座標系的原理涉及幾何學和橢球體的數學描述。簡單來說，轉換過程是透過將地球表面的點（經緯度和高度）轉換為以地球中心為原點的三維直角座標（X, Y, Z）。

轉換的原理核心在於透過幾何公式將橢球上的點投影到三維空間中的一個點。ECEF 座標系統與地球同步旋轉，因此適合用於全球定位系統（GPS）和其他涉及地球表面點的位置計算的系統。這種轉換能夠在全球範圍內提供一個統一的三維座標系，使得不同地點的測量資料能夠準確地相互比較和計算。

以下是這個轉換的詳細原理：

1. 橢球模型描述

  • WGS84 橢球模型：地球並非完美的球體，而是一個扁球體，赤道半徑大於極半徑。WGS84 使用一個橢球來近似地球的形狀：
    • 長半軸 (a)：赤道的半徑，大約為 6378137 米。
    • 短半軸 (b)：極點的半徑，大約為 6356752 米。
    • 扁率 (f)：描述了橢球的扁平程度，計算公式為 \(f = \frac{a - b}{a}\)。

2. 大地座標系

  • 緯度 (lat)：從地心到橢球表面點的法線與赤道平面的夾角。
  • 經度 (lon)：地球表面點所在的子午線與本初子午線（經度為 0° 的子午線）之間的夾角。
  • 高度 (h)：地球表面點沿法線方向相對於橢球面的距離。

3. 空間直角座標系 (ECEF)

  • ECEF 系統：這個座標系的原點位於地球中心，X 軸指向經度為 0° 的赤道上，Y 軸指向經度為 90° 的赤道上，Z 軸指向北極。

4. 轉換公式推導

為了將大地座標轉換為 ECEF 座標，需理解以下幾個概念：

曲率半徑 (N)
  • 對於緯度為 \(\text{lat}\) 的點，曲率半徑 \(N\) 表示的是橢球體表面上的一個點到橢球中心的距離，定義為：

\[ N = \frac{a}{\sqrt{1 - e^2 \sin^2(\text{lat})}} \]

其中，\(e^2 = \frac{2f - f^2}{1}\) 是橢球的第一偏心率的平方，描述了橢球的形狀。

ECEF 座標公式
  • 使用 \(N\) 和地心座標系的三維直角座標，X、Y、Z 計算如下：

    \[X = (N + h) \cos(\text{lat}) \cos(\text{lon}) \]

    \[Y = (N + h) \cos(\text{lat}) \sin(\text{lon}) \]

    \[Z = \left[\left(1 - e^2\right)N + h\right] \sin(\text{lat}) \]

這些公式從幾何上來講是透過將地球表面的一個點投影到地心座標系中來實現的。X 和 Y 是根據緯度和經度的餘弦值來確定，Z 則根據緯度的正弦值計算。

5. 轉換過程

  1. 計算 N：根據緯度來計算曲率半徑 \(N\)，這涉及到 WGS84 橢球的偏心率和緯度的三角函式。
  2. 計算 ECEF 座標：使用上面推匯出的公式，透過緯度、經度和高度來計算該點在地心座標系中的 X, Y, Z 座標。

6. 程式碼示例

void wgs84_to_ecef(double lat, double lon, double alt, double& x, double& y, double& z) {
    const double a = 6378137.0;             // WGS-84 橢球體的長半軸（赤道半徑）
    const double f = 1 / 298.257223563;     // 地球扁率
    const double e2 = 2 * f - f * f;        // 0.00669437999014;  WGS-84 橢球體的第一偏心率的平方

    double radLat = lat * M_PI / 180.0;     // 將緯度轉換為弧度
    double radLon = lon * M_PI / 180.0;     // 將經度轉換為弧度

    double N = a / std::sqrt(1 - e2 * std::sin(radLat) * std::sin(radLat));  // 橢球曲率半徑

    x = (N + alt) * std::cos(radLat) * std::cos(radLon);
    y = (N + alt) * std::cos(radLat) * std::sin(radLon);
    z = (N * (1 - e2) + alt) * std::sin(radLat);
}

逆轉:

void ecef_to_wgs84(double x, double y, double z, double& lat, double& lon, double& alt) {

    const double a = 6378137.0;                 // WGS-84 橢球體的長半軸（赤道半徑）
    const double f = 1.0 / 298.257223563;       // 扁率
    const double e2 = 2 * f - f * f;            //  0.00669437999014;  WGS-84 橢球體的第一偏心率的平方
    const double b = a * (1 - f);               // 6356752.314245; WGS-84 橢球體的短半軸（極半徑）


    double eps = 1e-12;   // 精度
    double p = std::sqrt(x * x + y * y);
    lon = std::atan2(y, x);    // 計算經度

    // 迭代計算緯度
    double theta = std::atan2(z * a, p * b);
    double sinTheta = std::sin(theta);
    double cosTheta = std::cos(theta);

    lat = std::atan2(z + e2 * b * sinTheta * sinTheta * sinTheta,
        p - e2 * a * cosTheta * cosTheta * cosTheta);

    double N = a / std::sqrt(1 - e2 * std::sin(lat) * std::sin(lat));
    alt = p / std::cos(lat) - N;

    lat = lat * 180.0 / M_PI;   // 轉換為度數
    lon = lon * 180.0 / M_PI;   // 轉換為度數
}

相關文章