Python求凸包及多边形面积教程

import sys
import math
import time
import random

#获取基准点的下标,基准点是p[k]
def get_leftbottompoint(p):
 k = 0
 for i in range(1, len(p)):
  if p[i][1] < p[k][1] or (p[i][1] == p[k][1] and p[i][0] < p[k][0]):
   k = i
 return k

#叉乘计算方法
def multiply(p1, p2, p0):
 return (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p2[0] - p0[0]) * (p1[1] - p0[1])

#获取极角，通过求反正切得出，考虑pi/2的情况
def get_arc(p1, p0):
 # 兼容sort_points_tan的考虑
 if (p1[0] - p0[0]) == 0:
  if ((p1[1] - p0[1])) == 0:
   return -1;
  else:
   return math.pi / 2
 tan = float((p1[1] - p0[1])) / float((p1[0] - p0[0]))
 arc = math.atan(tan)
 if arc >= 0:
  return arc
 else:
  return math.pi + arc

#对极角进行排序,排序结果list不包含基准点
def sort_points_tan(p, pk):
 p2 = []
 for i in range(0, len(p)):
  p2.append({"index": i, "arc": get_arc(p[i], pk)})
 #print('排序前:',p2)
 p2.sort(key=lambda k: (k.get('arc')))
 #print('排序后:',p2)
 p_out = []
 for i in range(0, len(p2)):
  p_out.append(p[p2[i]["index"]])
 return p_out

def convex_hull(p):
 p=list(set(p))
 #print('全部点:',p)
 k = get_leftbottompoint(p)
 pk = p[k]
 p.remove(p[k])
 #print('排序前去除基准点的所有点:',p,'基准点:',pk)

 p_sort = sort_points_tan(p, pk) #按与基准点连线和x轴正向的夹角排序后的点坐标
 #print('其余点与基准点夹角排序:',p_sort)
 p_result = [pk,p_sort[0]]

 top = 2
 for i in range(1, len(p_sort)):
  #####################################
  #叉乘为正,向前递归删点;叉乘为负,序列追加新点
  while(multiply(p_result[-2], p_sort[i],p_result[-1]) > 0):
   p_result.pop()
  p_result.append(p_sort[i]) 
 return p_result#测试

if __name__ == '__main__':
 pass
 test_data = [(220, -100), (0,0), (-40, -170), (240, 50), (-160, 150), (-210, -150)]
 print(test_data)

 result = convex_hull(test_data)
 print(result)
 t=0

import matplotlib.pyplot as plt
x1=[]
y1=[]
for i in range(len(test_data)):
 ri=test_data[i]
 #print(ri)
 x1.append(ri[0])
 y1.append(ri[1])

plt.plot(x1,y1,linestyle=' ',marker='.')


xx=[]
yy=[]
for i in range(len(result)):
 ri=result[i]
 #print(ri)
 xx.append(ri[0])
 yy.append(ri[1])

plt.plot(xx,yy,linestyle=' ',marker='*')

def GetAreaOfPolyGonbyVector(points):
 # 基于向量叉乘计算多边形面积
 area = 0
 if(len(points)<3):
  raise Exception("error")

 for i in range(0,len(points)-1):
  p1 = points[i]
  p2 = points[i + 1]

  triArea = (p1[0]*p2[1] - p2[0]*p1[1])/2
  #print(triArea)
  area += triArea

 fn=(points[-1][0]*points[0][1]-points[0][0]*points[-1][1])/2
 #print(fn)
 return abs(area+fn)

points = []
x = [1,3,2]
y = [1,2,2] 
#[(1,1),(3,1),(5,3),(3,5),(1,3)] 
# x=[1,3,5,3,1]
# y=[1,1,3,5,3]
for index in range(len(x)):
 points.append((x[index],y[index]))
area = GetAreaOfPolyGonbyVector(points)
print(area)
#print(math.ceil(area))

from geopy.distance import vincenty
import math
def HeronGetAreaOfPolyGonbyVector(points):
 # 基于海伦公式计算多边形面积
 area = 0
 if(len(points)<3):
  raise Exception("error")

 pb=((points[-1][0]+points[0][0])/2,(points[-1][1]+points[0][1])/2) #基准点选为第一个点和最后一个点连线边上的中点

 for i in range(0,len(points)-1):
  p1 = points[i]
  p2 = points[i + 1]

  db1 = vincenty(pb,p1).meters #根据维度转化成经纬度距离
  d12 = vincenty(p1,p2).meters
  d2b = vincenty(p2,pb).meters
  #print(db1,d12,d2b)

  hc = (db1+d12+d2b)/2 #db1是基准点和p1的距离，d12是p1和p2的距离，d2b是p2和基准点距离
  #print(hc, hc-db1, hc-d12, hc-d2b)
  triArea = math.sqrt(hc*(hc-db1)*(hc-d12)*(hc-d2b)) 
  #print(triArea)
  area += triArea

 return area


points = []
x = [1,3,2]
y = [1,2,2] 
#[(1,1),(3,1),(5,3),(3,5),(1,3)] 
# x=[1,3,5,3,1]
# y=[1,1,3,5,3]
for index in range(len(x)):
 points.append((x[index],y[index]))

area = HeronGetAreaOfPolyGonbyVector(points)
print(area)
#print(math.ceil(area))

p2=[{"index":0, "arc":get_arc(p[0],p[k])}，
 {"index":1, "arc":get_arc(p[1],p[k])}，
 ···
 {"index":k-1, "arc":get_arc(p[k-1],p[k])}，
 {"index":k+1, "arc":get_arc(p[k+1],p[k])}，
 ···
 {"index":n, "arc":get_arc(p[n],p[k])}]
#get_arc(p[0],p[k])即获得p[0]点与基准点p[k]连线的极角(与x轴正向夹角)
#根据p2的“arc”键的值从小到大排序，最后输出按该角度值排序对应顺序的各个点