0

从坐标列表开始，我正在尝试创建一个包含插值坐标的新列表。我错过了一些东西，它只是一遍又一遍地附加第一个和最后一个坐标。

问题出在 main 函数中，与用新创建的点替换原点有关。我只包括其他人，因为它们是必要的。如果您运行它，它将创建一个 .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]]

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

4

3 回答 3

0

我修好了它。

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-20T03:03:13.503 回答
0

没有检查结果，但它吐出更多的点。我改变了你的流程。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)

滑动窗口来自Python 中的滚动或滑动窗口迭代器

于 2014-06-20T02:59:15.213 回答
0

您的代码中的这一行：

        if running_dist < point_dist:

看起来running_dist是0在那一点上，我想point_dist会大于0，这意味着循环的一部分将永远不会运行。

于 2014-06-20T01:39:05.363 回答