在我的循环中需要另一组眼睛

从坐标列表开始，我试图创建一个包含插值坐标的新列表。我错过了一些东西，它只是反复追加第一个和最后一个坐标。在我的循环中需要另一组眼睛

问题出在主要功能上，并且与用新创建的点替换原点有关。我只包括其他人，因为他们是必要的。如果你运行这个，它会创建一个.kml文件，你可以在Google Earth中打开并查看问题。

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()

目前，新列表看起来像这样。你可以看到原点和目的地一遍又一遍地被追加，而没有附加新的点。

[[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]]

预期结果：以3米为间隔的原点和目的点之间的坐标列表（包括起点和终点）。

来源

2014-06-20 Will Luce

你要知道，什么是错的结果。要求其他人在您提供的脚本输出中启动Google地球并不是获得真正答案的非常有效的方法。 –

我也这么认为，但它不是要求某人盯着24套经纬度/经度对的列表，以查看它只是不断追加起源和目的地。我编辑它以显示两者。 –

当您在main（）中写入destination = coord_list [1]'时，是否指'destination = coord_list [x]'？ –

我固定它。

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

是需要：

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

来源

2014-06-20 03:03:13

这行代码：

 if running_dist < point_dist:

看起来running_dist是0在这一点上，和我想象point_dist会比0更大，这意味着环路的一部分将永远不会运行。

来源

2014-06-20 01:39:05

这并不完全正确。 running_dist从0开始并增加.003，直到达到point_dist。 –

没有检查结果，但它吐出了更多点。我改变了你的流程。 while循环处理似乎关闭。我重新使用滑动窗口，并移动while循环来仅用于原点和目标部分之间的插值。

interpolate_points()取代您的大部分main()功能。 我已经删除KML处理，所以我并不需要安装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)

滑动窗口是Rolling or sliding window iterator in Python

来源

2014-06-20 02:59:15 monkut

在我的循环中需要另一组眼睛

回答

相关问题