python计算高德地图距离和面积

因为项目中经常使用高德的距离和面积计算组件，但是高德并未公布计算逻辑，这就导致项目中数据出问题时不知道该如何去定位。因此花费了1天时间把距离计算和面积计算用python语言整理了出来。

距离计算公式:

from math import asin,sqrt,cos,pi

#point_a,point_b是经纬度,格式为[lng,lat]
def st_distance(point_a,point_b):
    #弧度
    radian = pi / 180.0
    #地球半径
    radius = 6378137
    
    f = point_a[1] * radian
    h = point_b[1] * radian
    k = 2 * radius
    d = point_b[0] * radian - point_a[0] * radian
    e = (1 - cos(h - f) + (1 - cos(d)) * cos(f) * cos(h)) / 2
    return k * asin(sqrt(e))
print(round(st_distance([116.434027, 39.941037],[116.461665, 39.941564]),1))

面积计算公式：

from math import cos,pi

def st_area(wkt_list):
    # 弧度
    radian = pi / 180.0
    # 地球半径
    radius = 6378137
    # 地球弧度
    earth_radian = radius * radian
    # 面积
    area = 0
    
    length = len(wkt_list)
    if length < 3 :
        return 0
    for index in range(0,length-1):
        front_point = wkt_list[index]
        rear_point = wkt_list[index+1]
        
        front_sector = front_point[0] * earth_radian * cos(front_point[1] * radian)
        front_line = front_point[1] * earth_radian
        rear_sector = rear_point[0] * earth_radian * cos(rear_point[1] * radian)
        area += (front_sector * rear_point[1] * earth_radian - rear_sector * front_line)
    wkt_rear = wkt_list[length-1]
    wkt_front = wkt_list[0]
    wkt_rear_sector = wkt_rear[0] * earth_radian * cos(wkt_rear[1] * radian)
    wkt_rear_line = wkt_rear[1] * earth_radian
    wkt_front_sector = wkt_front[0] * earth_radian * cos(wkt_front[1] * radian)
    area += wkt_rear_sector * wkt_front[1] * earth_radian - wkt_front_sector * wkt_rear_line
    return 0.5 * abs(area)/1000000

def format_wkt(wkt):
    point_list = wkt.split(',')
    wkt_list=[]
    for point in point_list:
        lng,lat = point.split(' ')
        wkt_list.append([float(lng),float(lat)])
    return wkt_list

#计算面积
wkt='121.543409 31.256203,121.526242 31.248425,121.497403 31.248572,121.48367 31.253268,121.476632 31.267355,121.48161 31.280853,121.492253 31.289508,121.511479 31.292736,121.530362 31.293176,121.544782 31.285401,121.551992 31.270876,121.555253 31.264274,121.543409 31.256203'
print(round(st_area(format_wkt(wkt)),2))