经纬度转换(go/python/rust)
生活随笔
收集整理的這篇文章主要介紹了
经纬度转换(go/python/rust)
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
golang
代碼
github.com/yanmengfei/coord
go get github.com/yanmengfei/coord@v0.0.1使用
package mainimport ("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/180const 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 mathdef main():# bd09轉gcj02bd_lat, bd_lon = 112, 33gcj_lat, gcj_lon = bd09_to_gcj02(bd_lat, bd_lon)# bd09轉wgs84bd_lat, bd_lon = 112, 33gcj_lat, gcj_lon = bd09_to_gcj02(bd_lat, bd_lon)wgs_lat, wgs_lon = gcj02_to_wgs84(gcj_lat, gcj_lon)# wgs84轉gcj02wgs_lat, wgs_lon = 112, 33gcj_lat, gcj_lon = wgs84_to_gcj02(wgs_lat, wgs_lon)# wgs84轉bd09wgs_lat, wgs_lon = 112, 33gcj_lat, gcj_lon = wgs84_to_gcj02(wgs_lat, wgs_lon)bd_lat, bd_lon = gcj02_to_bd09(gcj_lat, gcj_lon)# gcj02轉bd09gcj_lat, gcj_lon = 112, 33bd_lat, bd_lon = gcj02_to_bd09(gcj_lat, gcj_lon)# gcj02轉wgs84gcj_lat, gcj_lon = 112, 33wgs_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.0gcj_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.0065bd_lat = z * math.sin(theta) + 0.006bd_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.0result += (20.0 * math.sin(y * pi) + 40.0 * math.sin(y / 3.0 * pi)) * 2.0 / 3.0result += (160.0 * math.sin(y / 12.0 * pi) + 320 * math.sin(y * pi / 30.0)) * 2.0 / 3.0return resultdef 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.0result += (20.0 * math.sin(x * pi) + 40.0 * math.sin(x / 3.0 * pi)) * 2.0 / 3.0result += (150.0 * math.sin(x / 12.0 * pi) + 300.0 * math.sin(x / 30.0 * pi)) * 2.0 / 3.0return resultdlat = 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 * magicsqrtmagic = 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.01threshold = 0.0000000001dlat, dlon = initdelta, initdeltamlat = gcj_lat - dlatmlon = gcj_lon - dlonplat = gcj_lat + dlatplon = gcj_lon + dlonwgs_lat, wgs_lon, i = 0, 0, 0while True:wgs_lat = (mlat + plat) / 2wgs_lon = (mlon + plon) / 2tmplat, tmplon = wgs_corrected(wgs_lat, wgs_lon)tmplat, tmplon = float(tmplat), float(tmplon)dlat = tmplat - gcj_latdlon = tmplon - gcj_lonif (math.fabs(dlat) < threshold) and (math.fabs(dlon) < threshold):breakif dlat > 0:plat = wgs_latelse:mlat = wgs_latif dlon > 0:plon = wgs_lonelse:mlon = wgs_loni += 1if i > 10000:breakwgs_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.0result += (20.0 * math.sin(y * pi) + 40.0 * math.sin(y / 3.0 * pi)) * 2.0 / 3.0result += (160.0 * math.sin(y / 12.0 * pi) + 320 * math.sin(y * pi / 30.0)) * 2.0 / 3.0return resultdef 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.0result += (20.0 * math.sin(x * pi) + 40.0 * math.sin(x / 3.0 * pi)) * 2.0 / 3.0result += (150.0 * math.sin(x / 12.0 * pi) + 300.0 * math.sin(x / 30.0 * pi)) * 2.0 / 3.0return resultdlat = 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 * magicsqrtmagic = 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_lonif __name__ == '__main__':main()總結
以上是生活随笔為你收集整理的经纬度转换(go/python/rust)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: android 常用框架整理
- 下一篇: 有关于python的论文_有关pytho