UTM в WGS84 в C

Я пытаюсь преобразовать координаты из UTM в WGS84 / (LatLon) на языке C.

Я нашел эту похожую тему: Как преобразовать из UTM в LatLng в Python или Javascript

Однако при преобразовании кода широта верна, а долгота - нет.

Код

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

int main(int argc, char **argv)
{

// Input
double east = 724197;           // Easting
double north = 6175457;         // Northing
double utmZone = 32;            // UTM zone
double FalseNorth = 0;          // Northern hemisphere:  FalseNorth = 0;
                                // Southern hemisphere:  FalseNorth = 10000000;
// Constants (WGS84 datum)
double SemiMajor = 6378137;     // Ellipsoidal semi-major axis (Meters)
double FalseEast = 500000;      // UTM East bias (Meters)
double ScaleFactor = 0.9996;    // Scale at natural origin

double LngOrigin = utmZone * 6 - 183;   // Reference meridian of longitude

double Ecc = 0.081819190842622;         // Eccentricity
double EccSq = Ecc * Ecc;
double Ecc2Sq = EccSq / (1. - EccSq);
double Ecc2 = sqrt(Ecc2Sq);             // Secondary eccentricity
double E1 = ( 1 - sqrt(1-EccSq) ) / ( 1 + sqrt(1-EccSq) );
double E12 = E1 * E1;
double E13 = E12 * E1;
double E14 = E13 * E1;

// Output
double result_lat;            // Latitude in WGS84
double result_lon;            // Longitude in WGS84

// Calculate the Cassini projection parameters
double M1 = (north - FalseNorth) / ScaleFactor;
double Mu1 = M1 / ( SemiMajor * (1 - EccSq/4.0 - 3.0*EccSq*EccSq/64.0 -
                    5.0*EccSq*EccSq*EccSq/256.0) );

double Phi1 = Mu1   + (3.0*E1/2.0 - 27.0*E13/32.0)      * sin(2.0*Mu1)
                    + (21.0*E12/16.0 - 55.0*E14/32.0)   * sin(4.0*Mu1)
                    + (151.0*E13/96.0)                  * sin(6.0*Mu1)
                    + (1097.0*E14/512.0)                * sin(8.0*Mu1);

double sin2phi1 = sin(Phi1) * sin(Phi1);
double Rho1 = (SemiMajor * (1.0-EccSq) ) / pow(1.0-EccSq*sin2phi1,1.5);
double Nu1 = SemiMajor / sqrt(1.0-EccSq*sin2phi1);

// Compute parameters as defined in the POSC specification.  T, C and D
double T1 = tan(Phi1) * tan(Phi1);
double T12 = T1 * T1;
double C1 = Ecc2Sq * cos(Phi1) * cos(Phi1);
double C12 = C1 * C1;
double D  = (east - FalseEast) / (ScaleFactor * Nu1);
double D2 = D * D;
double D3 = D2 * D;
double D4 = D3 * D;
double D5 = D4 * D;
double D6 = D5 * D;

// Compute the Latitude and Longitude and convert to degrees
double lat = Phi1 - Nu1*tan(Phi1)/Rho1 *
                    ( D2/2.0 - (5.0 + 3.0*T1 + 10.0*C1 - 4.0*C12 - 9.0*Ecc2Sq)*D4/24.0
                    + (61.0 + 90.0*T1 + 298.0*C1 + 45.0*T12 - 252.0*Ecc2Sq - 3.0*C12)*D6/720.0);

result_lat = (lat * 180.0 / M_PI);

double lon = LngOrigin + 
                    ( D - (1.0 + 2.0*T1 + C1)*D3/6.0
                    + (5.0 - 2.0*C1 + 28.0*T1 - 3.0*C12 + 8.0*Ecc2Sq + 24.0*T12)*D5/120.0) / cos(Phi1);

result_lon = (lon * 180.0 / M_PI);

printf("result_lat = %lf\n", result_lat);
printf("result_lon = %lf\n", result_lon);

return 0;
}

Результат

result_lat = 55.673059
result_lon = 519.227573

Я ожидал, что result_lon будет близко к 12.56


c utm wgs84
person Magnus Olsen    schedule 14.07.2018    source источник
comment
printf %lf эквивалентно просто %f в C99. В C89 он имеет неопределенное поведение.   -  person melpomene    schedule 14.07.2018
comment
Какие у вас есть основания полагать, что исходный код, который вы перевели, был правильным? Вы запускали этот код JavaScript, чтобы проверить, дал ли он ожидаемые результаты?   -  person Eric Postpischil    schedule 14.07.2018

Ответы (1)


arrow_upward
0
arrow_downward

Вы также можете использовать классический PROJ API (сейчас устарел), показанный здесь:

#include <stdio.h>
#include <proj_api.h>

int main(int argc, char **argv)
{
    projPJ pj_utm, pj_latlong;
    char utm_s[50];
    double x = 724197.0;
    double y = 6175457.0;
    int zone = 32;
    sprintf(utm_s, "+proj=utm +zone=%d", zone);
    pj_utm = pj_init_plus(utm_s);
    pj_latlong = pj_init_plus("+proj=latlong");
    pj_transform(pj_utm, pj_latlong, 1, 1, &x, &y, NULL );
    x *= RAD_TO_DEG;
    y *= RAD_TO_DEG;
    printf("result (lat, lon) = (%f, %f)\n", y, x);
    pj_free(pj_utm);
    pj_free(pj_latlong);
    return 0;
}

Скомпилируйте, свяжите (при условии, что у вас где-то есть libproj.so) и запустите:

$ gcc code.c -lproj
$ ./a.out
result (lat, lon) = (55.673059, 12.565558)
person Mike T    schedule 10.08.2018