经纬度转换(go/python/rust)

golang

代码

github.com/yanmengfei/coord

go get github.com/yanmengfei/coord@v0.0.1

使用

package main

import (
    "fmt"
    "github.com/yanmengfei/coord"
)

func main() {
    var location = "115.668055,34.449162"
    var lon, lat, _ = coord.LocationToFloat64Coord(location)
    fmt.Println("GCJ02:", lon, lat)
    lon, lat = coord.GCJ02toWGS84(lon, lat)
    fmt.Println("WGS84:", lon, lat)
}

rust

pub mod coords_convert {
    use std::f64::consts::PI;
    // Pi*3000/180
    const XPI: f64 = 52.35987755982988;
    // 长半轴
    const AXIS: f64 = 6378245.0;
    // 偏心率平方
    const OFFSET: f64 = 0.00669342162296594323;
    // 百度坐标系经度偏移
    const BD_LON_OFFSET: f64 = 0.0065;
    // 百度坐标系纬度偏移
    const BD_LAT_OFFSET: f64 = 0.0060;
    const BD_Z_OFFSET: f64 = 0.00002;
    const BD_THETA_OFFSET: f64 = 0.000003;
    const ACCURACY_THRESHOLD: f64 = 0.0000000001;
    const WGS_OFFSET: f64 = 0.01;

    fn in_china(lon: f64, lat: f64) -> bool {
        (73.66 < lon && lon < 135.05) && (3.86 < lat && lat < 53.55)
    }

    fn convert(lon: f64, lat: f64) -> (f64, f64) {
        let lonlat = lon * lat;
        let abs_x = lon.abs().sqrt();
        let lon_pi = lon * PI;
        let lat_pi = lat * PI;
        let d = (lon_pi * 6.0).sin() * 20.0 + (lon_pi * 2.0).sin() * 20.0;
        let mut x = d;
        let mut y = d;
        x += 20.0 * lat_pi.sin() + (lat_pi / 3.0).sin() * 40.0;
        x += 160.0 * (lat_pi / 12.0).sin() + 320.0 * (lat_pi / 30.0).sin();

        y += 20.0 * lon_pi.sin() + (lon_pi / 3.0).sin() * 40.0;
        y += 150.0 * (lon_pi / 12.0).sin() + 300.0 * (lon_pi / 30.0).sin();
        let threshold: f64 = 2.0 / 3.0;
        x *= threshold;
        y *= threshold;
        x += 2.0 * lon + 3.0 * lat + 0.2 * lat * lat + 0.1 * lonlat + 0.2 * abs_x - 100.0;
        y += lon + 2.0 * lat + 0.1 * lon * lon + 0.1 * lonlat + 0.1 * abs_x + 300.0;
        return (x, y);
    }

    fn delta(lon: f64, lat: f64) -> (f64, f64) {
        let (de_lat, de_lon) = convert(lon - 105.0, lat - 35.0);
        let rad_lat = lat / 180.0 * PI;
        let magic = rad_lat.sin();
        let magic = 1.0 - OFFSET * magic * magic;
        let sqrtmagic = magic.sqrt();
        let de_lat = (de_lat * 180.0) / ((AXIS * (1.0 - OFFSET)) / (magic * sqrtmagic) * PI);
        let de_lon = (de_lon * 180.0) / (AXIS / sqrtmagic * rad_lat.cos() * PI);
        return (lon + de_lon, lat + de_lat);
    }

    pub fn wgs84_to_gcj02(lon: f64, lat: f64) -> (f64, f64) {
        if !in_china(lon, lat) {
            return (lon, lat);
        }
        return delta(lon, lat);
    }

    pub fn gcj02_to_bd09(lon: f64, lat: f64) -> (f64, f64) {
        let z = (lon * lon + lat * lat).sqrt() + BD_Z_OFFSET * (lat * XPI).sin();
        let theta = lat.atan2(lon) + BD_THETA_OFFSET * (lon * XPI).cos();
        return (
            z * theta.cos() + BD_LON_OFFSET,
            z * theta.sin() + BD_LAT_OFFSET,
        );
    }

    pub fn bd09_to_gcj02(lon: f64, lat: f64) -> (f64, f64) {
        let x = lon - BD_LON_OFFSET;
        let y = lat - BD_LAT_OFFSET;
        let z = (x * x + y * y).sqrt() - BD_Z_OFFSET * (y * XPI).sin();
        let theta = y.atan2(x) - BD_THETA_OFFSET * (x * XPI).cos();
        return (z * theta.cos(), z * theta.sin());
    }

    pub fn gcj02_to_wgs84(lon: f64, lat: f64) -> (f64, f64) {
        let mut mlon = lon - WGS_OFFSET;
        let mut mlat = lat - WGS_OFFSET;
        let mut plon = lon + WGS_OFFSET;
        let mut plat = lat + WGS_OFFSET;
        let mut dlon: f64;
        let mut dlat: f64;
        let mut wgs_lon: f64 = 0.0;
        let mut wgs_lat: f64 = 0.0;

        for _ in 0..10000 {
            wgs_lat = (mlat + plat) / 2.0;
            wgs_lon = (mlon + plon) / 2.0;
            let (tmp_lon, tmp_lat) = delta(wgs_lon, wgs_lat);
            dlon = tmp_lon - lon;
            dlat = tmp_lat - lat;
            if (dlat.abs() < ACCURACY_THRESHOLD) && (dlon.abs() < ACCURACY_THRESHOLD) {
                break;
            }
            if dlat > 0.0 {
                plat = wgs_lat;
            } else {
                mlat = wgs_lat;
            }
            if dlon > 0.0 {
                plon = wgs_lon;
            } else {
                mlon = wgs_lon;
            }
        }
        return (wgs_lon, wgs_lat);
    }
}

#[cfg(test)]
mod point_tests {
    use super::coords_convert;

    #[test]
    fn wgs84_to_gcj02_test() {
        assert_eq!(
            coords_convert::wgs84_to_gcj02(116.65383912049425, 39.992875254426366),
            (116.6597270000643, 39.99397599992571)
        )
    }

    #[test]
    fn gcj02_to_bd09_test() {
        assert_eq!(
            coords_convert::gcj02_to_bd09(113.56733, 34.81754),
            (113.5739274927449, 34.82327621798387)
        )
    }

    #[test]
    fn bd09_to_gcj02_test() {
        assert_eq!(
            coords_convert::bd09_to_gcj02(113.5739274927449, 34.82327621798387),
            (113.56732983450783, 34.81754111708383)
        )
    }

    #[test]
    fn gcj02_to_wgs84_test() {
        assert_eq!(
            coords_convert::gcj02_to_wgs84(113.5739274927449, 34.82327621798387),
            (113.56784261660992, 34.82437523040346)
        )
    }
}

python

#!/usr/bin/python
# coding:utf-8
# Filename:ConvertCoordinate.py
# Used To Convert Coordinate System.


import math


def main():
    # bd09转gcj02
    bd_lat, bd_lon = 112, 33
    gcj_lat, gcj_lon = bd09_to_gcj02(bd_lat, bd_lon)

    # bd09转wgs84
    bd_lat, bd_lon = 112, 33
    gcj_lat, gcj_lon = bd09_to_gcj02(bd_lat, bd_lon)
    wgs_lat, wgs_lon = gcj02_to_wgs84(gcj_lat, gcj_lon)

    # wgs84转gcj02
    wgs_lat, wgs_lon = 112, 33
    gcj_lat, gcj_lon = wgs84_to_gcj02(wgs_lat, wgs_lon)

    # wgs84转bd09
    wgs_lat, wgs_lon = 112, 33
    gcj_lat, gcj_lon = wgs84_to_gcj02(wgs_lat, wgs_lon)
    bd_lat, bd_lon = gcj02_to_bd09(gcj_lat, gcj_lon)

    # gcj02转bd09
    gcj_lat, gcj_lon = 112, 33
    bd_lat, bd_lon = gcj02_to_bd09(gcj_lat, gcj_lon)

    # gcj02转wgs84
    gcj_lat, gcj_lon = 112, 33
    wgs_lat, wgs_lon = gcj02_to_wgs84(gcj_lat, gcj_lon)


# BD09经纬度到GCJ02经纬度转换
def bd09_to_gcj02(bd_lat, bd_lon):
    x_pi = float(3.14159265358979324 * 3000.0 / 180.0)
    x, y = float(float(bd_lon) - 0.0065), float(float(bd_lat) - 0.006)
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)
    gcj_lon = z * math.cos(theta)
    gcj_lat = z * math.sin(theta)
    gcj_lon = format(float(gcj_lon), '.8f')
    gcj_lat = format(float(gcj_lat), '.8f')
    return gcj_lat, gcj_lon


# GCJ02经纬度到BD09经纬度转换
def gcj02_to_bd09(gcj_lat, gcj_lon):
    x_pi = 3.14159265358979324 * 3000.0 / 180.0
    gcj_lon, gcj_lat = float(gcj_lon), float(gcj_lat)
    z = math.sqrt(gcj_lon * gcj_lon + gcj_lat * gcj_lat) + 0.00002 * math.sin(gcj_lat * x_pi)
    theta = math.atan2(gcj_lat, gcj_lon) + 0.000003 * math.cos(gcj_lon * x_pi)
    bd_lon = z * math.cos(theta) + 0.0065
    bd_lat = z * math.sin(theta) + 0.006
    bd_lon = format(float(bd_lon), '.6f')
    bd_lat = format(float(bd_lat), '.6f')
    return bd_lat, bd_lon


# WGS84经纬度到GCJ02经纬度转换
def wgs84_to_gcj02(wgs_lat, wgs_lon):
    wgs_lat, wgs_lon = float(wgs_lat), float(wgs_lon)
    pi = 3.1415926535897932384626  # π
    a = 6378245.0  # 长半轴
    ee = 0.00669342162296594323  # 扁率

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

    def transformlon(x, y):
        result = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(math.fabs(x))
        result += (20.0 * math.sin(6.0 * x * pi) + 20.0 * math.sin(2.0 * x * pi)) * 2.0 / 3.0
        result += (20.0 * math.sin(x * pi) + 40.0 * math.sin(x / 3.0 * pi)) * 2.0 / 3.0
        result += (150.0 * math.sin(x / 12.0 * pi) + 300.0 * math.sin(x / 30.0 * pi)) * 2.0 / 3.0
        return result

    dlat = transformlat(wgs_lon - 105.0, wgs_lat - 35.0)
    dlon = transformlon(wgs_lon - 105.0, wgs_lat - 35.0)
    radlat = float(wgs_lat / 180.0 * pi)
    magic = float(math.sin(radlat))
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlon = (dlon * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    gcj_lat = format(float(wgs_lat + dlat), '.8f')
    gcj_lon = format(float(wgs_lon + dlon), '.8f')
    return gcj_lat, gcj_lon


# GCJ02经纬度到WGS84经纬度转换
def gcj02_to_wgs84(gcj_lat, gcj_lon):
    gcj_lat, gcj_lon = float(gcj_lat), float(gcj_lon)
    initdelta = 0.01
    threshold = 0.0000000001
    dlat, dlon = initdelta, initdelta
    mlat = gcj_lat - dlat
    mlon = gcj_lon - dlon
    plat = gcj_lat + dlat
    plon = gcj_lon + dlon
    wgs_lat, wgs_lon, i = 0, 0, 0
    while True:
        wgs_lat = (mlat + plat) / 2
        wgs_lon = (mlon + plon) / 2
        tmplat, tmplon = wgs_corrected(wgs_lat, wgs_lon)
        tmplat, tmplon = float(tmplat), float(tmplon)
        dlat = tmplat - gcj_lat
        dlon = tmplon - gcj_lon
        if (math.fabs(dlat) < threshold) and (math.fabs(dlon) < threshold):
            break

        if dlat > 0:
            plat = wgs_lat
        else:
            mlat = wgs_lat
        if dlon > 0:
            plon = wgs_lon
        else:
            mlon = wgs_lon
        i += 1
        if i > 10000:
            break
    wgs_lat, wgs_lon = format(float(wgs_lat), '.6f'), format(float(wgs_lon), '.6f')
    return wgs_lat, wgs_lon


# GPS纠偏算法，提高精确度
def wgs_corrected(wgs_lat, wgs_lon):
    """GPS纠偏算法，提高精确度"""
    pi = 3.1415926535897932384626  # π
    a = 6378245.0  # 长半轴
    ee = 0.00669342162296594323  # 扁率

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

    def transformlon(x, y):
        result = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(math.fabs(x))
        result += (20.0 * math.sin(6.0 * x * pi) + 20.0 * math.sin(2.0 * x * pi)) * 2.0 / 3.0
        result += (20.0 * math.sin(x * pi) + 40.0 * math.sin(x / 3.0 * pi)) * 2.0 / 3.0
        result += (150.0 * math.sin(x / 12.0 * pi) + 300.0 * math.sin(x / 30.0 * pi)) * 2.0 / 3.0
        return result

    dlat = transformlat(wgs_lon - 105.0, wgs_lat - 35.0)
    dlon = transformlon(wgs_lon - 105.0, wgs_lat - 35.0)
    radlat = float(wgs_lat / 180.0 * pi)
    magic = float(math.sin(radlat))
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlon = (dlon * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    gcj_lat = format(float(wgs_lat + dlat), '.8f')
    gcj_lon = format(float(wgs_lon + dlon), '.8f')

    return gcj_lat, gcj_lon


if __name__ == '__main__':
    main()