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