utm坐标和经纬度在线转换的C++实现

以下是UTM坐标和经纬度在线转换的C++实现示例：

#include <stdio.h>
#include <math.h>

#define PI 3.14159265

void UTMtoLL(double UTMNorthing, double UTMEasting, const char* UTMZone, double& Lat, double& Long)
{
    double k0 = 0.9996;
    double a = 6378137;
    double eccSquared = 0.00669438;
    double eccPrimeSquared;
    double e1 = (1 - sqrt(1 - eccSquared)) / (1 + sqrt(1 - eccSquared));
    double N1, T1, C1, R1, D, M;
    double LongOrigin;
    double mu, phi1, phi1Rad;
    double x, y;

    // remove 500,000 meter offset for longitude
    x = UTMEasting - 500000.0;

    // We must know somehow if we are in the Northern or Southern
    // hemisphere, this is the only time we use the letter so we can
    // assume that it is in North
    char hemisphere = 'N'; 

    // Make sure the zone number is between 1 and 60
    int zoneNumber = atoi(UTMZone);
    if ((zoneNumber < 1) || (zoneNumber > 60))
    {
        printf("Invalid UTM zone number: %d\n", zoneNumber);
        return;
    }

    LongOrigin = (zoneNumber - 1) * 6 - 180 + 3;  // +3 puts origin in middle of zone

    eccPrimeSquared = (eccSquared) / (1 - eccSquared);

    M = UTMNorthing / k0;
    mu = M / (a * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256));

    phi1Rad = 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);
    phi1 = phi1Rad * 180.0 / PI;

    N1 = a / sqrt(1 - eccSquared * sin(phi1Rad) * sin(phi1Rad));
    T1 = tan(phi1Rad) * tan(phi1Rad);
    C1 = eccPrimeSquared * cos(phi1Rad) * cos(phi1Rad);
    R1 = a * (1 - eccSquared) / pow(1 - eccSquared * sin(phi1Rad) * sin(phi1Rad), 1.5);
    D = x / (N1 * k0);

    Lat = phi1Rad - (N1 * tan(phi1Rad) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * eccPrimeSquared) * D * D * D * D / 24 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * eccPrimeSquared - 3 * C1 * C1) * D * D * D * D * D * D / 720);
    Lat = Lat * 180.0 / PI;

    Long = (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * eccPrimeSquared + 24 * T1 * T1) * D * D * D * D * D / 120) / cos(phi1Rad);
    Long = LongOrigin + Long * 180.0 / PI;

    printf("Latitude: %f Longitude: %f\n", Lat, Long);
}

void LLtoUTM(double Lat, double Long, double& UTMNorthing, double& UTMEasting, char* UTMZone)
{
    double k0 = 0.9996;
    double a = 6378137;
    double eccSquared = 0.00669438;
    double eccPrimeSquared;
    double N, T, C, A, M;
    double LatRad = Lat * PI / 180.0;
    double LongRad = Long * PI / 180.0;
    int zoneNumber;
    char hemisphere;

    if ((Long >= -180.0) && (Long < -174.0))
        zoneNumber = 1;
    else if ((Long >= -174.0) && (Long < -168.0))
        zoneNumber = 2;
    else if ((Long >= -168.0) && (Long < -162.0))
        zoneNumber = 3;
    else if ((Long >= -162.0) && (Long < -156.0))
        zoneNumber = 4;
    else if ((Long >= -156.0) && (Long < -150.0))
        zoneNumber = 5;
    else if ((Long >= -150.0) && (Long < -144.0))
        zoneNumber = 6;
    else if ((Long >= -144.0) && (Long < -138.0))
        zoneNumber = 7;
    else if ((Long >= -138.0) && (Long < -132.0))
        zoneNumber = 8;
    else if ((Long >= -132.0) && (Long < -126.0))
        zoneNumber = 9;
    else if ((Long >= -126.0) && (Long < -120.0))
        zoneNumber = 10;
    else if ((Long >= -120.0) && (Long < -114.0))
        zoneNumber = 11;
    else if ((Long >= -114.0) && (Long < -108.0))
        zoneNumber = 12;
    else if ((Long >= -108.0) && (Long < -102.0))
        zoneNumber = 13;
    else if ((Long >= -102.0) && (Long < -96.0))
        zoneNumber = 14;
    else if ((Long >= -96.0) && (Long < -90.0))
        zoneNumber = 15;
    else if ((Long >= -90.0) && (Long < -84.0))
        zoneNumber = 16;
    else if ((Long >= -84.0) && (Long < -78.0))
        zoneNumber = 17;
    else if ((Long >= -78.0) && (Long < -72.0))
        zoneNumber = 18;
    else if ((Long >= -72.0) && (Long < -66.0))
        zoneNumber = 19;
    else if ((Long >= -66.0) && (Long < -60.0))
        zoneNumber = 20;
    else if ((Long >= -60.0) && (Long < -54.0))
        zoneNumber = 21;
    else if ((Long >= -54.0) && (Long < -48.0))
        zoneNumber = 22;
    else if ((Long >= -48.0) && (Long < -42.0))
        zoneNumber = 23;
    else if ((Long >= -42.0) && (Long < -36.0))
        zoneNumber = 24;
    else if ((Long >= -36.0) && (Long < -30.0))
        zoneNumber = 25;
    else if ((Long >= -30.0) && (Long < -24.0))
        zoneNumber = 26;
    else if ((Long >= -24.0) && (Long < -18.0))
        zoneNumber = 27;
    else if ((Long >= -18.0) && (Long < -12.0))
        zoneNumber = 28;
    else if ((Long >= -12.0) && (Long < -6.0))
        zoneNumber = 29;
    else if ((Long >= -6.0) && (Long < 0.0))
        zoneNumber = 30;
    else if ((Long >= 0.0) && (Long < 6.0))
        zoneNumber = 31;
    else if ((Long >= 6.0) && (Long < 12.0))
        zoneNumber = 32;
    else if ((Long >= 12.0) && (Long < 18.0))
        zoneNumber = 33;
    else if ((Long >= 18.0) && (Long < 24.0))
        zoneNumber = 34;
    else if ((Long >= 24.0) && (Long < 30.0))
        zoneNumber = 35;
    else if ((Long >= 30.0) && (Long < 36.0))
        zoneNumber = 36;
    else if ((Long >= 36.0) && (Long < 42.0))
        zoneNumber = 37;
    else if ((Long >= 42.0) && (Long < 48.0))
        zoneNumber = 38;
    else if ((Long >= 48.0) && (Long < 54.0))
        zoneNumber = 39;
    else if ((Long >= 54.0) && (Long < 60.0))
        zoneNumber = 40;
    else
    {
        printf("Invalid longitude: %f\n", Long);
        return;
    }

    if (Lat >= 0)
        hemisphere = 'N';
    else
        hemisphere = 'S';

    eccPrimeSquared = (eccSquared) / (1 - eccSquared);

    N = a / sqrt(1 - eccSquared * sin(LatRad) * sin(LatRad));
    T = tan(LatRad) * tan(LatRad);
    C = eccPrimeSquared * cos(LatRad) * cos(LatRad);
    A = cos(LatRad) * (LongRad - (zoneNumber - 1) * 6 * PI / 180.0);
    M = a * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatRad
        - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(2 * LatRad)
        + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(4 * LatRad)
        - (35 * eccSquared * eccSquared * eccSquared / 3072) * sin(6 * LatRad));

    UTMEasting = (double)(k0 * N * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120) + 500000.0);
    UTMNorthing = (double)(k0 * (M + N * tan(LatRad) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A * A * A * A * A * A / 720)));

    if (hemisphere == 'S')
        UTMNorthing += 10000000.0; // 10000000 meter offset for southern hemisphere

    sprintf(UTMZone, "%d%c", zoneNumber, hemisphere);
    printf("UTM Northing: %f UTM Easting: %f UTM Zone: %s\n", UTMNorthing, UTMEasting, UTMZone);
}

使用示例：

int main()
{
    double UTMNorthing, UTMEasting, Lat, Long;
    char UTMZone[10];

    // UTM to LL
    UTMNorthing = 4383203.180;
    UTMEasting = 343409.350;
    strcpy(UTMZone, "32T");
    UTMtoLL(UTMNorthing, UTMEasting, UTMZone, Lat, Long);

    // LL to UTM
    Lat = 52.520007;
    Long = 13.404954;
    LLtoUTM(Lat, Long, UTMNorthing, UTMEasting, UTMZone);

    return 0;
}