经纬度度分秒转换、距离计算、墨卡托投影及平面直角坐标转换

一、度分秒转换

度转换为度分秒、度分秒转换为度示例如下：

// 度转换为度分秒
std::string decimalToDMS(double decimalDeg) {
    int degrees = static_cast<int>(decimalDeg);
    double minutesSeconds = (decimalDeg - degrees) * 60;
    int minutes = static_cast<int>(minutesSeconds);
    double seconds = (minutesSeconds - minutes) * 60;
    return std::to_string(degrees) + "°" + std::to_string(minutes) + "'" + std::to_string(seconds) + "\"";
}

// 度分秒转换为度
double DMSToDecimal(const std::string& dms) {
    size_t degreePos = dms.find("°");
    size_t minutePos = dms.find("'");
    size_t secondPos = dms.find("\"");

    int degrees = std::stoi(dms.substr(0, degreePos));
    int minutes = std::stoi(dms.substr(degreePos + 1, minutePos - degreePos - 1));
    double seconds = std::stod(dms.substr(minutePos + 1, secondPos - minutePos - 1));

    double decimalDeg = degrees + (minutes / 60.0) + (seconds / 3600.0);
    return decimalDeg;
}

 

二、经纬度间距离计算

通用计算方法为Haversine公式，转换为代码如下：

double calculateDistance(double lat1, double lon1, double lat2, double lon2) {
    double earthRadius = 6371.0; // 地球半径，单位：千米

    // 将经纬度转换为弧度
    double dLat = (lat2 - lat1) * M_PI / 180.0;
    double dLon = (lon2 - lon1) * M_PI / 180.0;

    // 应用Haversine公式计算直线距离
    double a = sin(dLat / 2) * sin(dLat / 2) + cos(lat1 * M_PI / 180.0) * cos(lat2 * M_PI / 180.0) * sin(dLon / 2) * sin(dLon / 2);
    double c = 2 * atan2(sqrt(a), sqrt(1 - a));
    double distance = earthRadius * c;

    return distance;
}

 

三、经纬度墨卡托投影

地理投影方法之一，可以将地球表面的经纬度坐标投影到平面坐标系上
原理参考：https://dandelioncloud.cn/article/details/1573648470961778689
转换为代码如下：

void mercatorProjection(double lat, double lon, double& x, double& y) {
    const double R = 6378137.0; // 地球半径，单位：米

    double lat_rad = lat * M_PI / 180.0; // 将纬度转换为弧度
    double lon_rad = lon * M_PI / 180.0; // 将经度转换为弧度

    x = R * lon_rad; // 经度直接转换为 X 坐标
    y = R * log(tan(M_PI / 4 + lat_rad / 2)); // 纬度转换为 Y 坐标
}

缺点：在墨卡托投影中，纬度值大致等于真实地球上的距离，而经度值仅作线性缩放。这种投影方法适用于大范围的地理数据处理。
 

四、经纬度平面直角坐标转换

假设地球是一个完美的球体，并将经纬度转换为直角坐标系中的X和Y坐标。然后，根据输入的目标航向角对坐标进行旋转，以得到参照物相对于目标的坐标偏移。
转换为代码如下：

const double EARTH_RADIUS = 6378137.0; // 地球半径，单位：米

struct Coordinate {
    double x; // 参照物在目标坐标系中的X坐标
    double y; // 参照物在目标坐标系中的Y坐标
};

Coordinate calculateCoordinates(double masterLatitude, double masterLongitude, double slaveLatitude, double slaveLongitude, double heading) {
    double masterRadLat = masterLatitude * M_PI / 180;
    double masterRadLon = masterLongitude * M_PI / 180;
    double slaveRadLat = slaveLatitude * M_PI / 180;
    double slaveRadLon = slaveLongitude * M_PI / 180;

    double latDifference = slaveRadLat - masterRadLat;
    double lonDifference = slaveRadLon - masterRadLon;

	//乘以 EARTH_RADIUS 来计算参照物相对于目标的东西方向（水平方向）的偏移量,由于纬度变化会导致经度方向上的弧长缩短，所以需要乘以 cos(masterRadLat)将纬度的影响纳入考虑，从而得到更准确的结果
    double dx = EARTH_RADIUS * lonDifference * std::cos(masterRadLat);
	//乘以 EARTH_RADIUS 来计算参照物相对于目标的南北方向（垂直方向）的偏移量
    double dy = EARTH_RADIUS * latDifference;

    // 考虑目标航向角，heading 为正,对坐标进行逆时针旋转
    double radHeading = heading * M_PI / 180;
	double rotatedX = dx * std::cos(radHeading) - dy * std::sin(radHeading);
    double rotatedY = dx * std::sin(radHeading) + dy * std::cos(radHeading);

	//  若顺时针旋转
	//double rotatedX = dx * std::cos(radHeading) + dy * std::sin(radHeading);
    //double rotatedY = -dx * std::sin(radHeading) + dy * std::cos(radHeading);
	
    Coordinate coord;
    coord.x = rotatedX;
    coord.y = rotatedY;

    return coord;
}

缺点：实际运用需对参照物运动模型进行预测和做滤波处理。