上一篇刚学会了如何获取地址所在经纬度以及经纬度对应地址,于是信心满满的准备在老板面前露一手,准备花点心思做个漂亮的可视化地图放在报告亮眼的位置。
当地图跑出来的那一刻顿时傻眼了,卧槽这些点定位的位置明显不对呀,酒店直定位到湖里了,这让老板看到你说后果严重不。
其实是因为百度地图的经纬度和高德地图经纬度所使用的坐标系编码不同,所以如果你制图的软件或平台是基于高德地图服务的话,那么就需要使用对应的高德标准的经纬度来进行绘图,同理如果你用的服务是基于百度地图的,必须使用百度标准的经纬度坐标,如果拿到的是高德坐标则必须经过算法转为百度标准。
今天这篇内容就给大家分享如何使用R和Python对百度坐标系和高德坐标系的经纬度进行互转,解决地图绘制最后一道坎儿。
内容同样是两个模块四小节,使用R语言和Python分别进行百度经纬度转高德经纬度以及高德标准转百度,这样我们在经纬度获取和转化处理上的技能掌握就比较系统了,不会再受制于工具和平台服务的标准差异而苦恼。
一、R语言方案
1)百度经纬度转腾讯&高德
library("leaflet")library('baidumap')library('ggplot2')library('ggmap')dt1 pipe(#百度转腾讯&高德bMapTransQQMap x_pi = 3.14159265358979324 * 3000.0 /180.0 x = lng - 0.0065 y = lat - 0.006 z = sqrt(x^2 + y^2) - 0.00002 * sin(y * x_pi) theta = atan2(y,x) - 0.000003 * cos(x * x_pi) lngs = z * cos(theta) lats = z * sin(theta) return(data.frame(lng,lat,lngs,lats)) }result
使用高德地图服务来呈现百度坐标系下的经纬度,可以看到偏移量非常明显。
leaflet(result) %>% setView(116.2938,40.00939, zoom = 13) %>% addTiles('http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',tileOptions(tileSize=256, minZoom=9, maxZoom=17),group="高德地图") %>% addCircleMarkers(lng = ~ lng,lat = ~lat,stroke = FALSE,fillOpacity = 1,radius =8)
使用高德地图服务来呈现经百度坐标系转换高德坐标系后的经纬度,可以看到位置基本已经还原了。
leaflet(result) %>% setView(116.2938,40.00939, zoom = 13) %>% addTiles('http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',tileOptions(tileSize=256, minZoom=9, maxZoom=17),group="高德地图") %>% addCircleMarkers(lng = ~ lngs,lat = ~lats,stroke = FALSE,fillOpacity = 1,radius =8)
2)高德坐标系转百度坐标系
# 腾讯&高德地图转百度qqMapTransBMap x_pi = 3.14159265358979324 * 3000.0 / 180.0 x = lngs y = lats z = sqrt(x^2 + y^2) + 0.00002 * sin(y * x_pi) theta = atan2(y,x) + 0.000003 * cos(x * x_pi) lng_b = z * cos(theta) + 0.0065 lat_b = z * sin(theta) + 0.006 return(data.frame(lngs,lats,lng_b,lat_b))}result_b
使用百度地图服务呈现高德坐标系下的经纬度,偏移情况仍然很严重。
options(baidumap.key = '***************')# 这里启用的baidumap服务,后台调用百度地图api,# 需要使用百度api开发秘钥,需自己申请bjMap <- getBaiduMap(location = c(116.284028,40.001732), width = 800, height =800, zoom = 16, scale = 2 , messaging = TRUE)ggmap(bjMap) + geom_point(data = result_b, aes(x = lngs, y = lats),shape = 21 , col = 'white',fill = 'red',size = 5) + theme_nothing()
使用百度地图服务呈现经高德坐标系转百度坐标系后的经纬度,可以看到为止基本已经还原到真实位置,误差相对较小。
ggmap(bjMap) + geom_point(data = result_b, aes(x = lng_b, y = lat_b),shape = 21 , col = 'white',fill = 'red',size = 5) + theme_nothing()
二、Python实现方案
3)百度坐标系转高德坐标系
import pandas as pdimport numpy as npimport foliumfrom folium import pluginsdf=pd.read_clipboard() def bMapTransQQMap(lng,lat): x_pi = 3.14159265358979324 * 3000.0 /180.0 x = np.array(lng) - 0.0065 y = np.array(lat) - 0.006 z = np.sqrt(np.power(x,2) + np.power(y,2)) - 0.00002 * np.sin(y * x_pi) theta = np.arctan2(y,x) - 0.000003 * np.cos(x * x_pi) lngs = z * np.cos(theta) lats = z * np.sin(theta) return(pd.DataFrame({'lng':lng,'lat':lat,'lngs':lngs,'lats':lats})) result = bMapTransQQMap(df['lon'],df['lat'])
使用folium库(底层调用leaflet服务)结合高德地图服务来呈现百度标准的的经纬度,误差依然很大。
lng = np.array(result["lng"],dtype=float)lat = np.array(result["lat"],dtype=float)data1 = [(lat[i],lng[i]) for i in range(len(result))]map_osm = folium.Map( location=[39.996710,116.281012], zoom_start=115, tiles='http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}', attr="© 高德地图" )marker_cluster = plugins.MarkerCluster().add_to(map_osm) for i in data1: folium.Marker(i).add_to(marker_cluster) display(map_osm)
百度坐标系经转化为高德坐标系之后,位置的精确度基本还原,明显改善。
lng = np.array(result["lngs"],dtype=float)lat = np.array(result["lats"],dtype=float)data2 = [(lat[i],lng[i]) for i in range(len(result))]map_osm = folium.Map( location=[39.996710,116.281012], zoom_start=115, tiles='http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}', attr="© 高德地图" )marker_cluster = plugins.MarkerCluster().add_to(map_osm) for i in data2: folium.Marker(i).add_to(marker_cluster) display(map_osm)
4)高德坐标系转百度坐标系
#腾讯&高德地图转百度def qqMapTransBMap(lngs,lats): x_pi = 3.14159265358979324 * 3000.0 / 180.0 x = np.array(lngs) y = np.array(lats ) z = np.sqrt(np.power(x,2) + np.power(y,2)) + 0.00002 * np.sin(y * x_pi) theta = np.arctan2(y,x) + 0.000003 * np.cos(x * x_pi) lng_b = z * np.cos(theta) + 0.0065 lat_b = z * np.sin(theta) + 0.006 return(pd.DataFrame({'lngs':lngs,'lats':lats,'lng_b':lng_b,'lat_b':lat_b}))result2 = qqMapTransBMap(result['lngs'],result['lats'])
至此,两种技术方案下的所有类型转换均以搞定,虽然至今还没明白算法里面的具体参数到底是啥意思,只是照葫芦画瓢把一篇博客上的java代码翻译了过来(再次感谢原作者提供的java代码),不过原理不懂没有关系,语法能看懂就OK了。
参考资料:
https://www.cnblogs.com/wrld/p/10845870.html