[python] wgs84转为gcj02坐标

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import math
import pandas as pd

EE = 0.00669342162296594323
EARTH_RADIUS = 6378.137 #地球近似半径(km)
M_PI = math.pi

def transformLat( x,  y) :
    ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(math.fabs(x));
    ret += (20.0 * math.sin(6.0 * x * M_PI) + 20.0 * math.sin(2.0 * x * M_PI)) * 2.0 / 3.0;
    ret += (20.0 * math.sin(y * M_PI) + 40.0 * math.sin(y / 3.0 * M_PI)) * 2.0 / 3.0;
    ret += (160.0 * math.sin(y / 12.0 * M_PI) + 320 * math.sin(y * M_PI / 30.0)) * 2.0 / 3.0;
    return ret;

def transformLon( x,  y) :
    ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(math.fabs(x));
    ret += (20.0 * math.sin(6.0 * x * M_PI) + 20.0 * math.sin(2.0 * x * M_PI)) * 2.0 / 3.0;
    ret += (20.0 * math.sin(x * M_PI) + 40.0 * math.sin(x / 3.0 * M_PI)) * 2.0 / 3.0;
    ret += (150.0 * math.sin(x / 12.0 * M_PI) + 300.0 * math.sin(x / 30.0 * M_PI)) * 2.0 / 3.0;
    return ret;
 
def Wgs84ToGcj02(lon, lat) :
    dLat = transformLat(lon - 105.0, lat - 35.0);
    dLon = transformLon(lon - 105.0, lat - 35.0);
    radLat = lat / 180.0 * M_PI;
    magic = math.sin(radLat);
    magic = 1 - magic * magic * EE;
    sqrtMagic = math.sqrt(magic);

    dLat = (dLat * 180.0) / ((EARTH_RADIUS * 1000 * (1 - EE)) / (magic * sqrtMagic) * M_PI);
    dLon = (dLon * 180.0) / (EARTH_RADIUS * 1000 / sqrtMagic * math.cos(radLat) * M_PI);
    y02 = lat + dLat;
    x02 = lon + dLon;
 
    return x02,y02

def convert_csv(in_csv, use_cols, xy_cols, to_gcj02, time1, time2, check_time, has_header, out_csv):
    if len(xy_cols)<2:
        print('xy_cols size error', xy_cols)
        return
    if len(use_cols)<3:
        print('use_cols size error, first col must timesatmp', xy_cols)
        return

    x_col_idx = use_cols.index(xy_cols[0])
    y_col_idx = use_cols.index(xy_cols[1])

    oth_data = None
    
    if has_header :
        oth_data = pd.read_csv(in_csv, index_col=None, usecols=use_cols, encoding='utf8')
    else:
        oth_data = pd.read_csv(in_csv, index_col=None, usecols=use_cols, header=None, encoding='utf8')

    cols = []
    for a in oth_data:
        cols.append(a)

    if not has_header :    
        cols[x_col_idx]='x'
        cols[y_col_idx]='y'

    ds = []
    for oth in oth_data.itertuples():
        ks = list(oth[1:])
        # print('1', ks)
        
        if check_time : 
            if ks[0]>time2 or ks[0]<time1:
                continue

        if to_gcj02:
            x = ks[x_col_idx]
            y = ks[y_col_idx]
            x02,y02 = Wgs84ToGcj02(x,y)
            ks[x_col_idx] = x02
            ks[y_col_idx] = y02
        # print('2', ks)

        ds.append(ks)

    df = pd.DataFrame(ds, columns=cols)
    print(df)
    df.to_csv(out_csv, index = False)


if __name__ == "__main__":
    
    in_file = r'/home/w/temp/20200608/G6/M8P/record/G9DATA.csv'
    out_file = r'/home/w/temp/20200608/G6/M8P/record/G9DATA_022.csv'
    use_cols = [0,5,6,7]
    xy_cols = [5,6]
    to_gcj02 = True
    time1 = 1
    time2 = 2
    check_time = False
    has_header = True
    # convert_csv(in_file, use_cols, xy_cols, to_gcj02, time1, time2, check_time, has_header, out_file)

    in_file = r'/home/w/temp/20200608/G6/M8P/rtk.txt'
    out_file = r'/home/w/temp/20200608/G6/M8P/rtk_02_slim1.txt'
    use_cols = [0,1,2,3,4,5,6,7]
    xy_cols = [1,2]
    to_gcj02 = True
    time1 = 1591564304940
    time2 = 1591566905740
    check_time = True
    has_header = False
    convert_csv(in_file, use_cols, xy_cols, to_gcj02, time1, time2, check_time, has_header, out_file)

    in_file = r'/home/w/temp/20200608/G6/M8P/rtk_02.txt'
    out_file = r'/home/w/temp/20200608/G6/M8P/rtk_02_slim.txt'
    use_cols = [0,1,2,3,4,5,6,7]
    xy_cols = [1,2]
    to_gcj02 = False
    time1 = 1591564304940
    time2 = 1591566905740
    check_time = True
    has_header = False
    # convert_csv(in_file, use_cols, xy_cols, to_gcj02, time1, time2, check_time, has_header, out_file)