# engine.LingYi.entitys.geo

地理信息计算库
1、求地面两点的水平距离
2、获取三维直线距离
3、获取直线距离
4、角度调整为0-360度以内
5、获取角度差
6、获取point1 指向 point2 的方位角
7、WGS84 坐标系(x, y, z) 转大地坐标系(lat, lon, alt)
8、大地坐标系(lat, lon, alt) 转 WGS84 坐标系(x, y, z)
9、从地球海拔水平上，选一角度出发一定距离后，获取新的点
10、获取地球上两点的差值点
11、获取多边形的中心点
12、判断一个点是否在该多边形内

View Source

      1# -*- coding:utf-8 -*-
      2
      3"""
      4地理信息计算库
      51、求地面两点的水平距离
      62、获取三维直线距离
      73、获取直线距离
      84、角度调整为0-360度以内
      95、获取角度差
     106、获取point1 指向 point2 的方位角
     117、WGS84 坐标系(x, y, z) 转大地坐标系(lat, lon, alt)
     128、大地坐标系(lat, lon, alt) 转 WGS84 坐标系(x, y, z)
     139、从地球海拔水平上，选一角度出发一定距离后，获取新的点
     1410、获取地球上两点的差值点
     1511、获取多边形的中心点
     1612、判断一个点是否在该多边形内
     17"""
     18import math
     19
     20radian2degree = 180 / math.pi
     21degree2radian = math.pi / 180.0
     22NM2KM = 1.852  # 海里转千米
     23
     24# WGS84 坐标系参数：
     25flattening = 1 / 298.257223563  # 扁率
     26EARTH_RADIUS_LONG = 6378137.0   # 地球长半轴 单位：米
     27EARTH_RADIUS_SHORT = EARTH_RADIUS_LONG - EARTH_RADIUS_LONG*flattening
     28# 6356752.3142451793
     29e_pow2 = 1 - (1 - flattening)*(1 - flattening)  # 第一偏心率平方
     30EARTH_RADIUS = 6371137  # 地球平均半径
     31PI = math.pi
     32
     33
     34def get_horizontal_distance(geopoint1, geopoint2):
     35    """
     36    1、求地面两点的水平距离   Haversine公式
     37    :param geopoint1: tuple, (lat, lon), 例：(40.9, 140.0)
     38    :param geopoint2: tuple, (lat, lon), 例：(40.9, 142.0)
     39    :return: float, KM
     40    """
     41    try:
     42        lat1 = geopoint1[0] * degree2radian
     43        lon1 = geopoint1[1] * degree2radian
     44        lat2 = geopoint2[0] * degree2radian
     45        lon2 = geopoint2[1] * degree2radian
     46    except Exception as e:
     47        print("distance error:{},{},{}".format(geopoint1, geopoint2, e))
     48        return 0.0
     49
     50    difference = lat1 - lat2
     51    mdifference = lon1 - lon2
     52    distance = 2 * math.asin(math.sqrt(math.pow(math.sin(difference / 2), 2)
     53                                       + math.cos(lat1) * math.cos(lat2)
     54                                       * math.pow(math.sin(mdifference / 2), 2)))
     55    distance = distance * EARTH_RADIUS / 1000
     56    return distance
     57
     58
     59def get_slant_distance(geopoint1, geopoint2):
     60    """
     61    2、获取三维直线距离, 点高需为海拔高度
     62    :param geopoint1: tuple, (lat, lon, alt), 例：(40.9, 140.0, 560.8)
     63    :param geopoint2: tuple, (lat, lon, alt), 例：(40.9, 142.0, 4560.8)
     64    :return: float, KM
     65    """
     66    hd = get_horizontal_distance(geopoint1, geopoint2)
     67    delta_alt = geopoint1[2] - geopoint2[2]
     68    return get_range(hd, delta_alt)
     69
     70
     71def get_range(range_km, delta_alt):
     72    """
     73    3、获取直线距离
     74    :param range_km: float, 水平距离，KM
     75    :param delta_alt: float, 垂直距离，m
     76    :return: float, KM
     77    """
     78    range_km *= 1000.0
     79    return math.sqrt((range_km * range_km + delta_alt * delta_alt)) / 1000.0
     80
     81
     82def normal_angle(angle):
     83    """
     84    4、角度调整为0-360度以内
     85    :param angle: float, 角度
     86    :return: float
     87    """
     88    if 0 <= angle < 360:
     89        return angle
     90    else:
     91        return angle % 360
     92
     93
     94def get_angle_span(angle1: float, angle2: float) -> float:
     95    """
     96    5、获取角度差
     97    :param angle1:
     98    :param angle2:
     99    :return:
    100    """
    101    angle_span = angle1-angle2 if angle1 >= angle2 else angle2-angle1
    102    if angle_span <= 180:
    103        return angle_span
    104    else:
    105        return 360-angle_span
    106
    107
    108def get_azimuth(geopoint1, geopoint2):
    109    """
    110    6、获取point1 指向 point2 的方位角
    111    :param geopoint1: tuple, (lat, lon), 例：(40.9, 140.0)
    112    :param geopoint2: tuple, (lat, lon), 例：(40.9, 142.0)
    113    :return: 角度 0-360, 正北：0， 正东:90, 顺时针旋转，正西：270
    114    """
    115    lat1 = geopoint1[0] * degree2radian
    116    lon1 = geopoint1[1] * degree2radian
    117    lat2 = geopoint2[0] * degree2radian
    118    lon2 = geopoint2[1] * degree2radian
    119    azimuth = 180 * math.atan2(math.sin(lon2 - lon1), math.tan(lat2) * math.cos(lat1) - math.sin(lat1) * math.cos(lon2 - lon1)) / PI
    120    return normal_angle(azimuth)
    121
    122
    123def worldxyz2lla(x, y, z):
    124    """
    125    7、WGS84 坐标系(x, y, z) 转 大地坐标系(lat, lon, alt)
    126    :param x: float, m
    127    :param y: float, m
    128    :param z: float, m
    129    :return: tuple, (lat, lon, alt), 单位：(degree,小数, 小数, m)
    130    """
    131    alt = 0.0
    132    if not x:
    133        if not y:
    134            if z > 0:
    135                lon = 0.0
    136                lat = 90.0
    137                alt = z - EARTH_RADIUS_SHORT
    138            else:
    139                lon = 0.0
    140                lat = -90.0
    141                alt = -z - EARTH_RADIUS_SHORT
    142            return lat, lon, alt
    143        elif y > 0:
    144            lon = PI / 2
    145        else:
    146            lon = - PI/2
    147    else:
    148        lon = math.atan(math.fabs(y / x))
    149    if not z:
    150        # 在赤道平面，纬度为0
    151        lat = 0.0
    152        alt = math.sqrt(x * x + y * y) - EARTH_RADIUS_LONG
    153    else:
    154        lat = math.atan(math.fabs(z / math.sqrt(x * x + y * y)))
    155        for i in range(5):
    156            N = EARTH_RADIUS_LONG / math.sqrt(1.0 - e_pow2 * math.sin(lat) * math.sin(lat))
    157            alt = z / math.sin(lat) - N * (1.0 - e_pow2)
    158            lat = math.atan(z * (N + alt) / (math.sqrt(x * x + y * y) * (N * (1.0 - e_pow2) + alt)))
    159    if x < 0 and y > 0:
    160        lon = PI - lon
    161    elif x > 0 and y < 0:
    162        lon = -1.0 * lon
    163    elif x < 0 and y < 0:
    164        lon = lon - PI
    165    lat = lat * radian2degree
    166    lon = lon * radian2degree
    167    return lat, lon, alt
    168
    169
    170def lla2worldxyz(lat, lon, alt):
    171    """
    172    8、大地坐标系(lat, lon, alt) 转 WGS84 坐标系(x, y, z)
    173    :param lat: float, 纬度
    174    :param lon: float, 经度
    175    :param alt: float, 高度，m
    176    :return: x, y, z
    177    """
    178    N = EARTH_RADIUS_LONG / math.sqrt(1 - e_pow2 * math.sin(lat * degree2radian) * math.sin(lat * degree2radian))
    179    x = (N + alt) * math.cos(lat * degree2radian) * math.cos(lon * degree2radian)
    180    y = (N + alt) * math.cos(lat * degree2radian) * math.sin(lon * degree2radian)
    181    z = (N * (1 - e_pow2) + alt) * math.sin(lat * degree2radian)
    182    return x, y, z
    183
    184
    185def get_geopoint_from_distance(geo_point, azimuth, distance_m):
    186    """
    187    9、从地球海拔水平上，选一角度出发一定距离后，获取新的点. 距离越远，精度越差
    188    :param geo_point: tuple, (float, float), (纬度, 经度)
    189    :param azimuth: float, 角度，0-360， 正北为0， 顺时针旋转360度
    190    :param distance_m: float, 距离，单位：m
    191    :return: tuple, (lat, lon)
    192    """
    193    lat = geo_point[0]
    194    lon = geo_point[1]
    195    a = EARTH_RADIUS_LONG
    196    b = EARTH_RADIUS_SHORT
    197    alpha1 = azimuth * PI / 180
    198    sinAlpha1 = math.sin(alpha1)
    199    cosAlpha1 = math.cos(alpha1)
    200
    201    tanU1 = (1 - flattening) * math.tan(lat*PI/180)
    202    cosU1 = 1 / math.sqrt((1 + tanU1*tanU1))
    203    sinU1 = tanU1 * cosU1
    204    sigma1 = math.atan2(tanU1, cosAlpha1)
    205    sinAlpha = cosU1 * sinAlpha1
    206    cosSqAlpha = 1 - sinAlpha * sinAlpha
    207    uSq = cosSqAlpha * (a*a - b*b) / (b*b)
    208    A = 1 + uSq/16384*(4096 + uSq*(-768 + uSq*(320 - 175*uSq)))
    209    B = uSq/1024*(256 + uSq*(-128 + uSq*(74 - 47*uSq)))
    210
    211    sigma = distance_m / (b * A)
    212    sigmaP = 2 * PI
    213    sinSigma = 0
    214    cosSigma = 0
    215    cos2SigmaM = 0
    216    for i in range(8):
    217        if math.fabs(sigma - sigmaP) < 1e-12:
    218            break
    219        cos2SigmaM = math.cos(2*sigma1 + sigma)
    220        sinSigma = math.sin(sigma)
    221        cosSigma = math.cos(sigma)
    222        deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma*(-1 + 2*cos2SigmaM*cos2SigmaM)
    223                                                                        - B/6*cos2SigmaM*(-3 + 4 * sinSigma*sinSigma) * (-3 + 4*cos2SigmaM*cos2SigmaM)))
    224        sigmaP = sigma
    225        sigma = distance_m / (b*A) + deltaSigma
    226
    227    tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1
    228    lat2 = math.atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1,
    229                      (1 - flattening) * math.sqrt(sinAlpha*sinAlpha + tmp*tmp))
    230    lon_span = math.atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1)
    231    C = flattening / 16 * cosSqAlpha*(4 + flattening*(4 - 3*cosSqAlpha))
    232    lon_diff = lon_span - (1 - C)*flattening*sinAlpha*(sigma + C*sinSigma*(cos2SigmaM + C*cosSigma*(-1 + 2*cos2SigmaM*cos2SigmaM)))
    233    return lat2*180/PI, lon+lon_diff*180/PI
    234
    235
    236def get_two_point_dvalue(point_start, point_to, d_rate):
    237    """
    238    10、获取地球上两点的差值点
    239    :param point_start: tuple, (lat, lon), 例：(40.9, 140.0)
    240    :param point_to: tuple, (lat, lon), 例：(40.9, 142.0)
    241    :param d_rate: float, 差值比例, 0: 返回point_start, 1: 返回point_to, 2: 往到达点2倍距离点
    242    :return: tuple, (lat, lon)
    243    """
    244    if d_rate < 0:
    245        d_rate *= -1
    246        point_change = point_to
    247        point_to = point_start
    248        point_start = point_change
    249    distance = get_horizontal_distance(point_start, point_to) * d_rate * 1000
    250    azimuth = get_azimuth(point_start, point_to)
    251    return get_geopoint_from_distance(point_start, azimuth, distance)
    252
    253
    254def get_polygon_center(polygon):
    255    """
    256    11、获取多边形的中心点
    257    :param polygon: list, 参考点列表,例:[(40, 39.0)]，其中纬度值在前，经度值在后
    258    :return: tuple, (lat, lon), 例：(40.9, 140.0)
    259    """
    260    if polygon is None:
    261        return None
    262    n_count = len(polygon)
    263    if n_count == 0:
    264        return None
    265    lon_r = 0
    266    lat_r = 0
    267    for (lat, lon) in polygon:
    268        lat_r += lat
    269        lon_r += lon
    270    lon_r = lon_r/n_count
    271    lat_r = lat_r/n_count
    272    return lat_r, lon_r
    273
    274
    275def IsPtInPoly(aLat, aLon, pointList):
    276    """
    277    12、判断一个点是否在该多边形内
    278    :param aLat: double 纬度
    279    :param aLon: double 经度
    280    :param pointList: list [(lat, lon)...] 多边形点的顺序需根据顺时针或逆时针，不能乱
    281    """
    282    iSum = 0
    283    iCount = len(pointList)
    284    if iCount < 3:
    285        return False
    286    for i in range(iCount):
    287        pLon1 = pointList[i][1]
    288        pLat1 = pointList[i][0]
    289        if i == iCount - 1:
    290            pLon2 = pointList[0][1]
    291            pLat2 = pointList[0][0]
    292        else:
    293            pLon2 = pointList[i + 1][1]
    294            pLat2 = pointList[i + 1][0]
    295
    296        if ((aLat >= pLat1) and (aLat < pLat2)) or ((aLat>=pLat2) and (aLat < pLat1)):
    297            if abs(pLat1 - pLat2) > 0:
    298                pLon = pLon1 - ((pLon1 - pLon2) * (pLat1 - aLat)) / (pLat1 - pLat2)
    299                if pLon < aLon:
    300                    iSum += 1
    301    return iSum % 2 != 0
    302
    303
    304if __name__ == "__main__":
    305    # 点1：(lat, lon, alt(单位:m))
    306    point1 = (0.0, 0.0, 0)
    307    # 点2：(lat, lon, alt(单位:m))
    308    point2 = (30.0, 42.6, 3000)
    309    point_list = [[31.2, 41.8, 375],
    310                    [31.2, 41.9, 958],
    311                    [31.2, 42.3, 306],
    312                    [30.8, 42.4, 1255],
    313                    [30.2, 42.7, 8412],
    314                    [30.1, 42.9, 467],
    315                    [30.0, 43.2, 412],
    316                    [29.6, 43.9, 10518],
    317                    [29.2, 44.4, 528],
    318                    [29.0, 45.9, 8563],
    319                    [29.0, 46.5, 2536]]
    320
    321    # 求水平距离
    322    d1 = get_horizontal_distance(point2, point_list[0])
    323    print('{},{} horizontal distance:{}KM'.format(point2, point_list[0], d1))
    324
    325    # 求斜面距离
    326    d1 = get_slant_distance(point2, point_list[0])
    327    print('{},{} distance:{}KM'.format(point2, point_list[0], d1))
    328
    329    # 求点1往东10km的点
    330    east_point = get_geopoint_from_distance(point1, 90, 10000)
    331    print("{} east 10km is {}".format(point1, east_point))
    332
    333    # 求点1往西北200km的点
    334    wn_point = get_geopoint_from_distance(point1, 315, 200000)
    335    print("{} northwest 200km is {}".format(point1, wn_point))
    336
    337    # 点1到点2指向角度 正北为0，顺时针旋转到360
    338    azimuth = get_azimuth(point1, point2)
    339    print("{} to {}, azimuth degree {} degree".format(point1, point2, azimuth))
    340
    341    # 点2到点3分成3份，求中间两个点
    342    th_point1 = get_two_point_dvalue(point2, point_list[0], 1/3)
    343    th_point2 = get_two_point_dvalue(point2, point_list[0], 2/3)
    344    print("{} and {} Divide into three parts is {} and {}".format(point2, point_list[0], th_point1, th_point2))

radian2degree = 57.29577951308232
degree2radian = 0.017453292519943295
NM2KM = 1.852
flattening = 0.0033528106647474805
EARTH_RADIUS_LONG = 6378137.0
EARTH_RADIUS_SHORT = 6356752.314245179
e_pow2 = 0.00669437999014133
EARTH_RADIUS = 6371137
PI = 3.141592653589793
def get_horizontal_distance(geopoint1, geopoint2): View Source

    35def get_horizontal_distance(geopoint1, geopoint2):
    36    """
    37    1、求地面两点的水平距离   Haversine公式
    38    :param geopoint1: tuple, (lat, lon), 例：(40.9, 140.0)
    39    :param geopoint2: tuple, (lat, lon), 例：(40.9, 142.0)
    40    :return: float, KM
    41    """
    42    try:
    43        lat1 = geopoint1[0] * degree2radian
    44        lon1 = geopoint1[1] * degree2radian
    45        lat2 = geopoint2[0] * degree2radian
    46        lon2 = geopoint2[1] * degree2radian
    47    except Exception as e:
    48        print("distance error:{},{},{}".format(geopoint1, geopoint2, e))
    49        return 0.0
    50
    51    difference = lat1 - lat2
    52    mdifference = lon1 - lon2
    53    distance = 2 * math.asin(math.sqrt(math.pow(math.sin(difference / 2), 2)
    54                                       + math.cos(lat1) * math.cos(lat2)
    55                                       * math.pow(math.sin(mdifference / 2), 2)))
    56    distance = distance * EARTH_RADIUS / 1000
    57    return distance

1、求地面两点的水平距离 Haversine公式 :param geopoint1: tuple, (lat, lon), 例：(40.9, 140.0) :param geopoint2: tuple, (lat, lon), 例：(40.9, 142.0) :return: float, KM

def get_slant_distance(geopoint1, geopoint2): View Source

    60def get_slant_distance(geopoint1, geopoint2):
    61    """
    62    2、获取三维直线距离, 点高需为海拔高度
    63    :param geopoint1: tuple, (lat, lon, alt), 例：(40.9, 140.0, 560.8)
    64    :param geopoint2: tuple, (lat, lon, alt), 例：(40.9, 142.0, 4560.8)
    65    :return: float, KM
    66    """
    67    hd = get_horizontal_distance(geopoint1, geopoint2)
    68    delta_alt = geopoint1[2] - geopoint2[2]
    69    return get_range(hd, delta_alt)

2、获取三维直线距离, 点高需为海拔高度 :param geopoint1: tuple, (lat, lon, alt), 例：(40.9, 140.0, 560.8) :param geopoint2: tuple, (lat, lon, alt), 例：(40.9, 142.0, 4560.8) :return: float, KM

def get_range(range_km, delta_alt): View Source

    72def get_range(range_km, delta_alt):
    73    """
    74    3、获取直线距离
    75    :param range_km: float, 水平距离，KM
    76    :param delta_alt: float, 垂直距离，m
    77    :return: float, KM
    78    """
    79    range_km *= 1000.0
    80    return math.sqrt((range_km * range_km + delta_alt * delta_alt)) / 1000.0

3、获取直线距离 :param range_km: float, 水平距离，KM :param delta_alt: float, 垂直距离，m :return: float, KM

def normal_angle(angle): View Source

    83def normal_angle(angle):
    84    """
    85    4、角度调整为0-360度以内
    86    :param angle: float, 角度
    87    :return: float
    88    """
    89    if 0 <= angle < 360:
    90        return angle
    91    else:
    92        return angle % 360

4、角度调整为0-360度以内 :param angle: float, 角度 :return: float

def get_angle_span(angle1: float, angle2: float) -> float: View Source

     95def get_angle_span(angle1: float, angle2: float) -> float:
     96    """
     97    5、获取角度差
     98    :param angle1:
     99    :param angle2:
    100    :return:
    101    """
    102    angle_span = angle1-angle2 if angle1 >= angle2 else angle2-angle1
    103    if angle_span <= 180:
    104        return angle_span
    105    else:
    106        return 360-angle_span

5、获取角度差 :param angle1: :param angle2: :return:

def get_azimuth(geopoint1, geopoint2): View Source

    109def get_azimuth(geopoint1, geopoint2):
    110    """
    111    6、获取point1 指向 point2 的方位角
    112    :param geopoint1: tuple, (lat, lon), 例：(40.9, 140.0)
    113    :param geopoint2: tuple, (lat, lon), 例：(40.9, 142.0)
    114    :return: 角度 0-360, 正北：0， 正东:90, 顺时针旋转，正西：270
    115    """
    116    lat1 = geopoint1[0] * degree2radian
    117    lon1 = geopoint1[1] * degree2radian
    118    lat2 = geopoint2[0] * degree2radian
    119    lon2 = geopoint2[1] * degree2radian
    120    azimuth = 180 * math.atan2(math.sin(lon2 - lon1), math.tan(lat2) * math.cos(lat1) - math.sin(lat1) * math.cos(lon2 - lon1)) / PI
    121    return normal_angle(azimuth)

6、获取point1 指向 point2 的方位角 :param geopoint1: tuple, (lat, lon), 例：(40.9, 140.0) :param geopoint2: tuple, (lat, lon), 例：(40.9, 142.0) :return: 角度 0-360, 正北：0， 正东:90, 顺时针旋转，正西：270

def worldxyz2lla(x, y, z): View Source

    124def worldxyz2lla(x, y, z):
    125    """
    126    7、WGS84 坐标系(x, y, z) 转 大地坐标系(lat, lon, alt)
    127    :param x: float, m
    128    :param y: float, m
    129    :param z: float, m
    130    :return: tuple, (lat, lon, alt), 单位：(degree,小数, 小数, m)
    131    """
    132    alt = 0.0
    133    if not x:
    134        if not y:
    135            if z > 0:
    136                lon = 0.0
    137                lat = 90.0
    138                alt = z - EARTH_RADIUS_SHORT
    139            else:
    140                lon = 0.0
    141                lat = -90.0
    142                alt = -z - EARTH_RADIUS_SHORT
    143            return lat, lon, alt
    144        elif y > 0:
    145            lon = PI / 2
    146        else:
    147            lon = - PI/2
    148    else:
    149        lon = math.atan(math.fabs(y / x))
    150    if not z:
    151        # 在赤道平面，纬度为0
    152        lat = 0.0
    153        alt = math.sqrt(x * x + y * y) - EARTH_RADIUS_LONG
    154    else:
    155        lat = math.atan(math.fabs(z / math.sqrt(x * x + y * y)))
    156        for i in range(5):
    157            N = EARTH_RADIUS_LONG / math.sqrt(1.0 - e_pow2 * math.sin(lat) * math.sin(lat))
    158            alt = z / math.sin(lat) - N * (1.0 - e_pow2)
    159            lat = math.atan(z * (N + alt) / (math.sqrt(x * x + y * y) * (N * (1.0 - e_pow2) + alt)))
    160    if x < 0 and y > 0:
    161        lon = PI - lon
    162    elif x > 0 and y < 0:
    163        lon = -1.0 * lon
    164    elif x < 0 and y < 0:
    165        lon = lon - PI
    166    lat = lat * radian2degree
    167    lon = lon * radian2degree
    168    return lat, lon, alt

7、WGS84 坐标系(x, y, z) 转 大地坐标系(lat, lon, alt) :param x: float, m :param y: float, m :param z: float, m :return: tuple, (lat, lon, alt), 单位：(degree,小数, 小数, m)

def lla2worldxyz(lat, lon, alt): View Source

    171def lla2worldxyz(lat, lon, alt):
    172    """
    173    8、大地坐标系(lat, lon, alt) 转 WGS84 坐标系(x, y, z)
    174    :param lat: float, 纬度
    175    :param lon: float, 经度
    176    :param alt: float, 高度，m
    177    :return: x, y, z
    178    """
    179    N = EARTH_RADIUS_LONG / math.sqrt(1 - e_pow2 * math.sin(lat * degree2radian) * math.sin(lat * degree2radian))
    180    x = (N + alt) * math.cos(lat * degree2radian) * math.cos(lon * degree2radian)
    181    y = (N + alt) * math.cos(lat * degree2radian) * math.sin(lon * degree2radian)
    182    z = (N * (1 - e_pow2) + alt) * math.sin(lat * degree2radian)
    183    return x, y, z

8、大地坐标系(lat, lon, alt) 转 WGS84 坐标系(x, y, z) :param lat: float, 纬度 :param lon: float, 经度 :param alt: float, 高度，m :return: x, y, z

def get_geopoint_from_distance(geo_point, azimuth, distance_m): View Source

    186def get_geopoint_from_distance(geo_point, azimuth, distance_m):
    187    """
    188    9、从地球海拔水平上，选一角度出发一定距离后，获取新的点. 距离越远，精度越差
    189    :param geo_point: tuple, (float, float), (纬度, 经度)
    190    :param azimuth: float, 角度，0-360， 正北为0， 顺时针旋转360度
    191    :param distance_m: float, 距离，单位：m
    192    :return: tuple, (lat, lon)
    193    """
    194    lat = geo_point[0]
    195    lon = geo_point[1]
    196    a = EARTH_RADIUS_LONG
    197    b = EARTH_RADIUS_SHORT
    198    alpha1 = azimuth * PI / 180
    199    sinAlpha1 = math.sin(alpha1)
    200    cosAlpha1 = math.cos(alpha1)
    201
    202    tanU1 = (1 - flattening) * math.tan(lat*PI/180)
    203    cosU1 = 1 / math.sqrt((1 + tanU1*tanU1))
    204    sinU1 = tanU1 * cosU1
    205    sigma1 = math.atan2(tanU1, cosAlpha1)
    206    sinAlpha = cosU1 * sinAlpha1
    207    cosSqAlpha = 1 - sinAlpha * sinAlpha
    208    uSq = cosSqAlpha * (a*a - b*b) / (b*b)
    209    A = 1 + uSq/16384*(4096 + uSq*(-768 + uSq*(320 - 175*uSq)))
    210    B = uSq/1024*(256 + uSq*(-128 + uSq*(74 - 47*uSq)))
    211
    212    sigma = distance_m / (b * A)
    213    sigmaP = 2 * PI
    214    sinSigma = 0
    215    cosSigma = 0
    216    cos2SigmaM = 0
    217    for i in range(8):
    218        if math.fabs(sigma - sigmaP) < 1e-12:
    219            break
    220        cos2SigmaM = math.cos(2*sigma1 + sigma)
    221        sinSigma = math.sin(sigma)
    222        cosSigma = math.cos(sigma)
    223        deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma*(-1 + 2*cos2SigmaM*cos2SigmaM)
    224                                                                        - B/6*cos2SigmaM*(-3 + 4 * sinSigma*sinSigma) * (-3 + 4*cos2SigmaM*cos2SigmaM)))
    225        sigmaP = sigma
    226        sigma = distance_m / (b*A) + deltaSigma
    227
    228    tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1
    229    lat2 = math.atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1,
    230                      (1 - flattening) * math.sqrt(sinAlpha*sinAlpha + tmp*tmp))
    231    lon_span = math.atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1)
    232    C = flattening / 16 * cosSqAlpha*(4 + flattening*(4 - 3*cosSqAlpha))
    233    lon_diff = lon_span - (1 - C)*flattening*sinAlpha*(sigma + C*sinSigma*(cos2SigmaM + C*cosSigma*(-1 + 2*cos2SigmaM*cos2SigmaM)))
    234    return lat2*180/PI, lon+lon_diff*180/PI

9、从地球海拔水平上，选一角度出发一定距离后，获取新的点. 距离越远，精度越差 :param geo_point: tuple, (float, float), (纬度, 经度) :param azimuth: float, 角度，0-360， 正北为0， 顺时针旋转360度 :param distance_m: float, 距离，单位：m :return: tuple, (lat, lon)

def get_two_point_dvalue(point_start, point_to, d_rate): View Source

    237def get_two_point_dvalue(point_start, point_to, d_rate):
    238    """
    239    10、获取地球上两点的差值点
    240    :param point_start: tuple, (lat, lon), 例：(40.9, 140.0)
    241    :param point_to: tuple, (lat, lon), 例：(40.9, 142.0)
    242    :param d_rate: float, 差值比例, 0: 返回point_start, 1: 返回point_to, 2: 往到达点2倍距离点
    243    :return: tuple, (lat, lon)
    244    """
    245    if d_rate < 0:
    246        d_rate *= -1
    247        point_change = point_to
    248        point_to = point_start
    249        point_start = point_change
    250    distance = get_horizontal_distance(point_start, point_to) * d_rate * 1000
    251    azimuth = get_azimuth(point_start, point_to)
    252    return get_geopoint_from_distance(point_start, azimuth, distance)

10、获取地球上两点的差值点 :param point_start: tuple, (lat, lon), 例：(40.9, 140.0) :param point_to: tuple, (lat, lon), 例：(40.9, 142.0) :param d_rate: float, 差值比例, 0: 返回point_start, 1: 返回point_to, 2: 往到达点2倍距离点 :return: tuple, (lat, lon)

def get_polygon_center(polygon): View Source

    255def get_polygon_center(polygon):
    256    """
    257    11、获取多边形的中心点
    258    :param polygon: list, 参考点列表,例:[(40, 39.0)]，其中纬度值在前，经度值在后
    259    :return: tuple, (lat, lon), 例：(40.9, 140.0)
    260    """
    261    if polygon is None:
    262        return None
    263    n_count = len(polygon)
    264    if n_count == 0:
    265        return None
    266    lon_r = 0
    267    lat_r = 0
    268    for (lat, lon) in polygon:
    269        lat_r += lat
    270        lon_r += lon
    271    lon_r = lon_r/n_count
    272    lat_r = lat_r/n_count
    273    return lat_r, lon_r

11、获取多边形的中心点 :param polygon: list, 参考点列表,例:[(40, 39.0)]，其中纬度值在前，经度值在后 :return: tuple, (lat, lon), 例：(40.9, 140.0)

def IsPtInPoly(aLat, aLon, pointList): View Source

    276def IsPtInPoly(aLat, aLon, pointList):
    277    """
    278    12、判断一个点是否在该多边形内
    279    :param aLat: double 纬度
    280    :param aLon: double 经度
    281    :param pointList: list [(lat, lon)...] 多边形点的顺序需根据顺时针或逆时针，不能乱
    282    """
    283    iSum = 0
    284    iCount = len(pointList)
    285    if iCount < 3:
    286        return False
    287    for i in range(iCount):
    288        pLon1 = pointList[i][1]
    289        pLat1 = pointList[i][0]
    290        if i == iCount - 1:
    291            pLon2 = pointList[0][1]
    292            pLat2 = pointList[0][0]
    293        else:
    294            pLon2 = pointList[i + 1][1]
    295            pLat2 = pointList[i + 1][0]
    296
    297        if ((aLat >= pLat1) and (aLat < pLat2)) or ((aLat>=pLat2) and (aLat < pLat1)):
    298            if abs(pLat1 - pLat2) > 0:
    299                pLon = pLon1 - ((pLon1 - pLon2) * (pLat1 - aLat)) / (pLat1 - pLat2)
    300                if pLon < aLon:
    301                    iSum += 1
    302    return iSum % 2 != 0

12、判断一个点是否在该多边形内 :param aLat: double 纬度 :param aLon: double 经度 :param pointList: list [(lat, lon)...] 多边形点的顺序需根据顺时针或逆时针，不能乱