
import numpy as np
x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626 
a = 6378245.0  
ee = 0.00669342162296594323  

def gcj02tobd09(lng, lat):
    """
    Convert coordinates from GCJ02 to BD09

    Parameters
    -------
    lng : Series or number
        Longitude
    lat : Series or number
        Latitude

    return
    -------
    lng : Series or number
        Longitude (Converted)
    lat : Series or number
        Latitude (Converted)
    """
    try:
        lng = lng.astype(float)
        lat = lat.astype(float)
    except:
        lng = float(lng)
        lat = float(lat)
    z = np.sqrt(lng * lng + lat * lat) + 0.00002 * np.sin(lat * x_pi)
    theta = np.arctan2(lat, lng) + 0.000003 * np.cos(lng * x_pi)
    bd_lng = z * np.cos(theta) + 0.0065
    bd_lat = z * np.sin(theta) + 0.006
    return bd_lng, bd_lat


def bd09togcj02(bd_lon, bd_lat):
    """
    Convert coordinates from BD09 to GCJ02

    Parameters
    -------
    lng : Series or number
        Longitude
    lat : Series or number
        Latitude

    return
    -------
    lng : Series or number
        Longitude (Converted)
    lat : Series or number
        Latitude (Converted)
    """
    try:
        bd_lon = bd_lon.astype(float)
        bd_lat = bd_lat.astype(float)
    except:
        bd_lon = float(bd_lon)
        bd_lat = float(bd_lat)
    x = bd_lon - 0.0065
    y = bd_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)
    gg_lng = z * np.cos(theta)
    gg_lat = z * np.sin(theta)
    return gg_lng, gg_lat


def wgs84togcj02(lng, lat):
    """
    Convert coordinates from WGS84 to GCJ02

    Parameters
    -------
    lng : Series or number
        Longitude
    lat : Series or number
        Latitude

    return
    -------
    lng : Series or number
        Longitude (Converted)
    lat : Series or number
        Latitude (Converted)
    """
    try:
        lng = lng.astype(float)
        lat = lat.astype(float)
    except:
        lng = float(lng)
        lat = float(lat)
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * 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 = lat + dlat
    mglng = lng + dlng
    return mglng, mglat


def gcj02towgs84(lng, lat):
    """
    Convert coordinates from GCJ02 to WGS84

    Parameters
    -------
    lng : Series or number
        Longitude
    lat : Series or number
        Latitude

    return
    -------
    lng : Series or number
        Longitude (Converted)
    lat : Series or number
        Latitude (Converted)
    """
    try:
        lng = lng.astype(float)
        lat = lat.astype(float)
    except:
        lng = float(lng)
        lat = float(lat)
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * 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 = lat + dlat
    mglng = lng + dlng
    return lng * 2 - mglng, lat * 2 - mglat

def wgs84tobd09(lon,lat):
    """
    Convert coordinates from WGS84 to BD09 

    Parameters
    -------
    lng : Series or number
        Longitude
    lat : Series or number
        Latitude

    return
    -------
    lng : Series or number
        Longitude (Converted)
    lat : Series or number
        Latitude (Converted)
    """
    try:
        lon = lon.astype(float)
        lat = lat.astype(float)
    except:
        lon = float(lon)
        lat = float(lat)
    lon,lat = wgs84togcj02(lon,lat)
    lon,lat = gcj02tobd09(lon,lat)
    return lon,lat

def bd09towgs84(lon,lat):
    """
    Convert coordinates from BD09 to WGS84

    Parameters
    -------
    lng : Series or number
        Longitude
    lat : Series or number
        Latitude

    return
    -------
    lng : Series or number
        Longitude (Converted)
    lat : Series or number
        Latitude (Converted)
    """
    try:
        lon = lon.astype(float)
        lat = lat.astype(float)
    except:
        lon = float(lon)
        lat = float(lat)
    lon,lat = bd09togcj02(lon,lat)
    lon,lat = gcj02towgs84(lon,lat)
    return lon,lat

def bd09mctobd09(x,y):
    """
    Convert coordinates from BD09MC to BD09

    Parameters
    -------
    x : Series or number
        x coordinates
    y : Series or number
        y coordinates

    return
    -------
    lng : Series or number
        Longitude (Converted)
    lat : Series or number
        Latitude (Converted)
    """
    MCBAND = [12890594.86, 8362377.87, 5591021, 3481989.83, 1678043.12, 0]
    MC2LL = [
      [1.410526172116255e-8, 0.00000898305509648872, -1.9939833816331, 200.9824383106796, -187.2403703815547, 91.6087516669843, -23.38765649603339, 2.57121317296198, -0.03801003308653, 17337981.2],
      [-7.435856389565537e-9, 0.000008983055097726239, -0.78625201886289, 96.32687599759846, -1.85204757529826, -59.36935905485877, 47.40033549296737, -16.50741931063887, 2.28786674699375, 10260144.86],
      [-3.030883460898826e-8, 0.00000898305509983578, 0.30071316287616, 59.74293618442277, 7.357984074871, -25.38371002664745, 13.45380521110908, -3.29883767235584, 0.32710905363475, 6856817.37],
      [-1.981981304930552e-8, 0.000008983055099779535, 0.03278182852591, 40.31678527705744, 0.65659298677277, -4.44255534477492, 0.85341911805263, 0.12923347998204, -0.04625736007561, 4482777.06],
      [3.09191371068437e-9, 0.000008983055096812155, 0.00006995724062, 23.10934304144901, -0.00023663490511, -0.6321817810242, -0.00663494467273, 0.03430082397953, -0.00466043876332, 2555164.4],
      [2.890871144776878e-9, 0.000008983055095805407, -3.068298e-8, 7.47137025468032, -0.00000353937994, -0.02145144861037, -0.00001234426596, 0.00010322952773, -0.00000323890364, 826088.5]
    ]
    y1 = y.iloc[0]
    for cD in range(len(MCBAND)):
        if y1 >= MCBAND[cD]:
            cE = MC2LL[cD]
            break
    cD = cE
    T = cD[0] + cD[1] * np.abs(x);
    cB = np.abs(y) / cD[9]
    cE = cD[2] + cD[3] * cB + cD[4] * cB * cB +\
    cD[5] * cB * cB * cB + cD[6] * cB * cB * cB * cB +\
    cD[7] * cB * cB * cB * cB * cB +\
    cD[8] * cB * cB * cB * cB * cB * cB
    return T,cE

def transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
        0.1 * lng * lat + 0.2 * np.sqrt(np.fabs(lng))
    ret += (20.0 * np.sin(6.0 * lng * pi) + 20.0 *
            np.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * np.sin(lat * pi) + 40.0 *
            np.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * np.sin(lat / 12.0 * pi) + 320 *
            np.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret

def transformlng(lng, lat):
    import numpy as np
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
        0.1 * lng * lat + 0.1 * np.sqrt(np.abs(lng))
    ret += (20.0 * np.sin(6.0 * lng * pi) + 20.0 *
            np.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * np.sin(lng * pi) + 40.0 *
            np.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * np.sin(lng / 12.0 * pi) + 300.0 *
            np.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret

def getdistance(lon1, lat1, lon2, lat2): 
    '''
    Input the origin/destination location in the sequence of [lon1, lat1, lon2, lat2] (in decimal) from DataFrame. The output is the distance (m).

    Parameters
    -------
    lon1 : Series or number
        Start longitude
    lat1 : Series or number
        Start latitude
    lon2 : Series or number
        End longitude
    lat2 : Series or number
        End latitude

    return
    -------
    distance : Series or number
        The distance
    '''
    try:
        lon1 = lon1.astype(float)
        lat1 = lat1.astype(float)
        lon2 = lon2.astype(float)
        lat2 = lat2.astype(float)
    except:
        lon1 = float(lon1)
        lat1 = float(lat1)
        lon2 = float(lon2)
        lat2 = float(lat2)
    lon1, lat1, lon2, lat2 = map(lambda r:r*pi/180, [lon1, lat1, lon2, lat2])
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
    c = 2 * np.arcsin(a**0.5) 
    r = 6371 # 地球平均半径，单位为公里
    return c * r * 1000

def transform_shape(gdf,method):
    '''
    Convert coordinates of all data. The input is the geographic elements’ DataFrame.

    Parameters
    -------
    gdf : GeoDataFrame
        Geographic elements
    method : function
        The coordinate converting function

    return
    -------
    gdf : GeoDataFrame
        The result of converting
    '''
    from shapely.ops import transform
    gdf1 = gdf.copy()
    gdf1['geometry'] = gdf1['geometry'].apply(lambda r:transform(method, r))
    return gdf1
[((99, 12, 99, 26), 'numpy.sin', 'np.sin', ({(99, 19, 99, 25): 'radlat'}, {}), '(radlat)', True, 'import numpy as np\n'), ((101, 16, 101, 30), 'numpy.sqrt', 'np.sqrt', ({(101, 24, 101, 29): 'magic'}, {}), '(magic)', True, 'import numpy as np\n'), ((136, 12, 136, 26), 'numpy.sin', 'np.sin', ({(136, 19, 136, 25): 'radlat'}, {}), '(radlat)', True, 'import numpy as np\n'), ((138, 16, 138, 30), 'numpy.sqrt', 'np.sqrt', ({(138, 24, 138, 29): 'magic'}, {}), '(magic)', True, 'import numpy as np\n'), ((32, 8, 32, 38), 'numpy.sqrt', 'np.sqrt', ({(32, 16, 32, 37): '(lng * lng + lat * lat)'}, {}), '(lng * lng + lat * lat)', True, 'import numpy as np\n'), ((33, 12, 33, 32), 'numpy.arctan2', 'np.arctan2', ({(33, 23, 33, 26): 'lat', (33, 28, 33, 31): 'lng'}, {}), '(lat, lng)', True, 'import numpy as np\n'), ((65, 8, 65, 30), 'numpy.sqrt', 'np.sqrt', ({(65, 16, 65, 29): '(x * x + y * y)'}, {}), '(x * x + y * y)', True, 'import numpy as np\n'), ((66, 12, 66, 28), 'numpy.arctan2', 'np.arctan2', ({(66, 23, 66, 24): 'y', (66, 26, 66, 27): 'x'}, {}), '(y, x)', True, 'import numpy as np\n'), ((67, 17, 67, 30), 'numpy.cos', 'np.cos', ({(67, 24, 67, 29): 'theta'}, {}), '(theta)', True, 'import numpy as np\n'), ((68, 17, 68, 30), 'numpy.sin', 'np.sin', ({(68, 24, 68, 29): 'theta'}, {}), '(theta)', True, 'import numpy as np\n'), ((235, 9, 235, 18), 'numpy.abs', 'np.abs', ({(235, 16, 235, 17): 'y'}, {}), '(y)', True, 'import numpy as np\n'), ((299, 12, 299, 29), 'numpy.arcsin', 'np.arcsin', ({(299, 22, 299, 28): '(a ** 0.5)'}, {}), '(a ** 0.5)', True, 'import numpy as np\n'), ((32, 51, 32, 69), 'numpy.sin', 'np.sin', ({(32, 58, 32, 68): '(lat * x_pi)'}, {}), '(lat * x_pi)', True, 'import numpy as np\n'), ((33, 46, 33, 64), 'numpy.cos', 'np.cos', ({(33, 53, 33, 63): '(lng * x_pi)'}, {}), '(lng * x_pi)', True, 'import numpy as np\n'), ((34, 17, 34, 30), 'numpy.cos', 'np.cos', ({(34, 24, 34, 29): 'theta'}, {}), '(theta)', True, 'import numpy as np\n'), ((35, 17, 35, 30), 'numpy.sin', 'np.sin', ({(35, 24, 35, 29): 'theta'}, {}), '(theta)', True, 'import numpy as np\n'), ((65, 43, 65, 59), 'numpy.sin', 'np.sin', ({(65, 50, 65, 58): '(y * x_pi)'}, {}), '(y * x_pi)', True, 'import numpy as np\n'), ((66, 42, 66, 58), 'numpy.cos', 'np.cos', ({(66, 49, 66, 57): '(x * x_pi)'}, {}), '(x * x_pi)', True, 'import numpy as np\n'), ((234, 24, 234, 33), 'numpy.abs', 'np.abs', ({(234, 31, 234, 32): 'x'}, {}), '(x)', True, 'import numpy as np\n'), ((298, 8, 298, 22), 'numpy.sin', 'np.sin', ({(298, 15, 298, 21): '(dlat / 2)'}, {}), '(dlat / 2)', True, 'import numpy as np\n'), ((321, 55, 321, 75), 'shapely.ops.transform', 'transform', ({(321, 65, 321, 71): 'method', (321, 73, 321, 74): 'r'}, {}), '(method, r)', False, 'from shapely.ops import transform\n'), ((103, 45, 103, 59), 'numpy.cos', 'np.cos', ({(103, 52, 103, 58): 'radlat'}, {}), '(radlat)', True, 'import numpy as np\n'), ((140, 45, 140, 59), 'numpy.cos', 'np.cos', ({(140, 52, 140, 58): 'radlat'}, {}), '(radlat)', True, 'import numpy as np\n'), ((244, 40, 244, 52), 'numpy.fabs', 'np.fabs', ({(244, 48, 244, 51): 'lng'}, {}), '(lng)', True, 'import numpy as np\n'), ((256, 40, 256, 51), 'numpy.abs', 'np.abs', ({(256, 47, 256, 50): 'lng'}, {}), '(lng)', True, 'import numpy as np\n'), ((298, 28, 298, 40), 'numpy.cos', 'np.cos', ({(298, 35, 298, 39): 'lat1'}, {}), '(lat1)', True, 'import numpy as np\n'), ((298, 43, 298, 55), 'numpy.cos', 'np.cos', ({(298, 50, 298, 54): 'lat2'}, {}), '(lat2)', True, 'import numpy as np\n'), ((298, 58, 298, 72), 'numpy.sin', 'np.sin', ({(298, 65, 298, 71): '(dlon / 2)'}, {}), '(dlon / 2)', True, 'import numpy as np\n'), ((245, 19, 245, 41), 'numpy.sin', 'np.sin', ({(245, 26, 245, 40): '(6.0 * lng * pi)'}, {}), '(6.0 * lng * pi)', True, 'import numpy as np\n'), ((246, 12, 246, 34), 'numpy.sin', 'np.sin', ({(246, 19, 246, 33): '(2.0 * lng * pi)'}, {}), '(2.0 * lng * pi)', True, 'import numpy as np\n'), ((247, 19, 247, 35), 'numpy.sin', 'np.sin', ({(247, 26, 247, 34): '(lat * pi)'}, {}), '(lat * pi)', True, 'import numpy as np\n'), ((248, 12, 248, 34), 'numpy.sin', 'np.sin', ({(248, 19, 248, 33): '(lat / 3.0 * pi)'}, {}), '(lat / 3.0 * pi)', True, 'import numpy as np\n'), ((249, 20, 249, 43), 'numpy.sin', 'np.sin', ({(249, 27, 249, 42): '(lat / 12.0 * pi)'}, {}), '(lat / 12.0 * pi)', True, 'import numpy as np\n'), ((250, 12, 250, 35), 'numpy.sin', 'np.sin', ({(250, 19, 250, 34): '(lat * pi / 30.0)'}, {}), '(lat * pi / 30.0)', True, 'import numpy as np\n'), ((257, 19, 257, 41), 'numpy.sin', 'np.sin', ({(257, 26, 257, 40): '(6.0 * lng * pi)'}, {}), '(6.0 * lng * pi)', True, 'import numpy as np\n'), ((258, 12, 258, 34), 'numpy.sin', 'np.sin', ({(258, 19, 258, 33): '(2.0 * lng * pi)'}, {}), '(2.0 * lng * pi)', True, 'import numpy as np\n'), ((259, 19, 259, 35), 'numpy.sin', 'np.sin', ({(259, 26, 259, 34): '(lng * pi)'}, {}), '(lng * pi)', True, 'import numpy as np\n'), ((260, 12, 260, 34), 'numpy.sin', 'np.sin', ({(260, 19, 260, 33): '(lng / 3.0 * pi)'}, {}), '(lng / 3.0 * pi)', True, 'import numpy as np\n'), ((261, 20, 261, 43), 'numpy.sin', 'np.sin', ({(261, 27, 261, 42): '(lng / 12.0 * pi)'}, {}), '(lng / 12.0 * pi)', True, 'import numpy as np\n'), ((262, 12, 262, 35), 'numpy.sin', 'np.sin', ({(262, 19, 262, 34): '(lng / 30.0 * pi)'}, {}), '(lng / 30.0 * pi)', True, 'import numpy as np\n')]