Wednesday, May 24, 2017

create a circle from a latitue/longitude point and distance

**** circle-from-latitude.py ****
#~/usr/bin/python

import sys
import math

pa = sys.argv[1].strip()
# pb = sys.argv[2].strip()
dst = float(sys.argv[2].strip())

pal = pa.split(",")
palt = float(pal[0])
paln = float(pal[1])

def getBpoint(co1, dst):
    R = 6371000
    lat1, lng1 =co1.split(",")
    # lat2, lng2 =co2.split(",")
    lat1 = float(lat1)
    lng1 = float(lng1)
    # lat2 = float(lat2)
    # lng2 = float(lng2)
    radianLat1 = math.radians(lat1)
    radianLat2 = math.radians(lat1)
    radianLng1 = math.radians(lng1)
    part1 = math.sin(dst/(2*R)) ** 2
    part2 = (math.sin((radianLat2/2 - radianLat1/2))) ** 2
    part3 = math.cos(radianLat1) * math.cos(radianLat2)
    part3I = (part1- part2)/part3
    part4 = math.asin(math.sqrt(part3I))
    # part3 = (math.sin((radianLng2/2 - radianLng1/2))) ** 2
    radianLng2 = (part4 + (radianLng1/2)) * 2

    degreeLat2 = math.degrees(radianLat2)
    degreeLng2 = math.degrees(radianLng2)
    return lat1, degreeLng2


def getXpoint(co1, nlat, dst):
    R = 6371000
    lat1, lng1 =co1.split(",")
    # lat2, lng2 =co2.split(",")
    lat1 = float(lat1)
    lng1 = float(lng1)
    # lat2 = float(lat2)
    # lng2 = float(lng2)
    radianLat1 = math.radians(lat1)
    radianLat2 = math.radians(nlat)
    radianLng1 = math.radians(lng1)
    part1 = math.sin(dst/(2*R)) ** 2
    part2 = (math.sin((radianLat2/2 - radianLat1/2))) ** 2
    part3 = math.cos(radianLat1) * math.cos(radianLat2)
    part3I = (part1- part2)/part3
    # print part3I
    # if part3I < 0:
    #     part3I *= -1
    part4 = math.asin(math.sqrt(part3I))

    # part3 = (math.sin((radianLng2/2 - radianLng1/2))) ** 2
    radianLng2 = (part4 + (radianLng1/2)) * 2
    degreeLat2 = math.degrees(radianLat2)
    degreeLng2 = math.degrees(radianLng2)
    # print degreeLat2, degreeLng2
    return degreeLng2


def mirroring(paln,list1):
    newlist = []
    if paln > list1[0][1]:
        for i in list1:
            dst = paln - i[1]
            newlg = paln + dst
            co = [i[0], newlg]
            newlist.append(co)
    else:
        for i in list1:
            dst = i[1] - paln
            newlg = paln - dst
            co = [i[0], newlg]
            newlist.append(co)
    return newlist




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

blat,blng = getBpoint(pa,dst)
perimeter = 2 * dst * math.pi
dst2 = perimeter/100
# getXpoint(pa,col2,dst,dst2)

radiusco = paln - blng
if radiusco < 0:
    radiusco *= -1
# print blat,blng
ratio = 5.00000001275/5.13450133924
radiusco = radiusco *ratio
ltdiameter  =radiusco * 2
# smallestLat = (palt - radiusco) * ratio
# biggestLat = (palt + radiusco) * ratio

smallestLat = (palt - radiusco)
biggestLat = (palt + radiusco)

step = ltdiameter/20
newlt = 0
lo = ["pa0" , str(paln), str(smallestLat)]
lo1 = ["pa1000" , str(paln), str(biggestLat)]
print ",".join(lo)
print ",".join(lo1)
rlist = []
for i in range(0,21):
    newlt = (i * step ) + smallestLat
    if i == 0:
        newlt = step/2 + smallestLat
        newlg = getXpoint(pa,newlt,dst)
        co = [newlt,newlg]
        rlist.append(co)

        newlt = step/6 + smallestLat
        newlg = getXpoint(pa,newlt,dst)
        co = [newlt,newlg]
        rlist.append(co)

        continue
    elif i == 20:
        newlt = (i -1 ) * step + smallestLat + step/2
        newlg = getXpoint(pa,newlt,dst)
        co = [newlt,newlg]
        rlist.append(co)

        newlt = (i -1 ) * step + smallestLat + (step -step/6)
        newlg = getXpoint(pa,newlt,dst)
        co = [newlt,newlg]
        rlist.append(co)

        continue

   
    # print newlt

    newlg = getXpoint(pa,newlt,dst)

    # print newlt, newlg
    co = [newlt,newlg]
    rlist.append(co)
    lo = ["pa" + str(i), str(newlg), str(newlt)]
    # print ",".join(lo)
getmirror = mirroring(paln,rlist)
newlist = rlist + getmirror

n = 1
for i in newlist:
    rlist.append(co)
    lo = ["pa" + str(n), str(i[1]), str(i[0])]
    print ",".join(lo)
    #if n == len(newlist):
    #    for m in range(1,len(newlist) + 1):
    #        print ",pa" + str(m)+",,,,,,,,,,"

    n+= 1





**** usage***

we want to get a list of latitude and longitude that form a circle from point x and with 7 km radius

$ python circle-from-latitudes.py "11.03122,105.0591" 7000
pa0,105.0591,10.9687625439
pa1000,105.0591,11.0936774561
pa1,105.080527726,10.9718854167
pa2,105.073180169,10.9698035015
pa3,105.087972657,10.9750082895
pa4,105.098111441,10.9812540351
pa5,105.105243388,10.9874997807
pa6,105.110632362,10.9937455263
pa7,105.114786656,10.9999912719
pa8,105.117968201,11.0062370176
pa9,105.12032878,11.0124827632
pa10,105.121960907,11.0187285088
pa11,105.122920456,11.0249742544
pa12,105.123237578,11.03122
pa13,105.122921812,11.0374657456
pa14,105.121963579,11.0437114912
pa15,105.120332683,11.0499572368
pa16,105.117973205,11.0562029824
pa17,105.114792573,11.0624487281
pa18,105.110638933,11.0686944737
pa19,105.105250252,11.0749402193
pa20,105.098118074,11.0811859649
pa21,105.08797818,11.0874317105
pa22,105.080532052,11.0905545833
pa23,105.073183112,11.0926364985
pa24,105.037672274,10.9718854167
pa25,105.045019831,10.9698035015
pa26,105.030227343,10.9750082895
pa27,105.020088559,10.9812540351
pa28,105.012956612,10.9874997807
pa29,105.007567638,10.9937455263
pa30,105.003413344,10.9999912719
pa31,105.000231799,11.0062370176
pa32,104.99787122,11.0124827632
pa33,104.996239093,11.0187285088
pa34,104.995279544,11.0249742544
pa35,104.994962422,11.03122
pa36,104.995278188,11.0374657456
pa37,104.996236421,11.0437114912
pa38,104.997867317,11.0499572368
pa39,105.000226795,11.0562029824
pa40,105.003407427,11.0624487281
pa41,105.007561067,11.0686944737
pa42,105.012949748,11.0749402193
pa43,105.020081926,11.0811859649
pa44,105.03022182,11.0874317105
pa45,105.037667948,11.0905545833
pa46,105.045016888,11.0926364985




- what it looks like in the maps

No comments:

Post a Comment