UTM <-> WGS84 변환 C 코드

아래 github 저장소의 UTM <-> WGS84를 변환하는 Python 코드가 있어 이를 C 코드로 변환하는 코드로 바꾸었다.

https://github.com/Turbo87/utm

 

시험용 코드

테스트 코드는 몇 개의 위치를 표본으로 뽑아서 시험했으며 다음과 같다.

#include <stdio.h>
#define _USE_MATH_DEFINES
#include "utm_lla_translator.h"
#include <math.h>
int main()
{
    wgs84_coordinate wgs84s[] = {
        {50.775350,   6.08389}, // Aachen_wgs     
        {40.714350, -74.00597}, // NewYork_wgs    
        {-41.28646, 174.77624}, // Wellington_wgs 
        {-33.92487,  18.42406}, // Capetown_wgs   
        {-32.89018, -68.84405}, // Mendoza_wgs    
        {64.837780,-147.71639}, // Fairbanks_wgs  
        {56.796800,  -5.00601}, // BenNevis_wgs   
        {84.000000,  -5.00601}, // Latitude84_wgs 
    };
    utm_coordinate utms[] = {
        {294409.0, 5628898.0, 32, 'U'}, // Aachen_utm     
        {583960.0, 4507523.0, 18, 'T'}, // NewYork_utm    
        {313784.0, 5427057.0, 60, 'G'}, // Wellington_utm 
        {261878.0, 6243186.0, 34, 'H'}, // Capetown_utm   
        {514586.0, 6360877.0, 19, 'h'}, // Mendoza_utm    
        {466013.0, 7190568.0,  6, 'W'}, // Fairbanks_utm  
        {377486.0, 6296562.0, 30, 'V'}, // BenNevis_utm   
        {476594.0, 9328501.0, 30, 'X'}, // Latitude84_utm 
    };
    
    for (int i = 0; i < 8; i++){
        utm_coordinate now_utm = from_latlon(wgs84s[i].lat, wgs84s[i].lon);
        printf("LLA %12.6f %12.6f -> ", wgs84s[i].lat, wgs84s[i].lon);
        printf("[NEW]  %10.3f %10.3f %2d %c /",  now_utm.easting, now_utm.northing, now_utm.zone_number, now_utm.zone_letter);
        printf("[TRUE] %10.3f %10.3f %2d %c\n", utms[i].easting, utms[i].northing, utms[i].zone_number, utms[i].zone_letter);
    }
    for (int i = 0; i < 8; i++){
        wgs84_coordinate now_wgs = to_latlon(utms[i].easting, utms[i].northing, utms[i].zone_number, utms[i].zone_letter);
        printf("utm  %10.3f %10.3f %2d %c -> ",  utms[i].easting, utms[i].northing, utms[i].zone_number, utms[i].zone_letter);
        printf("[NEW]  %12.6f %12.6f / ", now_wgs.lat, now_wgs.lon);
        printf("[TRUE] %12.6f %12.6f \n", wgs84s[i].lat, wgs84s[i].lon);
    }
 
    printf("test \n");
    int i = 0;
    scanf("%d",&i);
}
 

 

 

 

변환 함수 코드

접은 글을 펼쳐보자.

// Copyright (C) 2012 Tobias Bieniek <Tobias.Bieniek@gmx.de>
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
 
#include <stdio.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include <ctype.h>
 
#ifndef _MAX
#define max(a,b) ((a>b) ? a : b)
#endif
#ifndef _MIN
#define min(a,b) ((a<b) ? a : b)
#endif
 
double wgs84_radius = 6378137.0;
char ZONE_LETTERS[] = "CDEFGHJKLMNPQRSTUVWXX"; // 22
 
typedef struct{
    double lat;
    double lon;
}wgs84_coordinate;
 
typedef struct{
    double easting;
    double northing;
    int zone_number;
    char zone_letter;
}utm_coordinate;
 
double mod_angle(double value){
    // return (value + mathlib.pi) % (2 * mathlib.pi) - mathlib.pi
    int mod = (value + M_PI)/(2  * M_PI);
    double rem_out = (value + M_PI) - mod * (2 * M_PI) - M_PI;
    return rem_out;
}
int is_bounds(double value, double lower, double upper){
    // printf("%f < %f < %f ???\n",lower, value, upper);
    return (lower <= value) && (value <= upper);
}
// int check_valid_zone(int zone_number, char zone_letter){
// }
 
char latitude_to_zone_letter(double latitude){
    // If the input is a numpy array, just use the first element
    // User responsibility to make sure that all points are in one zone
    if ((-80.0 <= latitude) && (latitude <= 84.0)){
        int tmp = (latitude + 80.0)/ 8;
        return ZONE_LETTERS[(min(tmp ,22))];
    }
    else
        printf("latittude is out of range \n");
}
 
int latlon_to_zone_number(double latitude, double longitude){
    // If the input is a numpy array, just use the first element
    // User responsibility to make sure that all points are in one zone
    int lat_range1 = (56.0 <= latitude) && (latitude < 64.0);
    int lon_range1 = (3.0 <= longitude) && (longitude < 12.0);
    if (lat_range1 && lon_range1)
        return 32;
 
    int lat_range2 = (72 <= latitude) && (latitude <= 84);
    int lon_range2 = (longitude >= 0);
    if (lat_range2 && lon_range2)
        if (longitude < 9)
            return 31;
        else if (longitude < 21)
            return 33;
        else if (longitude < 33)
            return 35;
        else if (longitude < 42)
            return 37;
 
 
    int zone_number = ((longitude + 180.0) / 6) + 1;
    return zone_number;
}
 
double zone_number_to_central_longitude(int zone_number){
    double lon = ((zone_number - 1) * 6 - 180 + 3);
    return lon;
}
 
wgs84_coordinate to_latlon(double easting, double northing, int zone_number, char zone_letter){
    double K0 = 0.9996;
 
    double E = 0.00669438;
    double E2 = E * E;
    double E3 = E2 * E;
    double E_P2 = E / (1 - E);
 
    double SQRT_E = sqrt(1 - E);
    double _E = (1 - SQRT_E) / (1 + SQRT_E);
    double _E2 = _E * _E;
    double _E3 = _E2 * _E;
    double _E4 = _E3 * _E;
    double _E5 = _E4 * _E;
 
    double M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256);
    double M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024);
    double M3 = (15 * E2 / 256 + 45 * E3 / 1024);
    double M4 = (35 * E3 / 3072);
 
    double P2 = (3 / 2 * _E - 27 / 32 * _E3 + 269 / 512 * _E5);
    double P3 = (21 / 16 * _E2 - 55 / 32 * _E4);
    double P4 = (151 / 96 * _E3 - 417 / 128 * _E5);
    double P5 = (1097 / 512 * _E4);
 
    double R2D = 180.0 / M_PI;
    double D2R = M_PI / 180.0;
    // check_valid_zone
    int northern = (toupper(zone_letter) >= 'N');
 
    double x = easting - 500000.0;
    double y = northing;
 
    if (!northern)
        y -= 10000000;
 
    double m = y / K0;
    double mu = m / (wgs84_radius * M1);
 
    double p_rad = (mu +
                    P2 * sin(2 * mu) +
                    P3 * sin(4 * mu) +
                    P4 * sin(6 * mu) +
                    P5 * sin(8 * mu));
 
    double p_sin = sin(p_rad);
    double p_sin2 = p_sin * p_sin;
 
    double p_cos = cos(p_rad);
 
    double p_tan = p_sin / p_cos;
    double p_tan2 = p_tan * p_tan;
    double p_tan4 = p_tan2 * p_tan2;
 
    double ep_sin = 1 - E * p_sin2;
    double ep_sin_sqrt = sqrt(1 - E * p_sin2);
 
    double n = wgs84_radius / ep_sin_sqrt;
    double r = (1 - E) / ep_sin;
 
    double c = E_P2 * p_cos * p_cos;
    double c2 = c * c;
 
    double d = x / (n * K0);
    double d2 = d * d;
    double d3 = d2 * d;
    double d4 = d3 * d;
    double d5 = d4 * d;
    double d6 = d5 * d;
 
    double latitude =   (p_rad - (p_tan / r) *
                        (d2 / 2 -
                        d4 / 24 * (5 + 3 * p_tan2 + 10 * c - 4 * c2 - 9 * E_P2)) +
                        d6 / 720 * (61 + 90 * p_tan2 + 298 * c + 45 * p_tan4 - 252 * E_P2 - 3 * c2));
 
    double longitude = (d -
                        d3 / 6 * (1 + 2 * p_tan2 + c) +
                        d5 / 120 * (5 - 2 * c + 28 * p_tan2 - 3 * c2 + 8 * E_P2 + 24 * p_tan4)) / p_cos;
 
    longitude = mod_angle(longitude + D2R * (zone_number_to_central_longitude(zone_number)));
 
    wgs84_coordinate wgs84_value = {0.0,0.0};
    wgs84_value.lat = R2D*latitude;
    wgs84_value.lon = R2D*longitude;
 
    return wgs84_value;
}
 
utm_coordinate from_latlon(double latitude, double longitude){
 
    double K0 = 0.9996;
 
    double E = 0.00669438;
    double E2 = E * E;
    double E3 = E2 * E;
    double E_P2 = E / (1 - E);
 
    double SQRT_E = sqrt(1 - E);
    double _E = (1 - SQRT_E) / (1 + SQRT_E);
    double _E2 = _E * _E;
    double _E3 = _E2 * _E;
    double _E4 = _E3 * _E;
    double _E5 = _E4 * _E;
 
    double M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256);
    double M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024);
    double M3 = (15 * E2 / 256 + 45 * E3 / 1024);
    double M4 = (35 * E3 / 3072);
 
    double P2 = (3 / 2 * _E - 27 / 32 * _E3 + 269 / 512 * _E5);
    double P3 = (21 / 16 * _E2 - 55 / 32 * _E4);
    double P4 = (151 / 96 * _E3 - 417 / 128 * _E5);
    double P5 = (1097 / 512 * _E4);
 
    double R2D = 180.0 / M_PI;
    double D2R = M_PI / 180.0;
    if(!is_bounds(latitude, -80.0, 84.0))
        printf("latitude out of range (must be between 80 deg S and 84 deg N)\n");
    if(!is_bounds(longitude, -180.0, 180.0))
        printf("longitude out of range (must be between 180 deg W and 180 deg E)\n");
 
    
    double lat_rad = D2R * (latitude);
    double lat_sin = sin(lat_rad);
    double lat_cos = cos(lat_rad);
 
    double lat_tan = lat_sin / lat_cos;
    double lat_tan2 = lat_tan * lat_tan;
    double lat_tan4 = lat_tan2 * lat_tan2;
 
    char zone_number = latlon_to_zone_number(latitude, longitude);
    int zone_letter = latitude_to_zone_letter(latitude);
 
    double lon_rad = D2R * (longitude);
    double central_lon = zone_number_to_central_longitude(zone_number);
    double central_lon_rad = D2R * (central_lon);
 
    double n = wgs84_radius / sqrt(1 - E * lat_sin*lat_sin);
    double c = E_P2 * lat_cos*lat_cos;
 
    double a = lat_cos * mod_angle(lon_rad - central_lon_rad);
    double a2 = a * a;
    double a3 = a2 * a;
    double a4 = a3 * a;
    double a5 = a4 * a;
    double a6 = a5 * a;
 
    double m = wgs84_radius * (M1 * lat_rad -
             M2 * sin(2 * lat_rad) +
             M3 * sin(4 * lat_rad) -
             M4 * sin(6 * lat_rad));
 
    double easting = K0 * n * (a +
                        a3 / 6 * (1 - lat_tan2 + c) +
                        a5 / 120 * (5 - 18 * lat_tan2 + lat_tan4 + 72 * c - 58 * E_P2)) + 500000;
 
    double northing = K0 * (m + n * lat_tan * (a2 / 2 +
                                        a4 / 24 * (5 - lat_tan2 + 9 * c + 4 * c*c) +
                                        a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2)));
 
 
    if (latitude < 0.0)
        northing += 10000000;
 
    utm_coordinate utm_value = {0,0,0,'A'};
    utm_value.easting = easting;
    utm_value.northing = northing;
    utm_value.zone_letter = zone_letter;
    utm_value.zone_number = zone_number;
 
    return utm_value;
}

 

** EOF **