玩转GPS坐标系统

        在嵌入式开发中常会用到定位数据和地图功能，都会知道位置信息会有多个坐标系来表示，而不同地图厂商使用的坐标系也可能是不一样的，这就出现了在不同坐标系间转换坐标的需求。

应用需求：（1）一般GPS厂商提供的硬件设备，经过串口等工具获取的数据多为WGS84坐标系下的经纬度坐标。而我们使用的地图服务中，往往可能不是WGS84坐标系；（2）一些地图SDK给出的地图坐标转换接口的实现是在服务器进行的，那么这些接口的调用不但会发起网络请求而效率上也会非常的慢。以上两点，就需要本地实现坐标转换的算法来满足快速运算的需求。

1.GPS转换代码如下：

#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 >>> %.14lf-%.14lf\n", wgsLon, wgsLat, gcjLon, gcjLat);
 
    WGS84toBD09(wgsLat, wgsLon, &bdLat, &bdLon);
    printf("WGS84toBD09:%lf-%lf >>> %.14lf-%.14lf\n", wgsLon, wgsLat, bdLon, bdLat);
 
    GCJ02toBD09(gcjLat, gcjLon, &bdLat, &bdLon);
    printf("GCJ02toBD09:%lf-%lf >>> %.14lf-%.14lf\n", gcjLon, gcjLat, bdLon, bdLat);
 
    BD09toGCJ02(bdLat, bdLon, &gcjLat, &gcjLon);
    printf("BD09toGCJ02:%lf-%lf >>> %.14lf-%.14lf\n", bdLon, bdLat, gcjLon, gcjLat);
 
    return 0;
}

编译运行结果在线对比：

在线地图坐标系转换