python 凸包(经纬度)+面积［近似］

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

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

def getLTDis( A, B ):
    lon1, lat1, lon2, lat2 = map(radians, [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 * 1000.0

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


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 cross( vectorMinus(ans[ m -1 ] , ans [ m - 2 ]), 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 cross( vectorMinus(ans[ m -1 ] , ans [ m - 2 ]), 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 AREA(b)

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 += triangleAre( b[0] , x , y )
    return ans