Python geopy.distance 模块，great_circle() 实例源码

def get_search_points(self, cell_id):
        points = []

        # For cell level 15
        for c in Cell(CellId(cell_id)).subdivide():
            for cc in c.subdivide():
                latlng = LatLng.from_point(cc.get_center())
                point = (latlng.lat().degrees, latlng.lng().degrees)
                points.append(point)

        points[0], points[1] = points[1], points[0]
        points[14], points[15] = points[15], points[14]
        point = points.pop(2)
        points.insert(7, point)
        point = points.pop(13)
        points.insert(8, point)

        closest = min(points, key=lambda p: great_circle(self.bot.position, p).meters)
        index = points.index(closest)

        return points[index:] + points[:index]

def retrieve_neighbors(index_center, df, spatial_threshold, temporal_threshold):
    neigborhood = []

    center_point = df.loc[index_center]

    # filter by time 
    min_time = center_point['date_time'] - timedelta(seconds=temporal_threshold)
    max_time = center_point['date_time'] + timedelta(seconds=temporal_threshold)
    df = df[(df['date_time'] >= min_time) & (df['date_time'] <= max_time)]

    # filter by distance
    for index, point in df.iterrows():
        if index != index_center:
            distance = great_circle(
                (center_point['latitude'], center_point['longitude']),
                (point['latitude'], point['longitude'])).meters
            if distance <= spatial_threshold:
                neigborhood.append(index)

    return neigborhood

def calculate_area(info):
    """
    Calculate the area of a simple rectangle (bounding box)
    :param info - the XML string containing data
    :return: area of box
    """
    bbox = info.find("./bbox")
    if bbox is not None:
        n = bbox.find("./north").text
        s = bbox.find("./south").text
        e = bbox.find("./east").text
        w = bbox.find("./west").text
        ns = great_circle((float(n), 0.0), (float(s), 0.0)).kilometers
        ew = great_circle((0.0, float(w)), (0.0, float(e))).kilometers
        return ns * ew
    else:
        return 1

def smart_map_point(s):
  '''Map point s to closets point on route conscious of nextStopID'''
  #Calculate range of indexes of points to check
  half = int(len(points) / 3)
  lat, lng, nextStopID = s
  print(s)
  destination_index = stops[nextStopID]
  if (destination_index - half) >= 0: #chuck to check is in higher half
    check = [i for i in range(destination_index - half, destination_index)]
  else: #chunk to check dips below 0
    check = [i for i in range(0, destination_index)]
    delta = -(destination_index - half)
    check.extend( [i for i in range(len(points) - delta, len(points))])
  check.extend([i % len(points) for i in range(destination_index,
    destination_index + 3)])

  #Actual search
  running_min = float("inf")
  point_index = None #index of closest point on map
  for i in check:
    dist = distance(points[i], s).miles
    if dist < running_min:
      running_min = dist
      point_index = i
  return point_index

def update_cluster_distance(self, cluster):
        cluster["distance"] = great_circle(self.bot.position, cluster["center"]).meters

def get_distance(self, location, fort):
        return great_circle(location, (fort["latitude"], fort["longitude"])).meters

def get_distance(self, location, pokemon):
        return great_circle(location, (pokemon["latitude"], pokemon["longitude"])).meters

def get_alt(self, at_point=None):
        if at_point is None:
            at_point = self._last_pos
        if self._elevation_at_point:
            elevations = sorted([(great_circle(at_point, k).meters, v, k) for k, v in self._elevation_at_point.items()])

            if len(elevations) == 1:
                return elevations[0][1]
            else:
                (distance_to_p1, ep1, p1), (distance_to_p2, ep2, p2) = elevations[:2]
                distance_p1_p2 = great_circle(p1, p2).meters
                return self._get_relative_hight(ep1, ep2, distance_p1_p2, distance_to_p1, distance_to_p2)
        else:
            return None

def batch_pdist(data_slice):
    # Each data_slice has tuples consisting of two points that we need to
    # find the great circle distance between and their weight:
    partial_sum = 0
    for X, Y, weights in data_slice:
        dist = np.array([])
        zipped = zip(X, Y)
        for x, y in zipped:
            dist = np.append(dist, great_circle(x, y).km)
        partial_sum += np.sum(weights * dist )
    return partial_sum
    # return 10

def dfun(u, v):
        return great_circle(u, v).km

def dfun(a, b):  # example of dist() between 2 numbers
    # from math import sqrt
    # return sqrt(abs(a-b))
    from geopy.distance import great_circle
    return great_circle(a, b).km

def landed_check(coordinate_list, r):
    '''Check if the plane has not moved fast enough recently to be airborne.'''
    # first do naive check - two points distance / time
    naive_distance = great_circle(coordinate_list[0][0:1], coordinate_list[-1][0:1]).meters
    time_delta = coordinate_list[-1][3] - coordinate_list[0][3]
    if time_delta < 200:
        return 0
    elif (naive_distance / time_delta) < 8.94:
        # logger.debug("Something is moving too slow to be flying!")
        return airport_proximity_check((coordinate_list[-1][0], coordinate_list[-1][1]), r)
    else:
        return 0

def filterbylocation(self,user_location):
        user_lat_lon = (user_location[0], user_location[1])
        pok_loc = (float(self.latitude), float(self.longitude))
        return great_circle(user_lat_lon, pok_loc).km <= user_location[2]

def test_dist_itslef(self):
        """Compary cdist to GeoPy."""
        X = [[0, 10], [4, 2]]  # Just some points. I've no idea where on globe.
        c = cdist(X[0], X[1])
        string_geopy = '{}'.format(great_circle(X[0], X[1]))
        float_geopy = float(string_geopy[:-3])
        self.assertTrue(np.round(c) == np.round(float_geopy))

        X = [[34.0522, 118.2437],  # Lon Angeles
             [37.7749, 122.4194]]  # San Francisco
        c = cdist(X[0], X[1])
        string_geopy = '{}'.format(great_circle(X[0], X[1]))
        float_geopy = float(string_geopy[:-3])
        self.assertTrue(np.round(c) == np.round(float_geopy))

def geodistance( coords1 , coords2 ):
  lat1 , lon1 = coords1[ : 2]
  lat2 , lon2 = coords2[ : 2]
  try: return distance.vincenty( ( lat1 , lon1 ) , ( lat2 , lon2 ) ).meters / 1000.0
  except: return distance.great_circle( ( lat1 , lon1 ) , ( lat2 , lon2 ) ).meters / 1000.0

def gc(a, b, c, d):
    p1 = (b, a)
    p2 = (d, c)
    return great_circle(p1, p2).kilometers

def _ReadAndFilterData(filename, start_date, end_date, pos, distance):
  """Creates a filter function for the CSV reader, and reads csv data.

  Args:
    filename: The path of the file to read.
    start_date: Start date we should begin filtering the data.
    end_date: End date we should begin filtering the data.
    pos: Location we should be filtering the data for.
    distance: Distance in KM for which we include aftershocks.

  Returns:
    List of dictionaries of ISC data.
  """
  def _Filter(x):
    """Filter that we apply to all isc data."""
    try:
      # Remove the normal problems with the data.
      if not _IsRowValid(x): return False
      # Make sure the data point is in the date range specified.
      if not start_date <= x['date_time'] <= end_date: return False
      # Make sure the data point is within the distance specified.
      if vincenty((x['lat'], x['lon']), pos) > distance: return False
    except ValueError:
      # Vincenty doesn't always converge. If it fails to, then default to a
      # Great Circle distance.
      if great_circle((x['lat'], x['lon']), pos) > distance: return False
    return True
  return _ReadCsvFile(filename, data_validator=_Filter)

def check_median_absolute_deviation(data_points, median):
    """Calculate Median Absolute Deviation of a set of points."""
    distances = []
    for i in range(0, len(data_points)):
        try:
            distances.append(vincenty(median, data_points[i]).kilometers)
        except ValueError:
            # Vincenty doesn't always converge so fall back on great circle distance which is less accurate but always converges
            distances.append(great_circle(median, data_points[i]).kilometers)
    return(numpy.median(distances))

def denomsum(test_median, data_points):
    """Provides the denominator of the weiszfeld algorithm."""
    temp = 0.0
    for i in range(0, len(data_points)):
        try:
            temp += 1 / vincenty(test_median, data_points[i]).kilometers
        except ZeroDivisionError:
            continue  # filter points that equal the median out (otherwise no convergence)
        except ValueError:
            # Vincenty doesn't always converge so fall back on great circle distance which is less accurate but always converges
            temp += 1 / great_circle(test_median, data_points[i]).kilometers
    return temp

def numersum(test_median, data_point):
    """Provides the denominator of the weiszfeld algorithm depending on whether you are adjusting the candidate x or y."""
    try:
        return 1 / vincenty(test_median, data_point).kilometers
    except ZeroDivisionError:
        return 0  # filter points that equal the median out (otherwise no convergence)
    except ValueError:
        # Vincenty doesn't always converge so fall back on great circle distance which is less accurate but always converges
        return 1 / great_circle(test_median, data_point).kilometers

def objfunc(test_median, data_points):
    """This function calculates the sum of linear distances from the current candidate median to all points
    in the data set, as such it is the objective function that we are minimising.
    """
    temp = 0.0
    for i in range(0, len(data_points)):
        try:
            temp += vincenty(test_median, data_points[i]).kilometers
        except ValueError:
            # Vincenty doesn't always converge so fall back on great circle distance which is less accurate but always converges
            temp += great_circle(test_median, data_points[i]).kilometers
    return temp

def load_graph():
  '''Calculate graph of distances between route points
  Each index i holds distance from i to i+1'''
  dists = []
  for i in range(len(points)):
    d = distance(points[i], points[(i+1) % len(points)])
    dists.append(d.miles)
  return dists

def sum_graph(i, j):
  '''sum graph [i,j) to determine distance b/w indicies i and j
  Automatically perform wraparound'''
  if i < j: #normal case, no wraparound
    return sum(graph[i:j])
  elif i > j: #destination point is earlier in list, wrap around
    return sum(graph[i:]) + sum(graph[:j])
  else:
      return 0

def map_point(s): 
  '''Map point s to closest point on route, return index of route point'''
  running_min = float("inf")
  point_index = None #index of closest point on map
  for i, p in enumerate(points):
    dist = distance(p, s).miles
    if dist < running_min:
      running_min = dist
      point_index = i
  return point_index

def compare_coordinates(cls, coords1, coords2, *args, tolerance=None,
                            unit='km', **kwargs):
        """
        Check if a distance between the pairs of coordinates provided is
        within the specified tolerance.

        Return True if yes, otherwise return False.

        Use geopy (https://pypi.python.org/pypi/geopy).

        Try to use Vincenty formula, if error occurs use Great Circle formula.

        :param coords1: pair of coordinates - a tuple of two numbers
        :param coords2: pair of coordinates - a tuple of two numbers
        :param tolerance: number
        :param unit: str, one of: 'kilometers', 'km', 'meters',
                                  'm', 'miles', 'mi', 'feet',
                                  'ft', 'nautical', 'nm'
        :rtype: bool

        :Example:

        >>> a, b = (36.1332600, -5.4505100), (35.8893300, -5.3197900)
        >>> MatchBlock.compare_coordinates(a, b, tolerance=20, unit='mi')
        True

        >>> a, b = (36.1332600, -5.4505100), (35.8893300, -5.3197900)
        >>> MatchBlock.compare_coordinates(a, b, tolerance=1, unit='mi')
        False
        """

        if tolerance is None:
            tolerance = cls.coordinates_tolerance

        units = {'kilometers', 'km', 'meters', 'm', 'miles', 'mi',
                 'feet', 'ft', 'nautical', 'nm'}

        unit = unit.strip().lower()
        if unit not in units:
            raise ValueError('unsupported unit')

        try:
            length = getattr(vincenty(coords1, coords2, *args, **kwargs), unit)
        except ValueError as e:
            if 'Vincenty formula failed to converge!' in e.args:
                warnings.warn(
                    'vincenty formula failed, using great circle formula')
                length = getattr(great_circle(coords1, coords2, *args), unit)
            else:
                raise

        return length <= tolerance

def calculate_scores(predicted, gold, inspect=False, topocluster=False):
    """
    Given the predictions and the gold annotations, calculate precision, recall, F Score and accuracy.
    :param topocluster: Topocluster geoparser produces NON-STANDARD output so has to be treated differently
    :param inspect: If True, the differences between gold and predicted files will be printed
    :param predicted: path to the file with parser predictions
    :param gold: path to the file with gold annotations
    :return: a list of errors per toponym i.e how far away is each correctly identified toponym from
    the gold location. This is used to measure the accuracy of the geocoding part
    """
    tp, fp, fn = 0.0, 0.0, 0.0
    accuracy = {}
    wiki = True if "wiki" in predicted else False
    predictions_file = codecs.open(predicted)
    gold = codecs.open(gold)
    toponym_index = -1
    for predicted, gold in zip(predictions_file, gold):
        predicted_tops = predicted.split("||")[:-1]
        gold_tops = gold.split("||")[:-1]
        for gold_top in gold_tops[:]:
            toponym_index += 1
            gold_top_items = gold_top.split(",,")
            for predicted_top in predicted_tops[:]:
                predicted_top_items = predicted_top.split(",,")
                mean_g = (int(gold_top_items[4]) + int(gold_top_items[5])) / 2.0
                mean_p = (int(predicted_top_items[4]) + int(predicted_top_items[5])) / 2.0
                #  If the toponym position (its mean) is no more than 9 characters from gold AND the two
                #  strings are equal then it's a match. For reasons to do with UTF-8 encoding and decoding,
                #  the toponym indices may, in a few instances, be off by a few positions when using Web APIs.
                match = False  # A flag to establish whether this is a matching prediction
                if topocluster:  # Only match for the toponym name in this instance
                    if predicted_top_items[1].lower() == gold_top_items[1].lower():
                        match = True
                elif abs(mean_g - mean_p) < 10 and predicted_top_items[1].lower() == gold_top_items[1].lower():
                    match = True   # Change the number above to 0 for EXACT matches, 10 for INEXACT matches
                if match:
                    tp += 1
                    predicted_tops.remove(predicted_top)
                    gold_tops.remove(gold_top)
                    predicted_coord = (float(predicted_top_items[2]), float(predicted_top_items[3]))
                    gold_coord = (float(gold_top_items[2]), float(gold_top_items[3]))
                    accuracy[toponym_index] = numpy.log(1 + great_circle(predicted_coord, gold_coord).kilometers)
                    break
        if not wiki:
            fp += len(predicted_tops)
        fn += len(gold_tops)
        if inspect:
            if len(predicted_tops) > 0 or 0 < len(gold_tops):
                print "Predicted:", " - ".join(predicted_tops)
                print "Gold Tops:", " - ".join(gold_tops)
    f_score = (tp, fp, fn)
    output = {"f_score": f_score, "accuracy": accuracy}
    return output

def find_station(*args):
  stns=0
  if len(args) ==1:
    station_name=args[0]
    print("LOOKUP BY STATION NAME: ",station_name)
    station_name=station_name.upper()
    ghcnd_stations=gp.get_ghcnd_stations()

    stns=filter(lambda x: re.search(station_name,x), ghcnd_stations[:,5])
    print("GHCND ID          LAT        LON    ELEV  ST       STATION NAME")
    print("###############################################################")
    for station_counter in xrange(len(stns)):
      ghcnd_meta = ghcnd_stations[ghcnd_stations[:,5]== stns[station_counter]]
      print(ghcnd_meta[0][0],ghcnd_meta[0][1],ghcnd_meta[0][2],ghcnd_meta[0][3],ghcnd_meta[0][4],ghcnd_meta[0][5])

  elif len(args)==3:
    station_lat=args[0]
    station_lon=args[1]
    distance_limit=args[2]

    print("LOOKUP BY STATION LAT: ",station_lat," LON: ",station_lon, " DIST LIMIT (mi): ",distance_limit)

    target_latlon = (float(station_lat), float(station_lon))
    ghcnd_stations=gp.get_ghcnd_stations()

    print("GHCND ID          LAT        LON    ELEV  ST       STATION NAME")
    print("###############################################################")
    for ghcnd_counter in xrange(ghcnd_stations[:,0].size):
      candidate_latlon=(ghcnd_stations[ghcnd_counter][1], ghcnd_stations[ghcnd_counter][2])
      dist=great_circle(target_latlon, candidate_latlon).miles
      if dist <= distance_limit:
        print(ghcnd_stations[ghcnd_counter][0],ghcnd_stations[ghcnd_counter][1],ghcnd_stations[ghcnd_counter][2],ghcnd_stations[ghcnd_counter][3],ghcnd_stations[ghcnd_counter][4],ghcnd_stations[ghcnd_counter][5])

  else:
    print("USAGE\n  NAME or\n  LAT LON DIST")
    return None
  return None

#################################################
# MODULE: get_metadata
# Get Metadata From Station
# 2 sources
#    - ghcnd-stations.txt
#    - Historical Observing Metadata Repository
#      (HOMR)
#################################################

def compute_user_median(data_points, num_iter, csvwriter, current_uid):
    if len(data_points) < LIMIT_POINTS:  # Insufficient points for the user - don't record median
        if OUTPUT_ALL_USERS:
            csvwriter.writerow([current_uid, None])
    else:
        if SNAP_TO_USER_POINTS: # ensure median is one of the user's points
            lowest_dev = float("inf")
            for point in data_points:
                tmp_abs_dev = objfunc(point, data_points)
                if tmp_abs_dev < lowest_dev:
                    lowest_dev = tmp_abs_dev
                    test_median = point
        else:
            test_median = cand_median(data_points)  # Calculate centroid more or less as starting point
            if objfunc(test_median, data_points) != 0:  # points aren't all the same
                # iterate to find reasonable estimate of median
                for x in range(0, num_iter):
                    denom = denomsum(test_median, data_points)
                    next_lat = 0.0
                    next_lon = 0.0

                    for y in range(0, len(data_points)):
                        next_lat += (data_points[y][0] * numersum(test_median, data_points[y])) / denom
                        next_lon += (data_points[y][1] * numersum(test_median, data_points[y])) / denom

                    prev_median = test_median
                    test_median = (next_lat, next_lon)
                    try:
                        if vincenty(prev_median, test_median).meters < DISTANCE_THRESHOLD:
                            break
                    except:
                        if great_circle(prev_median, test_median).meters < DISTANCE_THRESHOLD:
                            break

                if x == num_iter - 1:
                    print('{0}: failed to converge. Last change between iterations was {1} meters.'.format(current_uid, great_circle(prev_median, test_median).meters))

        # Check if user points are under the limit median absolute deviation
        if check_median_absolute_deviation(data_points, test_median) <= LIMIT_MAD:
            csvwriter.writerow([current_uid, (round(test_median[0],6), round(test_median[1],6))])
        else:
            if OUTPUT_ALL_USERS:
                csvwriter.writerow([current_uid, None])