zoukankan      html  css  js  c++  java
  • 【615】国内国外经纬度坐标转换

    参考:python3实现GPS经纬度坐标(WGS84)国测局火星坐标(GCJ02)百度坐标(BD09)相互转换

    参考:地图,GPS位置地图坐标系:WGS-84(GPS)、GCJ-02(Google地图)、BD-09(百度地图),OpenGIS


      前言:下载的国外提供的数据集,对应的是国内的建筑数据,但是使用的是 GPS 坐标,通过 folium 加载国内地图显示出现位移,因此需要转为国内的 火星坐标,具体方法如下:

      在进行坐标转换前先来了解一下目前的的坐标体系分类:

    • 一是GPS坐标,也即WGS-84坐标是一个国际的标准,一般卫星导航,原始的GPS设备中的数据都是采用这一坐标系。国外的Google地图、OSM等采用的都是这一坐标。

    • 二是国测局坐标,国测局坐标GCJ-02坐标也叫火星坐标,是国家测绘局为了国家安全在原始坐标的基础上进行偏移得到的坐标,基本国内的电子地图、导航设备都是采用的这一坐标系,如:高德、搜搜、51地图MapABC地图,谷歌中国地图也是。

    • 三是百度坐标,百度坐标BD-09坐标是百度公司出于商业保护在国测局坐标基础上进行的二次加密。

      另外,还有一些其他公司也进行了二次加密形成了自己的坐标如:图吧坐标、搜狗坐标等。

      知道了上述坐标的基本情况后下面具体怎么进行相互转换,直接上代码:

    # -*- coding: utf-8 -*-
    import json
    import requests
    import math
    x_pi = 3.14159265358979324 * 3000.0 / 180.0
    pi = 3.1415926535897932384626  # π
    a = 6378245.0  # 长半轴
    ee = 0.00669342162296594323  # 偏心率平方
    class Geocoding:
        def __init__(self, api_key):
            self.api_key = api_key
        def geocode(self, address):
            """
            利用高德geocoding服务解析地址获取位置坐标
            :param address:需要解析的地址
            :return:
            """
            geocoding = {'s': 'rsv3',
                         'key': self.api_key,
                         'city': '全国',
                         'address': address}
            # geocoding = urllib.urlencode(geocoding)
            # ret = urllib.urlopen("http://restapi.amap.com/v3/geocode/geo{}".format(geocoding))
            url = "http://restapi.amap.com/v3/geocode/geo?"
            ret = requests.get(url, params=geocoding)
            if ret.status_code == 200:
                # res = ret.json()
                # json_obj = json.loads(res)
                json_obj = ret.json()
                if json_obj['status'] == '1' and int(json_obj['count']) >= 1:
                    geocodes = json_obj['geocodes'][0]
                    lng = float(geocodes.get('location').split(',')[0])
                    lat = float(geocodes.get('location').split(',')[1])
                    return [lng, lat]
                else:
                    return None
            else:
                return None
    def gcj02_to_bd09(lng, lat):
        """
        火星坐标系(GCJ-02)转百度坐标系(BD-09)
        谷歌、高德——>百度
        :param lng:火星坐标经度
        :param lat:火星坐标纬度
        :return:
        """
        z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi)
        theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi)
        bd_lng = z * math.cos(theta) + 0.0065
        bd_lat = z * math.sin(theta) + 0.006
        return [bd_lng, bd_lat]
    def bd09_to_gcj02(bd_lon, bd_lat):
        """
        百度坐标系(BD-09)转火星坐标系(GCJ-02)
        百度——>谷歌、高德
        :param bd_lat:百度坐标纬度
        :param bd_lon:百度坐标经度
        :return:转换后的坐标列表形式
        """
        x = bd_lon - 0.0065
        y = 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)
        gg_lng = z * math.cos(theta)
        gg_lat = z * math.sin(theta)
        return [gg_lng, gg_lat]
    def wgs84_to_gcj02(lng, lat):
        """
        WGS84转GCJ02(火星坐标系)
        :param lng:WGS84坐标系的经度
        :param lat:WGS84坐标系的纬度
        :return:
        """
        if out_of_china(lng, lat):  # 判断是否在国内
            return [lng, lat]
        dlat = _transformlat(lng - 105.0, lat - 35.0)
        dlng = _transformlng(lng - 105.0, lat - 35.0)
        radlat = lat / 180.0 * pi
        magic = math.sin(radlat)
        magic = 1 - ee * magic * magic
        sqrtmagic = math.sqrt(magic)
        dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
        dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
        mglat = lat + dlat
        mglng = lng + dlng
        return [mglng, mglat]
    def gcj02_to_wgs84(lng, lat):
        """
        GCJ02(火星坐标系)转GPS84
        :param lng:火星坐标系的经度
        :param lat:火星坐标系纬度
        :return:
        """
        if out_of_china(lng, lat):
            return [lng, lat]
        dlat = _transformlat(lng - 105.0, lat - 35.0)
        dlng = _transformlng(lng - 105.0, lat - 35.0)
        radlat = lat / 180.0 * pi
        magic = math.sin(radlat)
        magic = 1 - ee * magic * magic
        sqrtmagic = math.sqrt(magic)
        dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
        dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
        mglat = lat + dlat
        mglng = lng + dlng
        return [lng * 2 - mglng, lat * 2 - mglat]
    def bd09_to_wgs84(bd_lon, bd_lat):
        lon, lat = bd09_to_gcj02(bd_lon, bd_lat)
        return gcj02_to_wgs84(lon, lat)
    def wgs84_to_bd09(lon, lat):
        lon, lat = wgs84_to_gcj02(lon, lat)
        return gcj02_to_bd09(lon, lat)
    def _transformlat(lng, lat):
        ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 
              0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
        ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
                math.sin(2.0 * lng * pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(lat * pi) + 40.0 *
                math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
                math.sin(lat * pi / 30.0)) * 2.0 / 3.0
        return ret
    def _transformlng(lng, lat):
        ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 
              0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
        ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
                math.sin(2.0 * lng * pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(lng * pi) + 40.0 *
                math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
                math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
        return ret
    def out_of_china(lng, lat):
        """
        判断是否在国内,不在国内不做偏移
        :param lng:
        :param lat:
        :return:
        """
        return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)
    if __name__ == '__main__':
        lng = 116.382997
        lat = 39.915156
        result1 = gcj02_to_bd09(lng, lat)
        result2 = bd09_to_gcj02(lng, lat)
        result3 = wgs84_to_gcj02(lng, lat)
        result4 = gcj02_to_wgs84(lng, lat)
        result5 = bd09_to_wgs84(lng, lat)
        result6 = wgs84_to_bd09(lng, lat)
        g = Geocoding('apikey')  # 这里填写你的高德api的key
        result7 = g.geocode('广东省深圳市南山区')
    print(result1, result2, result3, result4, result5, result6, result7)
    

      对于我的需求,主要是用到了函数 wgs84_to_gcj02 来实现。 

  • 相关阅读:
    线性规划与网络流24题解题报告
    [bzoj4569][SCOI2016]萌萌哒-并查集+倍增
    [bzoj1002]轮状病毒-矩阵树定理
    [bzoj1005][HNOI2008]明明的烦恼-Prufer编码+高精度
    [bzoj3995][SDOI2015]道路修建-线段树
    [bzoj3993][SDOI2015]星际战争-二分+最大流
    [bzoj3994][SDOI2015]约数个数和-数论
    [bzoj3990][SDOI2015]排序-搜索
    [bzoj4518][Sdoi2016]征途-斜率优化
    [bzoj4515][Sdoi2016]游戏-树链剖分+李超线段树
  • 原文地址:https://www.cnblogs.com/alex-bn-lee/p/15032409.html
Copyright © 2011-2022 走看看