gps坐标系转换--C实现

#include <math.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <stdio.h>

const double pi = 3.14159265358979324;
const double a = 6378245.0;
const double ee = 0.00669342162296594323;
const double x_pi = 3.14159265358979324 * 3000.0 / 180.0;

#define port_abs(x)     ((x > 0) ? (x) : ((x) * -1))

bool outOfChina(double lat, double lon)
{
    if (lon < 72.004 || lon > 137.8347)
        return true;
    if (lat < 0.8293 || lat > 55.8271)
        return true;
    return false;
}

double transformLat(double x, double y)
{
    double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(port_abs(x));
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * sin(y * pi) + 40.0 * sin(y / 3.0 * pi)) * 2.0 / 3.0;
    ret += (160.0 * sin(y / 12.0 * pi) + 320 * sin(y * pi / 30.0)) * 2.0 / 3.0;
    return ret;
}

double transformLon(double x, double y)
{
    double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(port_abs(x));
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * sin(x * pi) + 40.0 * sin(x / 3.0 * pi)) * 2.0 / 3.0;
    ret += (150.0 * sin(x / 12.0 * pi) + 300.0 * sin(x / 30.0 * pi)) * 2.0 / 3.0;
    return ret;
}

/* 地球坐标转换为火星坐标，即WGS84（国际通用）转为GCJ02坐标系适用于腾讯地图、高德（阿里）地图或谷歌地图 */
void WGS84toGCJ02(double wgLat, double wgLon, double *gcjLat, double *gcjLon)
{
    if (outOfChina(wgLat, wgLon)) {
        *gcjLat = wgLat;
        *gcjLon = wgLon;
        return;
    }
    double dLat = transformLat(wgLon - 105.0, wgLat - 35.0);
    double dLon = transformLon(wgLon - 105.0, wgLat - 35.0);
    double radLat = wgLat / 180.0 * pi;
    double magic = sin(radLat);
    magic = 1 - ee * magic * magic;
    double sqrtMagic = sqrt(magic);
    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
    dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * pi);
    *gcjLat = wgLat + dLat;
    *gcjLon = wgLon + dLon;
}

/* 地球坐标转换为百度坐标，即WGS84（国际通用）坐标系转为BD09坐标系适用于百度地图 */
void WGS84toBD09(double wgLat, double wgLon, double *bdLat, double *bdLon)
{
    double gcjLat, gcjLon;
    WGS84toGCJ02(wgLat, wgLon, &gcjLat, &gcjLon);

    double z = sqrt(gcjLon * gcjLon + gcjLat * gcjLat) + 0.00002 * sin(gcjLat * x_pi);
    double theta = atan2(gcjLat, gcjLon) + 0.000003 * cos(gcjLon * x_pi);
    *bdLat = z * sin(theta) + 0.006;    // #0.006     #0.01205
    *bdLon = z * cos(theta) + 0.0062;   // #0.0065    #0.00370
}

/* 火星坐标转换为百度坐标 */
void GCJ02toBD09(double gcjLat, double gcjLon, double *bdLat, double *bdLon)
{
    double x = gcjLon, y = gcjLat;
    double z = sqrt(x * x + y * y) + 0.00002 * sin(y * x_pi);
    double theta = atan2(y, x) + 0.000003 * cos(x * x_pi);
    *bdLon = z * cos(theta) + 0.0065;
    *bdLat = z * sin(theta) + 0.006;
}

/* 百度转火星 */
void BD09toGCJ02(double bdLat, double bdLon, double *gcjLat, double *gcjLon)
{
    double x = bdLon - 0.0065, y = bdLat - 0.006;
    double z = sqrt(x * x + y * y) - 0.00002 * sin(y * x_pi);
    double theta = atan2(y, x) - 0.000003 * cos(x * x_pi);
    *gcjLon = z * cos(theta);
    *gcjLat = z * sin(theta);
}


/*
GS84坐标系	121.469522,31.231905
GCJ02坐标系 121.47405857079794,31.229974621114046
BD09坐标系  121.48055338931275,31.235941661557085
*/
int main(void)
{
    double wgsLat = 31.231905;
    double wgsLon = 121.469522;

    double gcjLat = 0, gcjLon = 0;
    double bdLat = 0, bdLon = 0;

    WGS84toGCJ02(wgsLat, wgsLon, &gcjLat, &gcjLon);
    printf("WGS84toGCJ02:%lf-%lf >>> %lf-%lf\n", wgsLon, wgsLat, gcjLon, gcjLat);

    WGS84toBD09(wgsLat, wgsLon, &bdLat, &bdLon);
    printf("WGS84toBD09:%lf-%lf >>> %lf-%lf\n", wgsLon, wgsLat, bdLon, bdLat);

    GCJ02toBD09(gcjLat, gcjLon, &bdLat, &bdLon);
    printf("GCJ02toBD09:%lf-%lf >>> %lf-%lf\n", gcjLon, gcjLat, bdLon, bdLat);

    BD09toGCJ02(bdLat, bdLon, &gcjLat, &gcjLon);
    printf("BD09toGCJ02:%lf-%lf >>> %lf-%lf\n", bdLon, bdLat, gcjLon, gcjLat);

    return 0;
}

测试结果：

WGS84toGCJ02:121.469522-31.231905 >>> 121.474059-31.229975
WGS84toBD09:121.469522-31.231905 >>> 121.480298-31.235903
GCJ02toBD09:121.474059-31.229975 >>> 121.480598-31.235903
BD09toGCJ02:121.480598-31.235903 >>> 121.474058-31.229975

实际地理位置：