Python實(shí)現(xiàn)常見坐標(biāo)系的相互轉(zhuǎn)換

更新時(shí)間：2023年02月10日 16:01:01 作者：hhhSir'blog

WGS84坐標(biāo)系、GCJ02坐標(biāo)系、BD09坐標(biāo)系和Web?墨卡托投影坐標(biāo)系是我們常見的四個(gè)坐標(biāo)系。這篇文章為大家整理了這四個(gè)坐標(biāo)系之間相互轉(zhuǎn)換的方法，需要的可以參考一下

一、背景

主流被使用的地理坐標(biāo)系并不統(tǒng)一，導(dǎo)致我們從不同平臺(tái)下載的數(shù)據(jù)由于坐標(biāo)系的差異往往對(duì)不齊。這個(gè)現(xiàn)象在多源數(shù)據(jù)處理的時(shí)候往往很常見，因此需要進(jìn)行坐標(biāo)轉(zhuǎn)換。

簡單介紹一下幾種常見的坐標(biāo)系：

WGS84坐標(biāo)系：即地球坐標(biāo)系（World Geodetic System），國際上通用的坐標(biāo)系。設(shè)備包含的GPS芯片或者北斗芯片獲取的經(jīng)緯度一般都是為WGS84地理坐標(biāo)系，目前谷歌地圖采用的是WGS84坐標(biāo)系（中國范圍除外）。

GCJ02坐標(biāo)系：GCJ-02是由中國國家測繪局（G表示Guojia國家，C表示Cehui測繪，J表示Ju局）制訂的地理信息系統(tǒng)的坐標(biāo)系統(tǒng)。由WGS84坐標(biāo)系經(jīng)加密后的坐標(biāo)系。谷歌中國采用的GCJ02地理坐標(biāo)系。也稱：火星坐標(biāo)系。

BD09坐標(biāo)系：即百度坐標(biāo)系，GCJ02坐標(biāo)系經(jīng)加密后的坐標(biāo)系。

Web 墨卡托投影坐標(biāo)系：也稱web墨卡托，是如今主流的Web地圖使用的坐標(biāo)系，如國外的 Google Maps，OpenStreetMap，Bing Map，ArcGIS 和 Heremaps 等，國內(nèi)的百度地圖、高德地圖、騰訊地圖和天地圖等也是基于Web墨卡托，與地理坐標(biāo)系不同，投影坐標(biāo)系的單位是m（由于國內(nèi)政策的原因，國內(nèi)地圖會(huì)有加密要求，一般有兩種情況，一種是在 Web墨卡托的基礎(chǔ)上經(jīng)過國家標(biāo)準(zhǔn)加密的國標(biāo)02坐標(biāo)系，熟稱“火星坐標(biāo)系”；另一種是在國標(biāo)的02坐標(biāo)系下進(jìn)一步進(jìn)行加密，如百度地圖的BD09坐標(biāo)系）。

墨卡托投影的“等角”特性，保證了對(duì)象的形狀的不變行，正方形的物體投影后不會(huì)變?yōu)殚L方形。“等角”也保證了方向和相互位置的正確性，因此在航海和航空中常常應(yīng)用，而Google們?cè)谟?jì)算人們查詢地物的方向時(shí)不會(huì)出錯(cuò)。

二、代碼

# -*- coding: utf-8 -*-

import math
import pandas as pd
import os

# WGS84、GCJ02（火星坐標(biāo)系）、BD09（百度坐標(biāo)系）以及百度地圖中保存矢量信息的web墨卡托
class LngLatTransfer():
    
    def __init__(self):
        self.x_pi = 3.14159265358979324 * 3000.0 / 180.0
        self.pi = math.pi  # π
        self.a = 6378245.0  # 長半軸
        self.es = 0.00669342162296594323  # 偏心率平方
        pass

    def GCJ02_to_BD09(self, gcj_lng, gcj_lat):
        """
        實(shí)現(xiàn)GCJ02向BD09坐標(biāo)系的轉(zhuǎn)換
        :param lng: GCJ02坐標(biāo)系下的經(jīng)度
        :param lat: GCJ02坐標(biāo)系下的緯度
        :return: 轉(zhuǎn)換后的BD09下經(jīng)緯度
        """
        z = math.sqrt(gcj_lng * gcj_lng + gcj_lat * gcj_lat) + 0.00002 * math.sin(gcj_lat * self.x_pi)
        theta = math.atan2(gcj_lat, gcj_lng) + 0.000003 * math.cos(gcj_lng * self.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(self, bd_lng, bd_lat):
        '''
        實(shí)現(xiàn)BD09坐標(biāo)系向GCJ02坐標(biāo)系的轉(zhuǎn)換
        :param bd_lng: BD09坐標(biāo)系下的經(jīng)度
        :param bd_lat: BD09坐標(biāo)系下的緯度
        :return: 轉(zhuǎn)換后的GCJ02下經(jīng)緯度
        '''
        x = bd_lng - 0.0065
        y = bd_lat - 0.006
        z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * self.x_pi)
        theta = math.atan2(y, x) - 0.000003 * math.cos(x * self.x_pi)
        gcj_lng = z * math.cos(theta)
        gcj_lat = z * math.sin(theta)
        return gcj_lng, gcj_lat


    def WGS84_to_GCJ02(self, lng, lat):
        '''
        實(shí)現(xiàn)WGS84坐標(biāo)系向GCJ02坐標(biāo)系的轉(zhuǎn)換
        :param lng: WGS84坐標(biāo)系下的經(jīng)度
        :param lat: WGS84坐標(biāo)系下的緯度
        :return: 轉(zhuǎn)換后的GCJ02下經(jīng)緯度
        '''
        dlat = self._transformlat(lng - 105.0, lat - 35.0)
        dlng = self._transformlng(lng - 105.0, lat - 35.0)
        radlat = lat / 180.0 * self.pi
        magic = math.sin(radlat)
        magic = 1 - self.es * magic * magic
        sqrtmagic = math.sqrt(magic)
        dlat = (dlat * 180.0) / ((self.a * (1 - self.es)) / (magic * sqrtmagic) * self.pi)
        dlng = (dlng * 180.0) / (self.a / sqrtmagic * math.cos(radlat) * self.pi)
        gcj_lng = lat + dlat
        gcj_lat = lng + dlng
        return gcj_lng, gcj_lat


    def GCJ02_to_WGS84(self, gcj_lng, gcj_lat):
        '''
        實(shí)現(xiàn)GCJ02坐標(biāo)系向WGS84坐標(biāo)系的轉(zhuǎn)換
        :param gcj_lng: GCJ02坐標(biāo)系下的經(jīng)度
        :param gcj_lat: GCJ02坐標(biāo)系下的緯度
        :return: 轉(zhuǎn)換后的WGS84下經(jīng)緯度
        '''
        dlat = self._transformlat(gcj_lng - 105.0, gcj_lat - 35.0)
        dlng = self._transformlng(gcj_lng - 105.0, gcj_lat - 35.0)
        radlat = gcj_lat / 180.0 * self.pi
        magic = math.sin(radlat)
        magic = 1 - self.es * magic * magic
        sqrtmagic = math.sqrt(magic)
        dlat = (dlat * 180.0) / ((self.a * (1 - self.es)) / (magic * sqrtmagic) * self.pi)
        dlng = (dlng * 180.0) / (self.a / sqrtmagic * math.cos(radlat) * self.pi)
        mglat = gcj_lat + dlat
        mglng = gcj_lng + dlng
        lng = gcj_lng * 2 - mglng
        lat = gcj_lat * 2 - mglat
        return lng, lat


    def BD09_to_WGS84(self, bd_lng, bd_lat):
        '''
        實(shí)現(xiàn)BD09坐標(biāo)系向WGS84坐標(biāo)系的轉(zhuǎn)換
        :param bd_lng: BD09坐標(biāo)系下的經(jīng)度
        :param bd_lat: BD09坐標(biāo)系下的緯度
        :return: 轉(zhuǎn)換后的WGS84下經(jīng)緯度
        '''
        lng, lat = self.BD09_to_GCJ02(bd_lng, bd_lat)
        return self.GCJ02_to_WGS84(lng, lat)


    def WGS84_to_BD09(self, lng, lat):
        '''
        實(shí)現(xiàn)WGS84坐標(biāo)系向BD09坐標(biāo)系的轉(zhuǎn)換
        :param lng: WGS84坐標(biāo)系下的經(jīng)度
        :param lat: WGS84坐標(biāo)系下的緯度
        :return: 轉(zhuǎn)換后的BD09下經(jīng)緯度
        '''
        lng, lat = self.WGS84_to_GCJ02(lng, lat)
        return self.GCJ02_to_BD09(lng, lat)


    def _transformlat(self, 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 * self.pi) + 20.0 *
                math.sin(2.0 * lng * self.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(lat * self.pi) + 40.0 *
                math.sin(lat / 3.0 * self.pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(lat / 12.0 * self.pi) + 320 *
                math.sin(lat * self.pi / 30.0)) * 2.0 / 3.0
        return ret


    def _transformlng(self, 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 * self.pi) + 20.0 *
                math.sin(2.0 * lng * self.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(lng * self.pi) + 40.0 *
                math.sin(lng / 3.0 * self.pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(lng / 12.0 * self.pi) + 300.0 *
                math.sin(lng / 30.0 * self.pi)) * 2.0 / 3.0
        return ret

    def WGS84_to_WebMercator(self, lng, lat):
        '''
        實(shí)現(xiàn)WGS84向web墨卡托的轉(zhuǎn)換
        :param lng: WGS84經(jīng)度
        :param lat: WGS84緯度
        :return: 轉(zhuǎn)換后的web墨卡托坐標(biāo)
        '''
        x = lng * 20037508.342789 / 180
        y = math.log(math.tan((90 + lat) * self.pi / 360)) / (self.pi / 180)
        y = y * 20037508.34789 / 180
        return x, y

    def WebMercator_to_WGS84(self, x, y):
        '''
        實(shí)現(xiàn)web墨卡托向WGS84的轉(zhuǎn)換
        :param x: web墨卡托x坐標(biāo)
        :param y: web墨卡托y坐標(biāo)
        :return: 轉(zhuǎn)換后的WGS84經(jīng)緯度
        '''
        lng = x / 20037508.34 * 180
        lat = y / 20037508.34 * 180
        lat = 180 / self.pi * (2 * math.atan(math.exp(lat * self.pi / 180)) - self.pi / 2)
        return lng, lat
    
if __name__=='__main__':
    fileName = r'F:\武漢軌跡數(shù)據(jù)\交通事故(2018年)\accidentFileLocations.csv'
    transData = pd.read_csv(fileName, engine='python')
    transData["WGS84lng"] = None
    transData["WGS84lat"] = None
    # 火星坐標(biāo)系 轉(zhuǎn)換為 wgs84坐標(biāo)系：GCJ02_to_WGS84 (lng, lat)
    handler = LngLatTransfer()
    transData[["WGS84lng", "WGS84lat"]] = transData.apply(lambda x : handler.GCJ02_to_WGS84(x["LON"], x["LAT"]), axis = 1, result_type="expand")
    os.chdir(r'F:\武漢軌跡數(shù)據(jù)\交通事故(2018年)')
    transData.to_csv("LoacationTransTest.csv", index = False)

直接貼個(gè)代碼，具體怎么實(shí)現(xiàn)和怎么使用的就很清楚了，不多言。代碼來源，而且真心實(shí)推GIS專業(yè)的學(xué)生多看看這個(gè)老哥的blog，大神。

上面代碼的邏輯可以用這張圖來表示，是不是更加清楚了。