广州精细地铁线路图轨迹图可视化.py

代码实现了在python下利用百度api获取地铁站点及地铁精细路线（轨迹）的墨卡托坐标数据并转化为WGS-84坐标系下的经纬度，并使用 plotly+mapbox 绘制地铁路线图

import numpy as np
import requests
import math
import time
import plotly.offline as py
import plotly.graph_objs as go
from plotly.offline import init_notebook_mode

init_notebook_mode（connected=True）


PI = math.pi

def _transformlat（coordinates）:
    lng = coordinates[ :  0] - 105
    lat = coordinates[ :  1] - 35
    ret = -100 + 2 * lng + 3 * lat + 0.2 * lat * lat + \
          0.1 * lng * lat + 0.2 * np.sqrt（np.fabs（lng））
    ret += （20 * np.sin（6 * lng * PI） + 20 *
            np.sin（2 * lng * PI）） * 2 / 3
    ret += （20 * np.sin（lat * PI） + 40 *
            np.sin（lat / 3 * PI）） * 2 / 3
    ret += （160 * np.sin（lat / 12 * PI） + 320 *
            np.sin（lat * PI / 30.0）） * 2 / 3
    return ret


def _transformlng（coordinates）:
    lng = coordinates[ :  0] - 105
    lat = coordinates[ :  1] - 35
    ret = 300 + lng + 2 * lat + 0.1 * lng * lng + \
          0.1 * lng * lat + 0.1 * np.sqrt（np.fabs（lng））
    ret += （20 * np.sin（6 * lng * PI） + 20 *
            np.sin（2 * lng * PI）） * 2 / 3
    ret += （20 * np.sin（lng * PI） + 40 *
            np.sin（lng / 3 * PI）） * 2 / 3
    ret += （150 * np.sin（lng / 12 * PI） + 300 *
            np.sin（lng / 30 * PI）） * 2 / 3
    return ret

def gcj02_to_wgs84（coordinates）:
    “““
    GCJ-02转WGS-84
    :param coordinates: GCJ-02坐标系的经度和纬度的numpy数组
    :returns: WGS-84坐标系的经度和纬度的numpy数组
    “““
    ee = 0.006693421622965943  # 偏心率平方
    a = 6378245  # 长半轴
    lng = coordinates[ :  0]
    lat = coordinates[ :  1]
    is_in_china = （lng > 73.66） & （lng < 135.05） & （lat > 3.86） & （lat < 53.55）
    _transform = coordinates[is_in_china]  #只对不在国内的坐标做偏移
    
    dlat = _transformlat（_transform）
    dlng = _transformlng（_transform）
    radlat = _transform[ :  1] / 180 * PI
    magic = np.sin（radlat）
    magic = 1 - ee * magic * magic
    sqrtmagic = np.sqrt（magic）
    dlat = （dlat * 180.0） / （（a * （1 - ee）） / （magic * sqrtmagic） * PI）
    dlng = （dlng * 180.0） / （a / sqrtmagic * np.cos（radlat） * PI）
    mglat = _transform[ :  1] + dlat
    mglng = _transform[ :  0] + dlng
    coordinates[is_in_china] = np.array（[_transform[ :  0] * 2 - mglng _transform[ :  1] * 2 - mglat]）.T
    return coordinates

def bd09_to_gcj02（coordinates）:
    “““
    BD-09转GCJ-02
    :param coordinates: BD-09坐标系的经度和纬度的numpy数组
    :returns: GCJ-02坐标系的经度和纬度的numpy数组
    “““
    x_pi = PI * 3000 / 180
    x = coordinates[ :  0] - 0.0065
    y = coordinates[ :  1] - 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）
    lng = z * np.cos（theta）
    lat = z * np.sin（theta）
    coordinates = np.array（[lng lat]）.T
    return coordinates

def bd09_to_wgs84（coordinates）:
    “““
    BD-09转WGS-84
    :param coordinates: BD-09坐标系的经度和纬度的numpy数组
    :returns: WGS-84坐标系的经度和纬度的numpy数组
    “““
    return gcj02_to_wgs84（bd09_to_gcj02（coordinates））

def mercator_to_bd09（mercator）:
    “““