因为项目中经常使用高德的距离和面积计算组件,但是高德并未公布计算逻辑,这就导致项目中数据出问题时不知道该如何去定位。因此花费了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))
-