分类 Python 开发 下的文章


1、GPS时间是由:周、秒组成

     

GPS时间系统是基于GPS星载原子钟和地面监控站原子钟的一种原子时基准,‌与国际原子时保持有19秒的常数差,‌并在GPS标准历元1980年1月6日零时与UTC保持一致。‌

GPS时间系统(‌GPST)‌的起点为1980年1月6日0时0分0秒,‌在起始时刻,‌GPST与UTC时对齐,‌两种时间系统给出的时间相同。‌由于UTC时存在跳秒现象,‌因此一段时间后,‌两种时间系统会相差n个整秒,‌其中n为这段时间内UTC积累的跳秒数。‌GPS时间系统的最大时间单位是周,‌一周等于604800秒。‌由于在储存周数的时候,‌用了10个比特,‌2的10次方是1024,‌所以每1024周(‌即7168天)‌为一循环周期。‌当WN从1023变为0时,‌就会发生GPS周数翻转,‌出现迎接新一周的说法。‌1024周对应到年上大概就是19.7年。‌从GPS系统时的起始时刻算起,‌上一次出现GPS周数翻转是1999年8月21日,‌这次就正好是2019年4月6日,‌2038年11月20日将会出现下一次GPS周数翻转。‌

这种周数翻转的现象类似于小朋友数数时每次数到100就又从0开始数的情形,‌而按10进制的计数规则,‌100以后是101,‌200以后是201…‌…‌以此类推。‌对于GPS接收机来说,‌如果在周数翻转时刻未升级,‌都将在周归零日错认为1980年1月6日。‌为了避免影响开发者使用,‌高德开放平台在GPS周数翻转时刻采用了系统时间和GPS对比的方法,‌以解决周数翻转带来的时间跳变问题。‌



import  sys
from datetime  import datetime, timedelta






def     gps_to_utc(gps_week, gps_seconds):
    # GPS 周起始时间
    gps_epoch = datetime(1980, 1, 6)
    # 计算GPS周对应的总天数
    total_days = gps_week * 7
    #计算 GPS秒对应的时间差,GPS时间与UTC时间在秒上相差18秒
    time_diff = timedelta(days=total_days, seconds=gps_seconds - 18)

    #计算UTC时间
    utc_time = gps_epoch + time_diff

    return  utc_time



def     GPS_To_UTC(argv):
    gps_weeks = int(argv[0])
    gps_secs  = float(argv[1])

    utc_time = gps_to_utc(gps_weeks, gps_secs)
    print(utc_time)




if __name__ == "__main__":
    if len(sys.argv)>3:
        print("Usage:\n  < GPS weeks >  < GPS seconds >")
        sys.exit()
    if 3== len(sys.argv):
        GPS_To_UTC(sys.argv[1:])
    else:
        gps_weeks = 2322
        gps_secs  = 112576.3 #output: 2024.07.08 07:15:58.3
        utc_time = gps_to_utc(gps_weeks, gps_secs)
        print(utc_time)




















Project description

China Coordinate Convertor

Actions StatusGitHub starsGitHub forksGitHub license

中国火星坐标转换命令行工具,用于WGS-84(未偏移坐标), GCJ-02(国家测绘局、高德、谷歌中国地图), BD-09(百度坐标系)三者之间的互相转换,支持文件格式:

  • ESRI Shapefile

  • GeoJSON

    图片名称

安装

推荐使用:pip install coord-convert

或者也可以从源码安装:

git clone https://github.com/sshuair/coord-convert.gitpip install -r requirements.txtpython setup.py install

依赖

  • python3
  • fiona
  • tqdm
  • click

使用方法

注意:火星坐标转换是针对经纬度的转换,因此在进行转换前需要将坐标转换成经纬度;比如web墨卡托(3857)等投影坐标系需要先转成经纬度坐标(4326)

Python API调用

在python程序中调用相应的坐标转换接口

from coord_convert.transform import wgs2gcj, wgs2bd, gcj2wgs, gcj2bd, bd2wgs, bd2gcj lon, lat = 120, 40gcj_lon, gcj_lat = wgs2gcj(lon, lat)bd_lon, bd_lat = wgs2bd(lon, lat)print(gcj_lon, gcj_lat) # the result should be: 120.00567568355486 40.0013047896019

命令行调用

火星坐标转换还支持命令行直接对shp、geojson等文件进行转换,比如

~/temp > coord_covert gcj2wgs tests/data/Polyline/polyline.geojson aa.geojson100%|██████████████████████████████████████████████████████| 219/219 [00:00<00:00, 550.93it/s]

更详细的用法:

convert input china coordinate to another.     Arguments:        convert_type {string} -- [coordinate convert type, e.g. wgs2bd]            wgs2gcj : convert WGS-84 to GCJ-02            wgs2bd  : convert WGS-84 to DB-09            gcj2wgs : convert GCJ-02 to WGS-84            gcj2bd  : convert GCJ-02 to BD-09            bd2wgs  : convert BD-09 to WGS-84            bd2gcj  : convert BD-09 to GCJ-02        src_path {string} -- [source file path]        dst_path {string} -- [destination file path]    Example:        coord_covert wgs2gcj ./tests/data/polygon/polygon_wgs.shp ~/temp/polygon_gcj.shp 



一、安装依赖

     1、folium, 地图显示

        pip3   install     folium

     2、coord-convert,坐标转换

      pip3 install coord-convert

   或源码安装

   git clone https://github.com/sshuair/coord-convert.git

   pip install -r requirements.txt

   python setup.py install


       

二、代码实现

     import os
import sys
import pandas as pd
import numpy as np

#import  matplotlib.pyplot as plt

import  folium
from pyproj import Transformer
from coord_convert.transform import wgs2gcj, wgs2bd, gcj2wgs, gcj2bd, bd2wgs, bd2gcj




def     Load_GNSS_Dataset(csv_file):
    data = pd.read_csv(csv_file)
    #print(data, type(data))
    frame = pd.DataFrame(data)
    t = frame.iloc[:, 0]
    lat = frame.iloc[:, 1]
    lon = frame.iloc[:, 2]
    alt = frame.iloc[:, 3]

    return  t.to_numpy(), lat.to_numpy(), lon.to_numpy(), alt.to_numpy()
    
    
    
def     Show_GPS_Trajectory_Map(argv):
    gnss_file = argv[0]
    t, wgs84_lats, wgs84_lons, wgs84_alts = Load_GNSS_Dataset(gnss_file)
    #lons, lats = wgs2gcj(wgs84_lons[0], wgs84_lats[0])
    Lats = []
    Lons = []
    for  i in range(len(wgs84_lats)):
        gcj_lon, gcj_lat = wgs2gcj(wgs84_lons[i], wgs84_lats[i])
        Lats.append(gcj_lat)
        Lons.append(gcj_lon)
    lats = np.array(Lats)
    lons = np.array(Lons)

    # create map
    Map = folium.Map(
            location=[lats[0], lons[0]], #map center
            zoom_start=12,
            tiles= 'https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scl=1&style=7',
            attr='高德-常规图',
            
            #tiles='http://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineCommunity/MapServer/tile/{z}/{y}/{x}',
            #attr='彩色版'

            #tiles =  'https://rt0.map.gtimg.com/tile?z={z}&x={x}&y={-y}',
            #attr='腾讯地图'
            )
    folium.PolyLine(locations=[list(zip(lats, lons))], color='blue', weight=5).add_to(Map)
    Map.save('GPS_track.html');



if __name__ == "__main__":
    if len(sys.argv)<2:
        print("Usage:\n  < GPS csv>")
        sys.exit()
    #print(sys.argv[1:], type(sys.argv[1:]))
    if 2== len(sys.argv):
        sys.exit(Show_GPS_Trajectory_Map(sys.argv[1:]))
    else:
        sys.exit()






















def     Show_Position_Trajectory(argv):
    gnss = argv[0]

    t, x, y, z, sigma_x, sigma_y, sigma_z = Load_GNSS_Dataset(gnss)


    #fig, (ax1, ax2, ax3, ax4) = plt.subplots(4)
    fig, axs = plt.subplots(2, 2)
    axs[0,0].plot(x, y, '-', color='red', label='GPS')
    axs[0,0].set_title('WGS84 latitude,longitude')
    axs[0,0].set(xlabel='E(degree)', ylabel='N(degree)')

    axs[0,1].plot(t, sigma_x, '-', color='red', label='x axis sigma')
    axs[0,1].set_title('GNSS x sigma[ENU]')
    axs[0,1].set(xlabel='t(s)', ylabel='m')

    axs[1,0].plot(t, sigma_y, '-', color='green', label='y axis sigma')
    axs[1,0].set_title('GNSS y sigma[ENU]')
    axs[1,0].set(xlabel='t(s)', ylabel='m')

    axs[1,1].plot(t, sigma_z, '-', color='blue', label='z axis sigma')
    axs[1,1].set_title('GNSS z sigma[ENU]')
    axs[1,1].set(xlabel='t(s)', ylabel='m')


    plt.legend()
    plt.show()
    return






def     Show_IMU_Data(argv):
    imu = argv[0]

    t, gyro_x, gyro_y, gyro_z, acce_x, acce_y, acce_z = Load_IMU_Dataset(imu)


    # gyro
    fig, axs = plt.subplots(3, 2)
    axs[0,0].plot(t, gyro_x, '-', color='red', label='Gyro-x')
    axs[0,0].set_title('Gyro')
    axs[0,0].set(xlabel='t(s)', ylabel='ang(rad/s)')
    axs[0,0].grid()

    axs[0,0].legend()


    axs[1,0].plot(t, gyro_y, '-', color='green', label='Gyro-y')
    axs[1,0].set_title('Gyro')
    axs[1,0].set(xlabel='t(s)', ylabel='ang(rad/s)')

    axs[1,0].grid()

    axs[1,0].legend()


    axs[2,0].plot(t, gyro_z, '-', color='blue', label='Gyro-z')
    axs[2,0].set_title('Gyro')
    axs[2,0].set(xlabel='t(s)', ylabel='ang(rad/s)')
  axs[2,0].grid()

    axs[2,0].legend()



    # acce
    axs[0,1].plot(t, acce_x, '-', color='red', label='Acce-x')
    axs[0,1].set_title('Acce')
    axs[0,1].set(xlabel='t(s)', ylabel='acc(m/s^2)')

    axs[0,1].grid()

    axs[0,1].legend()


    axs[1,1].plot(t, acce_y, '-', color='blue', label='Acce-y')
    axs[1,1].set_title('Acce')
    axs[1,1].set(xlabel='t(s)', ylabel='acc(m/s^2)')

    axs[1,1].grid()

    axs[1,1].legend()


    axs[2,1].plot(t, acce_z, '-', color='blue', label='Acce-z')
    axs[2,1].set_title('Acce')
    axs[2,1].set(xlabel='t(s)', ylabel='acc(m/s^2)')
    axs[2,1].grid()

    axs[2,1].legend()



    plt.legend()
    plt.show()
    return






def     Compare_GNSS_Trajectory(argv):
    gnss_gt  = argv[0]
    gnss_ekf = argv[1]

    tt, xt, yt, zt = Load_GNSS_Dataset(gnss_gt)
    t, X, Y, Z = Load_GNSS_Dataset(gnss_ekf)
    plt.plot(yt, xt, '-', color='blue', label='Ground Truth')

    #plt.scatter(Y, X, color='dodgerblue', label='ESKF')
    plt.plot(Y, X, '-', color='red', label='EKF/ESKF')

    fig, axs = plt.subplots(2, 1)
    axs[0].plot(tt, xt, '-', color='blue', label='Ground Truth')
    axs[0].plot(t, X, '-', color='red', label='EKF/ESKF')
    axs[0].set_title('Position[East]')
    axs[0].set(xlabel='t(s)', ylabel='East(m)')


    axs[1].plot(tt, yt, '-', color='blue', label='Ground Truth')
    axs[1].plot(t, Y, '-', color='red', label='EKF/ESKF')
    axs[1].set_title('Position[North]')
    axs[1].set(xlabel='t(s)', ylabel='North(m)')

    plt.legend()
    plt.show()
    return