0

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.

Will Luce
  • 1,781
  • 3
  • 20
  • 33
  • 1
    You shall know, what is wrong with your result. Asking other to start up Google Earth on the output of a script you have provided is not very effective method to get real answer. – Jan Vlcinsky Jun 20 '14 at 01:07
  • I thought so too, but it's either that or asking someone to stare at a list of 24 sets of lat/lon pairs to see that it is just constantly appending the origin and destination. I edited it to show both. – Will Luce Jun 20 '14 at 01:10
  • 1
    Did you mean `destination = coord_list[x]` when you wrote `destination = coord_list[1]` in `main()`? – R Sahu Jun 20 '14 at 01:15
  • Yes, thank you for noticing, but the problem remains. – Will Luce Jun 20 '14 at 01:20
  • What result do you expect ? – furas Jun 20 '14 at 01:25
  • I expect a list of coordinates (including the origin and destination) between the origin and destination at 3m intervals. – Will Luce Jun 20 '14 at 01:27
  • Add expected result to your question . You can use `print` in code to see what is going on with variables. – furas Jun 20 '14 at 01:29
  • The `Distance` object returned by distance.distance has a `km` attribute. There's no need to slice and convert to float, just access it directly. `point_dist = distance.distance(origin, destination).km` – monkut Jun 20 '14 at 01:33
  • When I wrote `Add expected result to your question` I thought about list with numbers :) – furas Jun 20 '14 at 01:36
  • Thanks monkut. I fixed it, but the problem remains. – Will Luce Jun 20 '14 at 01:36
  • I know you did furas, but if I could find a way to print the list of numbers I'm looking for, I wouldn't be asking the question. I know this is some silly loop thing that I'm missing. I'm a competent programmer, but I need another set of eyes on a dumb problem. – Will Luce Jun 20 '14 at 01:38

3 Answers3

0

This line in your code:

        if running_dist < point_dist:

It looks like running_dist is 0 at that point, and I imagine point_dist will be greater than 0, meaning that part of the loop will never run.

Ethan Furman
  • 63,992
  • 20
  • 159
  • 237
  • That's not quite true. running_dist starts at 0 and is incremented by .003 until it reaches point_dist. – Will Luce Jun 20 '14 at 01:45
0

Didn't check the result, but it's spitting out more points. I changed your flow around. The while loop handling seemed to be off. I reworked it to use a sliding window, and moved the while loop to be used for the interpolation between origin and destination portion only.

interpolate_points() replaces most of your main() function. I've removed kml processing so I don't need to install simplekml

import math
from geopy import distance
from itertools import islice

def window(seq, n=2):
    "Returns a sliding window (of width n) over data from the iterable"
    "   s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...                   "
    it = iter(seq)
    result = tuple(islice(it, n))
    if len(result) == n:
        yield result    
    for elem in it:
        result = result[1:] + (elem,)
        yield result


[...]

def interpolate_points(coord_list, meters_increment=3):
    point_list = []
    for origin, destination in window(coord_list, 2):
        distance_to_destination = distance.distance(origin, destination).km
        # points need to be tuples for bearing function
        origin = tuple(origin)
        destination = tuple(destination)        
        init_bearing = bearing(origin, destination)

        # process and create points between original orign and distance
        remaining_distance = 0
        while remaining_distance < distance_to_destination:            
            point_list.append(origin)
            lat, lon, azimuth = vinc_pt(origin[0], origin[1], init_bearing, meters_increment)
            origin = (lat, lon)
            # recheck distance
            remaining_distance = distance.distance(origin, destination).km
    # with sliding window desination becomes the origin
    # --> make sure last point gets added
    last_coord = tuple(coord_list[-1]) # convert to tuple to match others 
    if last_coord not in point_list:
        point_list.append(last_coord)
    return point_list

coord_list = [[40.081059133213, -105.28215], [40.081188699819, -105.28215], [40.081318266425, -105.28215]]
interpolate_points(coord_list)

sliding window is from Rolling or sliding window iterator in Python

Community
  • 1
  • 1
monkut
  • 42,176
  • 24
  • 124
  • 155
0

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
Will Luce
  • 1,781
  • 3
  • 20
  • 33