1. 程式人生 > >python2介面:WGS-84轉BD座標系

python2介面:WGS-84轉BD座標系

一、背景 關於python 解析 PL2303裝置的原始資料得到WGS座標系的經緯度的文章檢視上一篇;

二、直接上乾貨,倆種方法,若有其他的請賜教 1.依賴百度API(耗網操作)

    def wgs84tobd09ByNetwork(self, lon, lat):
    	#傳入你的ak
        queryStr = '/geoconv/v1/?coords={},{}&from=1&to=5&ak=yourak'.format(lon, lat)
        # 對queryStr進行轉碼,safe內的保留字元不轉換
        encodedStr = urllib.quote(queryStr, safe="/:=&?#+!$,;'@()*[]")
        # 在最後直接追加上yoursk
        rawStr = encodedStr + 'fFyAglwnmgxnYI6c---------------------------'
        # md5計算出的sn值
        my_sn = hashlib.md5(urllib.quote_plus(rawStr)).hexdigest()
        url = 'http://api.map.baidu.com' + queryStr + "&sn=" + my_sn
        res = requests.get(url)
        # get收到的內容
        json_str = res.content
        # print json_str
        dictData = json.loads(json_str)
        print(dictData)
        print (dictData["result"][0]["x"])
        print (dictData["result"][0]["y"])
        return dictData["result"][0]["x"], dictData["result"][0]["y"]

2.演算法(誤差極小)

    def transformlat(self,lng, lat):
        pi = 3.1415926535897932384626  # π

        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(self,lng, lat):
        pi = 3.1415926535897932384626  # π
        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(self,lng, lat):
        """
        判斷是否在國內,不在國內不做偏移
        :param lng:
        :param lat:
        :return:
        """
        if lng < 72.004 or lng > 137.8347:
            return True
        if lat < 0.8293 or lat > 55.8271:
            return True
        return False
        
    def wgs84tobd09(self,lng, lat):
        x_pi = 3.14159265358979324 * 3000.0 / 180.0
        pi = 3.1415926535897932384626   # π
        a = 6378245.0  # 長半軸
        ee = 0.00669342162296594323      # 扁率
        if self.out_of_china(lng, lat):  # 判斷是否在國內
            return lng, lat
        dlat = self.transformlat(lng - 105.0, lat - 35.0)
        dlng = self.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
        z = math.sqrt(mglng * mglng + mglat * mglat) + 0.00002 * math.sin(mglat * x_pi)
        theta = math.atan2(mglat, mglng) + 0.000003 * math.cos(mglng * x_pi)
        bd_lng = z * math.cos(theta) + 0.0065
        bd_lat = z * math.sin(theta) + 0.006
        print('轉換成百度座標:%s,%s'%(bd_lng,bd_lat))
        return bd_lng,bd_lat