【多点构成的面积估算】凸包算法python 实现 + 经纬度坐标 python 实现

做用户聚类，有一个feature涉及到面积计算，温故一下之前了解到的凸包算法的流程，实际上，这篇文章基于的原文章，达到的精度，满足我的需求。具体实现上代码参考，因为坐标系不是笛卡尔坐标系

先上笔者参考已有的实现稍作修改后作为静态工具类使用的代码

# the class that seals the algorithm that uses convex hull to calculate the points-distributed-area
from math import *


class cal_area(object):
    @staticmethod
    def deg2rad(deg):
        return deg * (pi / 180)

    @staticmethod
    def cross(A, B):
        return A[0] * B[1] - A[1] * B[0]

    @staticmethod
    def vectorMinus(a, b):
        return ((a[0] - b[0]) * 1000, (a[1] - b[1]) * 1000)

    @staticmethod
    def getLTDis(A, B):
        lon1, lat1, lon2, lat2 = map(cal_area.deg2rad, [A[0], A[1], B[0], B[1]])
        dlon = lon2 - lon1
        dlat = lat2 - lat1
        a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
        c = 2 * asin(sqrt(a))
        r = 6371.393
        # print A,B
        return c * r * 1.0  # 1000.0

    @staticmethod
    def triangleAre(A, B, C):
        x, y, z = cal_area.getLTDis(A, B), cal_area.getLTDis(B, C), cal_area.getLTDis(C, A)
        c = (x + y + z) / 2
        return sqrt((c) * (c - y) * (c - z) * (c - x))

    @staticmethod
    def grahamScanArea(data):
        data.sort(key=lambda x: (x[0], x[1]), reverse=False)
        ans = [0] * (len(data) * 2)
        m = 0
        for item in data:
            top = len(item)
            while (m > 1 and cal_area.cross(cal_area.vectorMinus(ans[m - 1], ans[m - 2]),
                                            cal_area.vectorMinus(item, ans[m - 2])) <= 0): m = m - 1
            ans[m] = item
            m = m + 1
        k = m
        flag = True
        data.reverse()
        for item in data:
            if flag:
                flag = False
                continue
            while (m > k and cal_area.cross(cal_area.vectorMinus(ans[m - 1], ans[m - 2]),
                                            cal_area.vectorMinus(item, ans[m - 2])) <= 0): m = m - 1
            ans[m] = item
            m = m + 1
        m = m - 1
        b = [ans[i] for i in range(0, m)]
        if len(b) < 3: return 0
        # if DEBUG : print b
        return cal_area.AREA(b)

    @staticmethod
    def AREA(b):
        ans = 0.0
        for i in range(len(b)):
            if i == 0 or i + 1 >= len(b): continue
            x, y = b[i], b[i + 1]
            ans += cal_area.triangleAre(b[0], x, y)
        return ans


## 使用方式
# calculate coverage of the business that customer has been to
# form the user_points table
user_visited_area=[]
for user_id in user_business.keys():
    points=[]
    visited_business=user_business[user_id]
    for business_id in visited_business:
        business_