Wednesday, May 17, 2017

get a list of latitude and longitude of a straight between two coordination point

*** script ***
#~/usr/bin/python

import sys
import math

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

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

pbl = pb.split(",")
pblt = float(pbl[0])
pbln = float(pbl[1])

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 equationxp(a,b):
    x2 = a ** 2
    x = a * b * 2
    c = b ** 2
    eq = [x2,x,c]
    return eq
# line equation y  = mx + b
slop = (pbln - paln)/(pblt - palt)
b = pbln -(slop * pblt)
# print slop
# print b
eq1 = equationxp(1,(-1)*palt)
eq2 = equationxp(slop,b-paln)
# Quadratic equation
feq = [eq1[0] + eq2[0], eq1[1]+eq2[1], eq1[2]+eq2[2] -dst]
# print feq
sqtB2ac4 =((feq[1] ** 2) - (4 * feq[0] * feq[2]))**(0.5)
x1 = (((-1) * feq[1]) - ( sqtB2ac4)) /(2 * feq[0])
x2 = (((-1) * feq[1]) + (sqtB2ac4)) /(2 * feq[0])

if x2 > x1:
    npx = x2
    npy = (slop * npx) - b
else:
    npx = x1
    npy = (slop * npx) - b

# print npx, npy
d = getDistance2(pa,pb)

daxis = palt - pblt
if daxis < 0:
    daxis *=  -1

nung = round(d/dst,0)
# print nung
distg = round(daxis/nung,6)

for i in range(1,int(nung)):
    if palt > pblt:
        newp = pblt
    else:
        newp =palt
    adding = i * distg
   
    newx = newp + adding
    newy = slop * newx + b
    print str(newx)+","+str(newy)



*** example

$ python line-between-latitudes.py "13.0019,103.201" "13.0973,102.939" 5000
13.0178,103.157333333
13.0337,103.113666667
13.0496,103.07
13.0655,103.026333333
13.0814,102.982666667

No comments:

Post a Comment