如何使用Boost库的pj_geodetic_to_geocentric函数及传递Range参数？

这两个函数是Boost Geometry从PROJ4移植的坐标转换工具，用于**大地坐标（纬度、经度、椭球高）和地心坐标（X、Y、Z）**之间的相互转换，以下是具体使用细节：

参数解析 #

• T const& a：椭球的长半轴（单位通常为米，比如WGS84椭球的长半轴是6378137.0）
• T const& es：椭球的第一偏心率平方（WGS84椭球的该值为0.00669437999014，也可通过公式(f*(2-f))计算，其中f是扁率，WGS84扁率为1/298.257223563）
• Range &range：输入输出参数，用于传递待转换的坐标，转换后结果会直接覆盖原容器内的值。只要是符合Boost Range概念的容器都可以用，比如：
  • std::array<T, 3>（推荐，固定大小更安全）
  • std::vector<T>（需确保初始大小为3）
  • 原生C数组（如double[3]）

坐标顺序规则 #

• pj_geodetic_to_geocentric（大地转地心）：
  输入range的元素顺序为：[纬度（弧度）, 经度（弧度）, 椭球高（米）]
  转换后range的元素变为：[地心X, 地心Y, 地心Z]（单位与长半轴一致，通常为米）
• pj_geocentric_to_geodetic（地心转大地）：
  输入range的元素顺序为：[地心X, 地心Y, 地心Z]
  转换后range的元素变为：[纬度（弧度）, 经度（弧度）, 椭球高（米）]

代码示例 #

示例1：大地坐标转地心坐标 #

#include <array>
#include <iostream>
#include <cmath>
#include <boost/geometry/srs/projections/impl/pj_transform.hpp>

int main() {
    // WGS84椭球参数
    const double a = 6378137.0;
    const double es = 0.00669437999014;

    // 初始化大地坐标：40°N，116°E，椭球高100米（需转成弧度）
    std::array<double, 3> coords = {
        40.0 * M_PI / 180.0, 
        116.0 * M_PI / 180.0, 
        100.0
    };

    // 执行转换，返回0表示成功
    int status = pj_geodetic_to_geocentric(a, es, coords);

    if (status == 0) {
        std::cout << "转换后的地心坐标：" << std::endl;
        std::cout << "X: " << coords[0] << " m" << std::endl;
        std::cout << "Y: " << coords[1] << " m" << std::endl;
        std::cout << "Z: " << coords[2] << " m" << std::endl;
    } else {
        std::cerr << "转换失败，错误码：" << status << std::endl;
    }

    return 0;
}

示例2：地心坐标转大地坐标 #

#include <array>
#include <iostream>
#include <cmath>
#include <boost/geometry/srs/projections/impl/pj_transform.hpp>

int main() {
    const double a = 6378137.0;
    const double es = 0.00669437999014;

    // 初始化地心坐标（对应示例1的转换结果）
    std::array<double, 3> coords = {
        4482563.47, 
        4480450.85, 
        4078729.73
    };

    int status = pj_geocentric_to_geodetic(a, es, coords);

    if (status == 0) {
        // 将弧度转成角度输出
        std::cout << "转换后的大地坐标：" << std::endl;
        std::cout << "纬度: " << coords[0] * 180.0 / M_PI << "°" << std::endl;
        std::cout << "经度: " << coords[1] * 180.0 / M_PI << "°" << std::endl;
        std::cout << "椭球高: " << coords[2] << " m" << std::endl;
    } else {
        std::cerr << "转换失败，错误码：" << status << std::endl;
    }

    return 0;
}

注意事项 #

1. 输入的纬度、经度必须是弧度值，如果是角度需要先转换（角度 * M_PI / 180.0）
2. Range容器必须包含且仅包含3个元素，否则转换会出错
3. 返回值为0表示转换成功，非0值对应PROJ4定义的错误类型（比如坐标超出有效范围）

内容的提问来源于stack exchange，提问作者yuan