首页 > 代码库 > python 凸包(经纬度) + 面积[近似]

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.0def 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

 

python 凸包(经纬度) + 面积[近似]