python 坐标转换

用于地图显示，坐标转换

import math
from typing import Tuple
import numpy as np
import numbers

PI = math.pi
R = 20037508.34
K0 = 0.9996
FE = 500000
FN = 0
a = 6378137
b = 6356752.3142
e = 0.0818192
e2 = pow(e, 2)
e4 = pow(e, 4)
e6 = pow(e, 6)
x_pi = 3.14159265358979324 * 3000.0 / 180.0


def convert_to_mercator(lon, lat):
    """
    经纬度坐标转换墨卡托坐标

    Args:
        lon: 经度
        lat: 纬度

    Returns: 墨卡托坐标

    """

    x = lon * R / 180
    y = np.log(np.tan((90 + lat) * PI / 360)) / (PI / 180)
    y = y * R / 180
    if isinstance(x, numbers.Real) and isinstance(y, numbers.Real):
        return int(x), int(y)
    else:
        x = x.astype(int)
        y = y.astype(int)
        return x, y


def convert_to_coordinate(x, y):
    """
    墨卡托坐标转经纬度坐标

    Args:
        x: 墨卡托坐标x
        y: 墨卡托坐标y

    Returns: 经纬度坐标

    """
    dx = x / float(R) * 180
    dy = y / float(R) * 180
    dy = 180 / float(PI) * (2 * np.arctan(np.exp(dy * float(PI) / 180)) - float(PI) / 2)
    return dx, dy


def bd09_to_gcj02(bd09_lon, bd09_lat):
    """
    bd09转gcj02

    Args:
        bd09_lon: bd09坐标系经度
        bd09_lat: bd09坐标系纬度

    Returns: gcj02坐标系

    """
    x = bd09_lon - 0.0065
    y = bd09_lat - 0.006
    z = np.sqrt(x * x + y * y) - 0.00002 * np.sin(y * x_pi)
    theta = np.arctan2(y, x) - 0.000003 * np.cos(x * x_pi)
    gcj02_lon = z * np.cos(theta)
    gcj02_lat = z * np.sin(theta)
    return gcj02_lon, gcj02_lat


def gcj02_to_wgs84(gcj02_lon, gcj02_lat):
    """
    gcj02转wgs84

    Args:
        gcj02_lon: gcj02坐标系经度
        gcj02_lat: gcj02坐标系纬度

    Returns: wgs84坐标系

    """
    if isinstance(gcj02_lon, numbers.Real) and isinstance(gcj02_lat, numbers.Real):
        lon = float(gcj02_lon)
        lat = float(gcj02_lat)
    else:
        lon = gcj02_lon.astype(float)
        lat = gcj02_lat.astype(float)

    a = 6378245.0
    ee = 0.00669342162296594323

    x = lon - 105.0
    y = lat - 35.0
    # 经度
    d_lon = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * np.sqrt(abs(x))
    d_lon += (20.0 * np.sin(6.0 * x * PI) + 20.0 * np.sin(2.0 * x * PI)) * 2.0 / 3.0
    d_lon += (20.0 * np.sin(x * PI) + 40.0 * np.sin(x / 3.0 * PI)) * 2.0 / 3.0
    d_lon += (150.0 * np.sin(x / 12.0 * PI) + 300.0 * np.sin(x / 30.0 * PI)) * 2.0 / 3.0
    # 纬度
    d_lat = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * np.sqrt(abs(x))
    d_lat += (20.0 * np.sin(6.0 * x * PI) + 20.0 * np.sin(2.0 * x * PI)) * 2.0 / 3.0
    d_lat += (20.0 * np.sin(y * PI) + 40.0 * np.sin(y / 3.0 * PI)) * 2.0 / 3.0
    d_lat += (160.0 * np.sin(y / 12.0 * PI) + 320 * np.sin(y * PI / 30.0)) * 2.0 / 3.0
    rad_lat = lat / 180.0 * PI
    magic = np.sin(rad_lat)
    magic = 1 - ee * magic * magic
    sqrt_magic = np.sqrt(magic)
    d_lat = (d_lat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * PI)
    d_lon = (d_lon * 180.0) / (a / sqrt_magic * np.cos(rad_lat) * PI)
    wgs84_lon = lon - d_lon
    wgs84_lat = lat - d_lat
    if isinstance(gcj02_lon, numbers.Real) and isinstance(y, numbers.Real):
        return float(wgs84_lon), float(wgs84_lat)
    else:
        wgs84_lon = wgs84_lon.astype(float)
        wgs84_lat = wgs84_lat.astype(float)
        return wgs84_lon, wgs84_lat

可用于dataframe，某一列批量处理，提高速度