Python实现GIS坐标系的转换

曾经终败给现在 2022-10-22 02:47 229阅读 0赞

#!/usr/bin/env python3.6.5
    # -*- coding: UTF-8 -*-
    """ Author: ReddingtonLin Date: 2020/4/1 14:23 docs: GPS坐标转换： WGS-84：是国际标准，GPS坐标（Google Earth使用、或者GPS模块） GCJ-02：中国坐标偏移标准，Google Map、高德、腾讯使用 BD-09：百度坐标偏移标准，Baidu Map使用 """
    
    import math
    
    
    def transformLat(x, y):
        ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
        return ret
    
    
    def transformLon(x, y):
        ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
        return ret
    
    
    def delta(lng, lat):
        a = 6378245.0
        # a: 卫星椭球坐标投影到平面地图坐标系的投影因子
        ee = 0.00669342162296594323
        # ee: 椭球的偏心率
        dLat = transformLat(lng - 105.0, lat - 35.0)
        dLon = transformLon(lng - 105.0, lat - 35.0)
        radLat = lat / 180.0 * math.pi
        magic = math.sin(radLat)
        magic = 1 - ee * magic * magic
        sqrtMagic = math.sqrt(magic)
        dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * math.pi)
        dLon = (dLon * 180.0) / (a / sqrtMagic * math.cos(radLat) * math.pi)
        return dLon, dLat
    
    
    def wgs2gcj(wgsLng, wgsLat):
        """ WGS-84转成GCJ-02 """
        if outOfChina(wgsLng, wgsLat):
            print("The latitude or longitude is out of China!")
            return wgsLng, wgsLat
        lng, lat = delta(wgsLng, wgsLat)
        return wgsLng + lng, wgsLat + lat
    
    
    def gcj2wgs_rough(gcjLon, gcjLat):
        """ GCJ-02 转 WGS-84 粗略版 """
        if outOfChina(gcjLon, gcjLat):
            print("The latitude or longitude is out of China!")
            return gcjLon, gcjLat
        lng, lat = delta(gcjLon, gcjLat)
        return gcjLon - lng, gcjLat - lat
    
    
    def gcj2wgs_accurate(gcjLon, gcjLat):
        """ GCJ-02 转 WGS-84 精确版 """
        initDelta = 0.01
        threshold = 0.000000001
        dLat = initDelta
        dLon = initDelta
        mLat = gcjLat - dLat
        mLon = gcjLon - dLon
        pLat = gcjLat + dLat
        pLon = gcjLon + dLon
        wgsLat = 0
        wgsLon = 0
        i = 0
        while 1:
            wgsLat = (mLat + pLat) / 2
            wgsLon = (mLon + pLon) / 2
            lon, lat = gcj2wgs_rough(wgsLon, wgsLat)
            dLat = lat - gcjLat
            dLon = lon - gcjLon
            if (abs(dLat) < threshold) and (abs(dLon) < threshold):
                break
            if dLat > 0:
                pLat = wgsLat
            else:
                mLat = wgsLat
            if dLon > 0:
                pLon = wgsLon
            else:
                mLon = wgsLon
            if ++i > 10000:
                break
        return wgsLon, wgsLat
    
    
    def gcj2bd(gcjLon, gcjLat):
        """ GCJ-02 转 BD-09 """
        x_pi = math.pi * 3000.0 / 180.0
        x = gcjLon
        y = gcjLat
        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)
        bdLon = z * math.cos(theta) + 0.0065
        bdLat = z * math.sin(theta) + 0.006
        return bdLon, bdLat
    
    
    def bd2gcj(bdLon, bdLat):
        """ BD-09 转 GCJ-02 """
        x_pi = math.pi * 3000.0 / 180.0
        x = bdLon - 0.0065
        y = bdLat - 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)
        gcjLon = z * math.cos(theta)
        gcjLat = z * math.sin(theta)
        return gcjLon, gcjLat
    
    
    def wgs2mercator(wgsLon, wgsLat):
        """ WGS-84 to Web mercator mercatorLat -> y mercatorLon -> x """
        x = wgsLon * 20037508.34 / 180.
        y = math.log(math.tan((90. + wgsLat) * math.pi / 360)) / (math.pi / 180)
        y = y * 20037508.34 / 180.
        return x, y
    
    
    def mercator2wgs(mercatorLon, mercatorLat):
        """ Web mercator to WGS-84 mercatorLat -> y mercatorLon -> x """
        x = mercatorLon / 20037508.34 * 180
        y = mercatorLat / 20037508.34 * 180
        y = 180 / math.pi * (2 * math.atan(math.exp(y * math.pi / 180.)) - math.pi / 2)
        return x, y
    
    
    def outOfChina(lng, lat):
        """ 判断是否在中国范围外 """
        if lng < 72.004 or lng > 137.8347:
            return True
        if lat < 0.8293 or lat > 55.8271:
            return True
        return False
    
    
    def haversine(lon1, lat1, lon2, lat2):
        """ :param: 纬度1，经度1，纬度2，经度2（十进制度数） :return: 二个坐标之间的距离（单位米） Calculate the great circle distance between two points on the earth (specified in decimal degrees) """
        # 将十进制度数转化为弧度
        lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
        # haversine公式
        dlon = lon2 - lon1
        dlat = lat2 - lat1
        a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
        c = 2 * math.asin(math.sqrt(a))
        r = 6371  # 地球平均半径，单位为公里
        return c * r * 1000
    
    
    def wgs2bd(lon, lat):
        newlon, newlat = wgs2gcj(lon, lat)
        return gcj2bd(newlon, newlat)
    
    
    if __name__ == '__main__':
        a = [(116.32232,39.79082),
            (116.328371,39.790688),
            (116.329444,39.782592),
            (116.322449,39.783251),
            (116.32232,39.79082)]
        r = []
        for i in a:
            lon, lat = i
            lon_wgs84, lat_wgs84 = gcj2wgs_rough(lon, lat)
            r.append((lon_wgs84, lat_wgs84))
        print(r)