Friday, May 26, 2017

calculate area of a polygon on maps

***** area-polygon.py**********

#~/usr/bin/python

import sys
import math


def getDistance2(co1, co2):
    R = 6371000
    lat1,lng1 =co1.split(",")
    lat2, lng2 =co2.split(",")
    lat1 = float(lat1)
    lng1 = float(lng1)
    lat2 = float(lat2)
    lng2 = float(lng2)
    rlat1 = math.radians(lat1)
    rlat2 = math.radians(lat2)
    rlng1 = math.radians(lng1)
    rlng2 = math.radians(lng2)
    llat = math.radians(lat2 - lat1)
    llng = math.radians(lng2 - lng1)

    d = 2 * R * math.asin(math.sqrt((math.sin((rlat2/2 - rlat1/2))) ** 2 + math.cos(rlat1) * math.cos(rlat2) *(math.sin((rlng2/2 - rlng1/2))) ** 2 ) )

    return d

def slopb(co1, co2):

    lat1,lng1 = co1.split(",")
    lat2, lng2 =co2.split(",")
    lat1 = float(lat1)
    lng1 = float(lng1)
    lat2 = float(lat2)
    lng2 = float(lng2)
    slop = (lng1 - lng2)/(lat1 - lat2)
    b = lng1 -(slop * lat1)

    return slop, b

def ifonline(co1, slop, b):
    # print co1
    lat1,lng1 =co1.split(",")
    lat1 = float(lat1)
    lng1 = float(lng1)
    onLine = False
    if lng1 == (slop * lat1 ) + b:
        onLine = True

    return onLine




def equationxp(a,b):
    x2 = a ** 2
    x = a * b * 2
    c = b ** 2
    eq = [x2,x,c]
    return eq


def beginLineendLine(list1):
    basePoint = list1[0]
    bl = 0
    el = 0
    n = 1
    for i in list1[1:]:
        # print i
        if i == basePoint:
            break
        slop, b = slopb(basePoint, i)
        if ifonline(list1[n+ 1],slop, b):
            bl = n + 1
        else:
            break
        n += 1
    m = len(list1)
    revList1 = list1[:]
    revList1.reverse()
    for i in revList1[1:]:
        if i == basePoint:
            break
        slop, b = slopb(basePoint, i)
        if ifonline(list1[m -1],slop, b):
            el = m -2
        else:
            break
        m -= 1

    return bl, el

def getperpendicular(co1, slop):
    lat1,lng1 =co1.split(",")
    lat1 = float(lat1)
    lng1 = float(lng1)
    newslop = -1/slop
    b = lng1 -(newslop * lat1)

    return newslop, b

def getIntersec(slop1,slop2,b1,b2):

    newx =(b2-b1)/ (slop1 - slop2)
    newy = slop1 * newx + b1

    returnCo = str(newx) + "," + str(newy)
    # print returnCo
    return returnCo


bl,el = beginLineendLine(latlnList)
basePoint = latlnList[0]
if el == 0:
    endpoint = len(latlnList) - 1
else:
    endpoint = el
if bl == 0:
    startpoint = 1
else:
    startpoint = el
n = startpoint
area = 0
for i in latlnList[startpoint:endpoint]:
    # print latlnList[n+1]
    dist1 = getDistance2(basePoint,i)
    slop1, b1 = slopb(basePoint,i)
    slop2, b2 = getperpendicular(latlnList[n+1],slop1)

    intersecPoint = getIntersec(slop1, slop2, b1, b2)
    # print intersecPoint
    dist2 = getDistance2(latlnList[n+1], intersecPoint)
    # print dist1, dist2
    area += (dist1 * dist2) /2


    n+=1

print area
print area/1000000, "km"





**** usages ***
$ python area-polygon.py "11.59607,104.85324" "11.59724,104.86921" "11.58927,104.86984" "11.58871,104.85343" "11.58871,104.85343" "11.58926761448781,104.83902633190155"
2148804.81147
2.14880481147 km



$ python area-polygon.py "11.59607,104.85324" "11.59724,104.86921" "11.58927,104.86984" "11.58871,104.85343" "11.58926761448781,104.83902633190155"
2148804.81147
2.14880481147 km


$ python area-polygon.py "11.59607,104.85324" "11.59724,104.86921" "11.58927,104.86984" "11.58871,104.85343"
1507431.22992
1.50743122992 km

No comments:

Post a Comment