GPS数据：从读取、解析到转换

最近因为一个小课题开始研究GPS模块数据的处理，要想拿到最后实际可用的GPS数据信息，从数据读取、数据解析以及数据转换都需要进行相应的研究。
在本文中将简要如何通过STM32读取GPS模块原始数据，并对数据进行解析处理得到可视的GPS数据、最终通过坐标系转换得到在地图上可用的GPS信息。

GPS模块介绍

GPS接收器

颜色	功能
红色	VCC
黑色	GND
绿色	TxD
黄色	RxD

GPS数据介绍

GPS上电后，每隔一定的时间就会返回一定格式的数据，该数据遵循NMEA0183标准，
该协议采用ASCII码，其串行通信默认参数为：波特率=9600bps，数据位=8bit，开始位=1bit，停止位=1bit，无奇偶校验。
GPS数据格式为：

$aaccc,ddd,ddd,...,ddd*hh<CR><LF>

$: 帧命令起始位
aaccc: 地址域，前两位为识别符，后三位为语句名
ddd,...,ddd: 数据，数据间以逗号分隔，若有部分数据缺失，数据对应位置留空，保留逗号分隔符。
*: 校验和前缀
hh: 校验和，$与*之间所有字符ASCII码的校验和
<CR><LF>: 帧结束标志，回车符和换行符

以下是示例数据：

$GPRMC,092927.000,A,2235.9058,N,11400.0518,E,0.000,74.11,151216,,D*49
$GPVTG,74.11,T,,M,0.000,N,0.000,K,D*0B
$GPGGA,092927.000,2235.9058,N,11400.0518,E,2,9,1.03,53.1,M,-2.4,M,0.0,0*6B
$GPGSA,A,3,29,18,12,25,10,193,32,14,31,,,,1.34,1.03,0.85*31
$GPGSV,3,1,12,10,77,192,17,25,59,077,42,32,51,359,39,193,49,157,36*48
$GPGSV,3,2,12,31,47,274,25,50,46,122,37,18,45,158,37,14,36,326,18*70
$GPGSV,3,3,12,12,24,045,45,26,17,200,18,29,07,128,38,21,02,174,*79

GPGGA GPS定位数据

Global Positioning System Fix Data

$GPGGA,<1>,<2>,<3>,<4>,<5>,<6>,<7>,<8>,<9>,M,<10>,M,<11>,<12>*hh

<1> UTC时间，格式为hhmmss.sss(时分秒)
<2> 纬度，格式为ddmm.mmmm(首位补零)
<3> 纬度半球，N=北纬，S=南纬
<4> 经度，格式为dddmm.mmmm(首位补零)
<5> 经度半球，E=东经，W=西经
<6> 定位质量指示，0=定位无效，1=定位有效
<7> 卫星数量，00~12(首位补零)
<8> 水平精确度，0.5~99.9
<9> 天线离海平面的高度，-9999.9~9999.9米，M指单位米
<10> 大地水准面高度，-9999.9~9999.9米，M指单位米
<11> 差分GPS数据期限(RTCMSC-104)，最后设立RTCM传送的秒数量
<12> 差分参考基站标号，0000~1023(首位补零)

GPRMC 最小定位信息

Recommended Minimum Specific GPS/TRANSIT Data

$GPRMC,<1>,<2>,<3>,<4>,<5>,<6>,<7>,<8>,<9>,<10>,<11>,<12>*hh

<1> UTC时间，格式为hhmmss.sss(时分秒)
<2> 定位状态，A=有效定位，V=无效定位
<3> 纬度，格式为ddmm.mmmm(度分)(首位补零)
<4> 纬度半球，N=北纬，S=南纬
<5> 经度，格式为dddmm.mmmm(度分)(首位补零)
<6> 经度半球，E=东经，W=西经
<7> 地面速度，000.0~999.9节(首位补零)
<8> 地面航向，000.0~359.9度，以真北为参考基准(首位补零)
<9> UTC日期，格式为ddmmyy(日月年)
<10> 磁偏角，000.0~180.0度(首位补零)
<11> 磁偏角方向，E=东，W=西
<12> 模式指示(仅NMEA0183-V3.00输出，A=自主定位，D=差分，E=估算，N=数据无效)

其他信息

• GPVTG 地面速度信息
• GPGSV 可视卫星状态
• GPGSA 当前卫星信息

GPS数据读取

配置串口

STM32CubeMX自动生成配置代码

UART_HandleTypeDef huart2;
void MX_USART2_UART_Init(void)
{
  huart2.Instance = USART2;
  huart2.Init.BaudRate = 9600;
  huart2.Init.WordLength = UART_WORDLENGTH_8B;
  huart2.Init.StopBits = UART_STOPBITS_1;
  huart2.Init.Parity = UART_PARITY_NONE;
  huart2.Init.Mode = UART_MODE_TX_RX;
  huart2.Init.HwFlowCtl = UART_HWCONTROL_NONE;
  huart2.Init.OverSampling = UART_OVERSAMPLING_16;
  if (HAL_UART_Init(&huart2) != HAL_OK)
  {
    _Error_Handler(__FILE__, __LINE__);
  }
}

void HAL_UART_MspInit(UART_HandleTypeDef* uartHandle)
{
  GPIO_InitTypeDef GPIO_InitStruct;
  if(uartHandle->Instance==USART2)
  {
    __HAL_RCC_USART2_CLK_ENABLE();
  
    GPIO_InitStruct.Pin = GPIO_PIN_2|GPIO_PIN_3;
    GPIO_InitStruct.Mode = GPIO_MODE_AF_PP;
    GPIO_InitStruct.Pull = GPIO_PULLUP;
    GPIO_InitStruct.Speed = GPIO_SPEED_FREQ_VERY_HIGH;
    GPIO_InitStruct.Alternate = GPIO_AF7_USART2;
    HAL_GPIO_Init(GPIOA, &GPIO_InitStruct);

    HAL_NVIC_SetPriority(USART2_IRQn, 0, 0);
    HAL_NVIC_EnableIRQ(USART2_IRQn);
  }
}

void HAL_UART_RxCpltCallback(UART_HandleTypeDef *huart)
{
    if(huart->Instance == USART2)
    {
        length = GPS_MsgRecv(aRxBuffer);
        if(length)  //这一部分也可以直接放在main函数中进行判断
        {
            memset(printBuffer,0,80);
            memcpy(printBuffer,rmcBuffer,length);
            decodeFlag = length;
        }
        HAL_UART_Receive_IT(&huart2, &aRxBuffer, 1);
    }
}

GPS数据解析

//获取当前数据中指定逗号后面的字符在该数据中的位置，用于定位具体数据块的位置
uint8_t GetComma(char* str,uint8_t num)
{
    uint8_t cnt=0,tmp=0,comma=0;
    while(comma!=num)
    {
        tmp=str[cnt++];
        if(tmp=='\0') return 0;
        if(tmp==',') comma++;
    }
    return cnt;
}

//获取指定数据块的数据并转换成浮点数
double GetDoubleNumber(char *s)
{
    char buf[128];
    double rev;
    int i = GetComma(s,1);
    strncpy(buf,s,i);
    buf[i] = '\0';
    rev = atof(buf);
    return rev;
}

//将度分秒格式转换为度数
double DMM2Degree(double dmm)
{
    return (int)dmm + (dmm - (int)dmm)/60.0*100.0;
}

void UTC2BTC(GPSMsgType* gps)
{
    #define isLeap(year) (((year%4==0&&year/100!=0)||(year%400==0)))
    static int day_list[2][12]={{31,28,31,30,31,30,31,31,30,31,30,31},{31,29,31,30,31,30,31,31,30,31,30,31}};
    gps->hour += 8;
    if(gps->hour > 23)
    {
        gps->hour -= 24;
        gps->day++;
        if(gps->day > day_list[isLeap(gps->year)][gps->month-1])
        {
            gps->day -= day_list[isLeap(gps->year)][gps->month-1];
            gps->month++;
            if(gps->month > 12)
            {
                gps->month -= 12;
                gps->year++;
            }
        }
    }
}

typedef data struct
{
    double  latitude; //经度
    double  longitude; //纬度
    float   speed; //速度
    float   direction; //航向
    float   height; //海拔高度
    int     satellite;
    char    NS;
    char    EW;
    int     year; 
    int     month; 
    int     day;
    int     hour;
    int     minute;
    int     second;
}GPS_INFO;

GPS数据转换

地图坐标

由于与GPS相关的坐标系有若干种，在开发与地点相关的应用时，如果对坐标系不加处理，很容易出现几种形式的地图之间出现坐标点偏移的情况。因此在坐标系之间进行转换非常重要。

• WGS-84：世界大地测量系统，国际通用标准
  • GPS模块、国际版Google Map使用
• GCJ-02：中国坐标偏移标准，俗称火星坐标系
  • 遵循国内测绘法，基于WGS-84系统加入随机漂移
  • 国内版谷歌地图、高德地图、腾讯地图使用
• BD-09：百度坐标偏移标准
  • 基于GCJ-02多增加一次变换
  • 百度地图使用

C语言代码

基于网上开源内容整理了C语言实现GPS坐标转换的代码。

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

#define PI      3.14159265358979324
#define AXIS    6378245.0                   //卫星椭球坐标投影到平面地图坐标系的投影因子
#define OFFSET  0.00669342162296594323      //椭球的偏心率
#define X_PI    PI*3000.0/180.0

typedef struct
{
    double lat;
    double lon;
}LLType;

//Area out of China don't need axis transfer
int outOfChina(LLType tmpll)
{
    if (tmpll.lon < 72.004 || tmpll.lon > 137.8347)  return 1;
    if (tmpll.lat < 0.8293 || tmpll.lat > 55.8271)   return 1;
    return 0;
}

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(fabs(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(fabs(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;  
}

LLType delta(LLType tmpll)
{
    LLType latlng;
    double dLat = transformLat(tmpll.lon - 105.0, tmpll.lat - 35.0);  
    double dLon = transformLon(tmpll.lon - 105.0, tmpll.lat - 35.0);  
    double radLat = tmpll.lat / 180.0 * PI;  
    double magic = sin(radLat);
    magic = 1 - OFFSET * magic * magic;  
    double sqrtMagic = sqrt(magic);  
    latlng.lat = (dLat * 180.0) / ((AXIS * (1 - OFFSET)) / (magic * sqrtMagic) * PI);
    latlng.lon = (dLon * 180.0) / (AXIS / sqrtMagic * cos(radLat) * PI);
    return latlng;
}

LLType WGS2GCJ(LLType wgsll)
{
    if(outOfChina(wgsll))   return wgsll;

    LLType gcjll, deltaD = delta(wgsll);
    gcjll.lat = wgsll.lat + deltaD.lat;
    gcjll.lon = wgsll.lon + deltaD.lon;
    return gcjll;
}

LLType GCJ2WGS(LLType gcjll)
{
    if(outOfChina(gcjll))   return gcjll;

    LLType wgsll, deltaD = delta(gcjll);
    wgsll.lat = gcjll.lat - deltaD.lat;
    wgsll.lon = gcjll.lon - deltaD.lon;
    return wgsll;
}

//使用二分极限法
LLType GCJ2WGSExactly(LLType gcjll)
{
    double initDelta = 0.01;
    double threshold = 0.000000001;
    double dLat = initDelta, dLon = initDelta;
    double mLat = gcjll.lat - dLat, mLon = gcjll.lon - dLon;
    double pLat = gcjll.lat + dLat, pLon = gcjll.lon + dLon;
    double i = 0;
    LLType wgsll;
    while (1) {
        wgsll.lat = (mLat + pLat) / 2;
        wgsll.lon = (mLon + pLon) / 2;
        LLType tmp = WGS2GCJ(wgsll);
        dLat = tmp.lat - gcjll.lat;
        dLon = tmp.lon - gcjll.lon;
        if ((fabs(dLat) < threshold) && (fabs(dLon) < threshold))   break;

        if (dLat > 0) pLat = wgsll.lat; else mLat = wgsll.lat;
        if (dLon > 0) pLon = wgsll.lon; else mLon = wgsll.lon;

        if (++i > 10000) break;
    }
    return wgsll;
}

LLType GCJ2BD(LLType gcjll)
{
    LLType bdll;
	double x = gcjll.lon, y = gcjll.lat;  
	double z = sqrt(x * x + y * y) + 0.00002 * sin(y * X_PI);  
	double theta = atan2(y, x) + 0.000003 * cos(x * X_PI);  
	bdll.lat = z * sin(theta) + 0.006;
	bdll.lon = z * cos(theta) + 0.0065;  
	return bdll;
}

LLType BD2GCJ(LLType bdll)
{
    LLType gcjll;
    double x = bdll.lon - 0.0065;
    double y = bdll.lat - 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);  
    gcjll.lat = z * sin(theta);
    gcjll.lon = z * cos(theta);  
    return gcjll;
}

LLType WGS2BD(LLType wgsll)
{
    return GCJ2BD(WGS2GCJ(wgsll));
}

LLType BD2WGS(LLType bdll)
{
    return GCJ2WGS(BD2GCJ(bdll));
}

double Distance(LLType all, LLType bll)
{
    int earthR = 6371000;
    double x = cos(all.lat * PI/180) * cos(bll.lat * PI/180) * cos((all.lon - bll.lon) * PI/180);
    double y = sin(all.lat * PI/180) * sin(bll.lat * PI/180);
    double s = x + y;
    if (s > 1)  s = 1;
    if (s < -1) s = -1;
    double alpha = acos(s);
    double distance = alpha * earthR;
    return distance;
}

int main()
{
    LLType first,middle,last;
    first.lat = 22.502412986242;
    first.lon = 113.93832783228;

    // WGS->GCJ->BD
    middle = WGS2GCJ(first);
    last = GCJ2BD(middle);
    printf("WGS84: %f,%f\r\n",first.lat,first.lon);
    printf("GCJ02: %f,%f\r\n",middle.lat,middle.lon);
    printf("BD09:  %f,%f\r\n",last.lat,last.lon);

    //BD->GCJ->WGS
    middle = BD2GCJ(last);
    first = GCJ2WGS(middle); // first = GCJ2WGSExactly(middle);
    printf("GCJ02: %f,%f\r\n",middle.lat,middle.lon);
    printf("WGS84: %f,%f\r\n",first.lat,first.lon);

    return 0;
}

可以使用这个在线工具进行测试：C 在线工具 | 菜鸟工具

参考文献

GPS坐标互转：WGS-84(GPS)、GCJ-02(Google地图)、BD-09(百度地图) | JS,PHP
gps各种地图坐标系转换 | Java