标签：GIS Python sin lat return wgsLat pi 坐标系 math

#!/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)

标签：GIS,Python,sin,lat,return,wgsLat,pi,坐标系,math
来源： https://blog.csdn.net/linzi1994/article/details/115409765

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。