Starting with a list of coordinates, I'm trying to create a new list with interpolated coordinates included. I'm missing something and it just appends the first and last coordinate over and over again.
The problem is in the main function and has something to do with replacing the origin point with the newly created point. I've only included the others because they're necessary. If you run this, it will create a .kml file that you can open in google earth and see the problem.
import math
from geopy import distance
import simplekml
def build_kml_points(filename, coord_list):
kml = simplekml.Kml()
name = 1
for coord_pair in coord_list:
kml.newpoint(name=str(name), coords=[(coord_pair[1], coord_pair[0])])
name += 1
kml.save(str(filename))
def bearing(pointA, pointB):
# Calculates the bearing between two points.
#
# :Parameters:
# - `pointA: The tuple representing the latitude/longitude for the
# first point. Latitude and longitude must be in decimal degrees
# - `pointB: The tuple representing the latitude/longitude for the
# second point. Latitude and longitude must be in decimal degrees
#
# :Returns:
# The bearing in degrees
#
# :Returns Type:
# float
# if (type(pointA) != tuple) or (type(pointB) != tuple):
# raise TypeError("Only tuples are supported as arguments")
lat1 = math.radians(pointA[0])
lat2 = math.radians(pointB[0])
diffLong = math.radians(pointB[1] - pointA[1])
x = math.sin(diffLong) * math.cos(lat2)
y = math.cos(lat1) * math.sin(lat2) - (math.sin(lat1)
* math.cos(lat2) * math.cos(diffLong))
initial_bearing = math.atan2(x, y)
# Now we have the initial bearing but math.atan2 return values
# from -180 to + 180 which is not what we want for a compass bearing
# The solution is to normalize the initial bearing as shown below
initial_bearing = math.degrees(initial_bearing)
compass_bearing = (initial_bearing + 360) % 360
return compass_bearing
# Vincenty's Direct formulae
def vinc_pt(phi1, lembda1, alpha12, s ) :
"""
Returns the lat and long of projected point and reverse azimuth
given a reference point and a distance and azimuth to project.
lats, longs and azimuths are passed in decimal degrees
Returns ( phi2, lambda2, alpha21 ) as a tuple
f = flattening of the ellipsoid: 1/298.277223563
a = length of the semi-major axis (radius at equator: 6378137.0)
phi1 = latitude of the starting point
lembda1 = longitude of the starting point
alpha12 = azimuth (bearing) at the starting point
s = length to project to next point
"""
f = 1/298.277223563
a = 6378137.0
piD4 = math.atan( 1.0 )
two_pi = piD4 * 8.0
phi1 = phi1 * piD4 / 45.0
lembda1 = lembda1 * piD4 / 45.0
alpha12 = alpha12 * piD4 / 45.0
if ( alpha12 < 0.0 ) :
alpha12 = alpha12 + two_pi
if ( alpha12 > two_pi ) :
alpha12 = alpha12 - two_pi
# length of the semi-minor axis (radius at the poles)
b = a * (1.0 - f)
TanU1 = (1-f) * math.tan(phi1)
U1 = math.atan( TanU1 )
sigma1 = math.atan2( TanU1, math.cos(alpha12) )
Sinalpha = math.cos(U1) * math.sin(alpha12)
cosalpha_sq = 1.0 - Sinalpha * Sinalpha
u2 = cosalpha_sq * (a * a - b * b ) / (b * b)
A = 1.0 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * \
(320 - 175 * u2) ) )
B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2) ) )
# Starting with the approx
sigma = (s / (b * A))
last_sigma = 2.0 * sigma + 2.0 # something impossible
# Iterate the following 3 eqs unitl no sig change in sigma
# two_sigma_m , delta_sigma
while ( abs( (last_sigma - sigma) / sigma) > 1.0e-9 ) :
two_sigma_m = 2 * sigma1 + sigma
delta_sigma = B * math.sin(sigma) * ( math.cos(two_sigma_m) \
+ (B/4) * (math.cos(sigma) * \
(-1 + 2 * math.pow( math.cos(two_sigma_m), 2 ) - \
(B/6) * math.cos(two_sigma_m) * \
(-3 + 4 * math.pow(math.sin(sigma), 2 )) * \
(-3 + 4 * math.pow( math.cos (two_sigma_m), 2 ))))) \
last_sigma = sigma
sigma = (s / (b * A)) + delta_sigma
phi2 = math.atan2 ( (math.sin(U1) * math.cos(sigma) + math.cos(U1) * math.sin(sigma) * math.cos(alpha12) ), \
((1-f) * math.sqrt( math.pow(Sinalpha, 2) +
pow(math.sin(U1) * math.sin(sigma) - math.cos(U1) * math.cos(sigma) * math.cos(alpha12), 2))))
lembda = math.atan2( (math.sin(sigma) * math.sin(alpha12 )), (math.cos(U1) * math.cos(sigma) -
math.sin(U1) * math.sin(sigma) * math.cos(alpha12)))
C = (f/16) * cosalpha_sq * (4 + f * (4 - 3 * cosalpha_sq ))
omega = lembda - (1-C) * f * Sinalpha * \
(sigma + C * math.sin(sigma) * (math.cos(two_sigma_m) +
C * math.cos(sigma) * (-1 + 2 * math.pow(math.cos(two_sigma_m),2) )))
lembda2 = lembda1 + omega
alpha21 = math.atan2 ( Sinalpha, (-math.sin(U1) * math.sin(sigma) +
math.cos(U1) * math.cos(sigma) * math.cos(alpha12)))
alpha21 = alpha21 + two_pi / 2.0
if ( alpha21 < 0.0 ) :
alpha21 = alpha21 + two_pi
if ( alpha21 > two_pi ) :
alpha21 = alpha21 - two_pi
phi2 = phi2 * 45.0 / piD4
lembda2 = lembda2 * 45.0 / piD4
alpha21 = alpha21 * 45.0 / piD4
return phi2, lembda2, alpha21
def main():
coord_list = [[40.081059133213, -105.28215], [40.081188699819, -105.28215], [40.081318266425, -105.28215]]
point_list = []
x = 1
running_dist = 0
while x < 3:
origin = coord_list[x-1]
destination = coord_list[x]
# append the point from the original list
point_list.append(origin)
point_dist = distance.distance(origin, destination).km
point_dist = float(point_dist[:-3])
init_bearing = bearing(origin, destination)
if running_dist < point_dist:
new_point = vinc_pt(origin[0], origin[1], init_bearing, 3)
point_list.append([new_point[0], new_point[1]])
running_dist += .003
else:
x += 1
running_dist = 0
point_list.append(destination)
build_kml_points('Test.kml', point_list)
main()
Currently, the new list looks like this. You can see that the origin and destination are appended over and over again without appending new points.
[[40.081059133213, -105.28215], [40.08108615142624, -105.28215], [40.081059133213, -105.28215], [40.08108615142624, -105.28215], [40.081059133213, -105.28215], [40.08108615142624, -105.28215], [40.081059133213, -105.28215], [40.08108615142624, -105.28215], [40.081059133213, -105.28215], [40.08108615142624, -105.28215], [40.081059133213, -105.28215], [40.081188699819, -105.28215], [40.081318266425, -105.28215], [40.08129124821285, -105.28215], [40.081318266425, -105.28215], [40.08129124821285, -105.28215], [40.081318266425, -105.28215], [40.08129124821285, -105.28215], [40.081318266425, -105.28215], [40.08129124821285, -105.28215], [40.081318266425, -105.28215], [40.08129124821285, -105.28215], [40.081318266425, -105.28215], [40.081188699819, -105.28215]]
Expected result: a list of coordinates (including the origin and destination) between the origin and destination at 3m intervals.
I fixed it.
if running_dist < point_dist:
new_point = vinc_pt(origin[0], origin[1], init_bearing, 3)
point_list.append([new_point[0], new_point[1]])
running_dist += .003
else:
x += 1
running_dist = 0
needed to be:
while running_dist < point_dist:
new_point = vinc_pt(origin[0], origin[1], init_bearing, 3)
print new_point
origin = [new_point[0], new_point[1]]
point_list.append([new_point[0], new_point[1]])
running_dist += .003
print running_dist
x += 1
running_dist = 0