Python geopy.distance.vincenty() Examples

def calculate_transition_cost(self, source, target):
        max_route_distance = self.calculate_max_route_distance(
            source.measurement, target.measurement)
        try:
            _, route_distance = road_routing.road_network_route(
                (source.edge, source.location),
                (target.edge, target.location),
                self.get_road_edges,
                max_path_cost=max_route_distance)
        except shortest_path.PathNotFound as err:
            # Not reachable
            return -1
        # Geodesic distance based on WGS 84 spheroid
        great_circle_distance = vincenty(
            (source.measurement.lat, source.measurement.lon),
            (target.measurement.lat, target.measurement.lon)).meters
        delta = abs(route_distance - great_circle_distance)
        return delta / self.beta 

def get_fort_ids_within_range(session, forts, range, lat, lon):

    for cache in DBCache.fort_ids_within_range:
        if cache.range == range and cache.lat == lat and cache.lon == lon:
            return cache.ids

    if forts is None:
        forts = get_forts(session)

    ids_with_range = []
    for fort in forts:
        distance = vincenty((fort.lat, fort.lon), (lat, lon)).meters
        if distance <= range:
            ids_with_range.append([distance, fort.id])

    session.commit()

    ids_with_range = sorted(ids_with_range, key=lambda x: x[0], reverse=False)
    ids = [obj[1] for obj in ids_with_range]
    cache_object = DBCacheFortIdsWithinRange(range, lat, lon, ids)
    DBCache.fort_ids_within_range.append(cache_object)

    return ids 

def closest_point(location, location_dict):
    """ take a tuple of latitude and longitude and 
        compare to a dictonary of locations where
        key = location name and value = (lat, long)
        returns tuple of (closest_location , distance) """
    closest_location = None
    for city in location_dict.keys():
        distance = vincenty(location, location_dict[city]).kilometers
        if closest_location is None:
            closest_location = (city, distance)
        elif distance < closest_location[1]:
            closest_location = (city, distance)
    return closest_location

#test = (39.524325, -122.293592) #likely 'Willows'
#closest_point(test, city_coords)


#run number 2 to determine both the nearest city, and then
	#also the nearest city with 1million people (subset the original dict) 

def get_circle_centers(b1, b2, radius):
    """
    the function covers the area within the bounds with circles

    :param b1: south-west bounds [lat, lng]
    :param b2: north-east bounds [lat, lng]
    :param radius: specified radius, adapt for high density areas
    :return: list of circle centers that cover the area between lower/upper
    """

    sw = Point(b1)
    ne = Point(b2)

    # north/east distances
    dist_lat = vincenty(Point(sw[0], sw[1]), Point(ne[0], sw[1])).meters
    dist_lng = vincenty(Point(sw[0], sw[1]), Point(sw[0], ne[1])).meters

    circles = cover_rect_with_cicles(dist_lat, dist_lng, radius)
    cords = [
        VincentyDistance(meters=c[0])
        .destination(
            VincentyDistance(meters=c[1])
            .destination(point=sw, bearing=90),
            bearing=0
        )[:2]
        for c in circles
    ]

    return cords 

def get_distance(self):
        position = Pokeconfig.get().position
        distance = vincenty(position, self.position)
        if Pokeconfig.get().distance_unit == 'meters':
            return distance.meters
        else:
            return distance.miles 

def calculate_transition_costs(self, source, targets):
        if not targets:
            return []

        # All measurements in targets should be the same, since they
        # are grouped by measurement ID
        target_measurement = targets[0].measurement

        max_route_distance = self.calculate_max_route_distance(
            source.measurement, target_measurement)
        route_results = road_routing.road_network_route_many(
            (source.edge, source.location),
            [(tc.edge, tc.location) for tc in targets],
            self.get_road_edges,
            max_path_cost=max_route_distance)

        # Geodesic distance based on WGS 84 spheroid
        great_circle_distance = vincenty(
            (source.measurement.lat, source.measurement.lon),
            (target_measurement.lat, target_measurement.lon)).meters

        costs = []
        for target, (path, route_distance) in zip(targets, route_results):
            if route_distance < 0:
                # Not reachable
                costs.append(-1)
                continue
            target.path[source] = path
            delta = abs(route_distance - great_circle_distance)
            costs.append(math.exp(-delta / self.beta) / self.beta)

        # They should be equivalent (only for testing):
        # for cost, target in zip(costs, targets):
        #     single_cost = self.calculate_transition_cost(source, target)
        #     assert abs(cost - single_cost) < 0.0000001

        return costs 

def calculate_transition_costs(self, source, targets):
        if not targets:
            return []

        # All measurements in targets should be the same, since they
        # are grouped by measurement ID
        target_measurement = targets[0].measurement

        max_route_distance = self.calculate_max_route_distance(
            source.measurement, target_measurement)
        route_results = road_routing.road_network_route_many(
            (source.edge, source.location),
            [(tc.edge, tc.location) for tc in targets],
            self.get_road_edges,
            max_path_cost=max_route_distance)

        # Geodesic distance based on WGS 84 spheroid
        great_circle_distance = vincenty(
            (source.measurement.lat, source.measurement.lon),
            (target_measurement.lat, target_measurement.lon)).meters

        costs = []
        for target, (path, route_distance) in zip(targets, route_results):
            if route_distance < 0:
                # Not reachable
                costs.append(-1)
                continue
            target.path[source] = path
            delta = abs(route_distance - great_circle_distance)
            costs.append(delta / self.beta)

        # They should be equivalent (only for testing):
        # for cost, target in zip(costs, targets):
        #     single_cost = self.calculate_transition_cost(source, target)
        #     assert abs(cost - single_cost) < 0.0000001

        return costs 

def points_between(start, end, numpoints, constantvalue=False):
    distance = []
    latitude = []
    longitude = []

    lat0 = start.latitude
    lon0 = start.longitude
    lat1 = end.latitude
    lon1 = end.longitude

    if constantvalue and np.isclose(lat0, lat1):
        latitude = np.ones(numpoints) * lat0
        longitude = np.linspace(lon0, lon1, num=numpoints)
        for lon in longitude:
            distance.append(vincenty(start, geopy.Point(lat0, lon)).km)
        if lon1 > lon0:
            b = 90
        else:
            b = -90
    elif constantvalue and np.isclose(lon0, lon1):
        latitude = np.linspace(lat0, lat1, num=numpoints)
        longitude = np.ones(numpoints) * lon0
        for lat in latitude:
            distance.append(vincenty(start, geopy.Point(lat, lon0)).km)
        if lat1 > lat0:
            b = 0
        else:
            b = 180
    else:
        total_distance = vincenty(start, end).km
        distance = np.linspace(0, total_distance, num=numpoints)
        b = bearing(lat0, lon0, lat1, lon1)

        for d in distance:
            p = VincentyDistance().destination(start, b, d)
            latitude.append(p.latitude)
            longitude.append(p.longitude)

    return list(map(np.array, [distance, latitude, longitude, b])) 

def get_distance_in_meter(self, node):
        return vincenty(self, node).meters 

def test_get_circle_centers():
    # test if circles fully cover the rect
    for sw, ne, w, h, r, circles in generate_testcases():
        # test with 1000 random points
        for tst in range(1000):
            # choose random point within rect
            p = (random.uniform(0,w), random.uniform(0,h))
            # translate to lat/lng
            pp = VincentyDistance(meters=p[0]).destination(
                VincentyDistance(meters=p[1])
                .destination(point=sw, bearing=90),
                bearing=0
            )
            # check if point is contained in any of the calculated circles
            assert any([vincenty(pp, Point(c[0], c[1])).meters <= r for c in circles]) 

def subset_wells_by_distance(self):
		# obtain longitude and latitudes from user
		while  len(self.userLocation) != 2:
			while True:
				try:
					self.userLocation = raw_input('\nDefine the center of your radius in Latitude (WGS84), and Longitude (WGS84) (separate by comma): ')
					self.userLocation = [x.strip() for x in self.userLocation.split(',')]
					self.userLocation = [float(x) for x in self.userLocation]
				except ValueError:
					print 'Please enter numbers'
					continue
				else:
					break

		# obtain the selection radius from user
		while True:
			try:
				userRadius = float(raw_input('\nDefine the radius within which you will keep all nearby wells (in miles): '))
			except ValueError:
				print 'Please enter numbers'
				continue
			else:
				break	

		# add vicintiy column to data set
		dist = np.zeros(len(self.wellDF['API/UWI']))
		for i,(lat,lon) in enumerate(zip(self.wellDF['Surface Latitude (WGS84)'], self.wellDF['Surface Longitude (WGS84)'])):
			dist[i] = vincenty([lat, lon], self.userLocation).miles

		self.wellDF['vicinity'] = dist

		# keep only wells withing the user selected radius
		self.wellDF = self.wellDF.loc[self.wellDF['vicinity'] <= userRadius]

		# notify user of changes to current selection
		print '%i wells selected' %(len(set(self.wellDF['API/UWI'])))

		return 

def gps_distance_elevation(fname):
    segment = gpxpy.parse(open(fname + '.gpx', 'r')).tracks[0].segments[0]
    elevation = []
    loc = []
    for p in segment.points:
        elevation.append(p.elevation)
        lat, lon = p.latitude, p.longitude
        loc.append((lat, lon))

    distance = np.array([0] + [vincenty(loc[i], loc[i-1]).meters for i in range(len(loc)-1)]).cumsum()
    plt.plot(distance, elevation, label=fname)
    plt.savefig(fname + '.png')
    plt.clf()
    return distance, elevation 

def main():
    """main"""
    args = get_args()
    infile = args.data
    sample_ids = re.split(r'\s*,\s*', args.sample_ids) if args.sample_ids else []
    print(sample_ids)
    
    if not os.path.isfile(infile):
        print('"{}" is not a file'.format(infile))
        sys.exit(1)

    records = []
    with open(infile) as fh:
        reader = csv.DictReader(fh, delimiter='\t')
        for rec in reader:
            records.append(rec)
    print('# rec = {}'.format(len(records)))

    combos = combinations(range(len(records)), 2)
    for i, j in combos:
        s1, s2 = records[i], records[j]
        dist = vincenty((s1['latitude'], s1['longitude']), 
                        (s2['latitude'], s2['longitude']))
        lat1, long1 = s1['latitude'], s1['longitude']
        print('{} -> {} = {}'.format(s1['sample_id'], s2['sample_id'], dist))
        print(s1)
        print(s2)



# -------------------------------------------------- 

def energy_generated(self):
        """
        The total energy of u_i located at l_i is:
        G(u_i,l_i) = -1 * sum of sij * g<u_i,u_j> over all friends.
        """
        #building dr, making ten bins based on social similarity
        distances_by_social_similarity = defaultdict(list)
        for user in self.users_with_location:
            location_user = self.mention_network.node_data(user)
            for neighbor in self.mention_network.neighbor_iter(user):
                if user == neighbor:
                    continue
                location_neighbor = self.mention_network.node_data(neighbor)
                social_similarity_rounded = round(self.sij[user][neighbor],1) #rounded to one significant figure
                distance = vincenty(location_user,location_neighbor)
                distances_by_social_similarity[social_similarity_rounded].append(distance)
        for social_similarity in distances_by_social_similarity:
            distances = distances_by_social_similarity[social_similarity]
            self.dr[social_similarity] = sum(distances)/len(distances)


        for user in self.users_with_location:
            location_user = self.mention_network.node_data(user)
            for neighbor in self.mention_network.neighbors_iter(user):
                if not user in self.users_with_location:
                    continue
                location_neighbor = self.mention_network.node_data(neighbor)

                social_similarity = self.sij[user][neighbor]
                #the exponent term, g<u_i,u_j> = -e^(-|l_i - l_j|/d_r)
                x = - vincenty(location_user,location_neighbor) / self.dr[round(social_similarity,1)]
                #summing sij * g<u_i,u_j> over all friends
                #I've factored out a -1 from np.exp(x) and cancelled it with
                #the leading -1 in the summation.
                self.g_ui[user] += social_similarity* np.exp(x) 

def training_build_rc_rg(self):
        """
        TODO
        """
        for user in self.users_with_location:
            rgs = []
            rcs = []
            for neighbor in self.mention_network.neighbors(user):
                if not neighbor in self.users_with_location:
                    continue
                location = self.mention_network.node_data(neighbor)
                rg = self.g_ui[neighbor]
                rc = self.c[neighbor]
                rgs.append((rg,location))
                rcs.append((rc,location))
            rgs = sorted(rgs)
            rcs = sorted(rcs)

            user_location = self.mention_network.node_data(user)
            best_distance = float('inf')
            for neighbor_location in rg:
                if vincenty(neighbor_location[1],user_location).km < best_distance:
                    self.best_neighbor[user] = neighbor_location[1]

            for rg in rgs:
                #storing the ranking index number r_G, and r_C, and the location
                self.rg_rc[user].append((rgs[1].index(rg[1]) + 1,rcs[1].index(rg[1]) + 1,rgs[1])) 

def reverse_geocode(self, lat, lon):
        """
        Returns the closest city name given a latitude and lonitude.
        """
        cities = self.reverse_geocoder[(round(lat,2),round(lon,2))]
        closest_location = float("inf")
        best_location = None
        for city in cities:
            lat2 = city[0]
            lon2 = city[1]
            distance = vincenty((lat,lon),(lat2,lon2)).km
            if distance < closest_location:
                closest_location = distance
                best_location = city[2]

        #if the direct search fails to find a location within 20km of the requested location,
        #the lat/lon is assumed to be either not available or a fringe case. In the latter
        #we check all the surrounding boxes in the reverse geocoder.
        if closest_location > 20:
            cities = self.reverse_geocoder[(round(lat+0.01,2),round(lon+0.01,2))] \
                 + self.reverse_geocoder[(round(lat+0.01,2),round(lon-0.01,2))] + self.reverse_geocoder[(round(lat-0.01,2),round(lon+0.01,2))] \
                 + self.reverse_geocoder[(round(lat-0.01,2),round(lon-0.01,2))] + self.reverse_geocoder[(round(lat,2),round(lon+0.01,2))] \
                 + self.reverse_geocoder[(round(lat,2),round(lon-0.01,2))] + self.reverse_geocoder[(round(lat+0.01,2),round(lon,2))] \
                 + self.reverse_geocoder[(round(lat-0.01,2),round(lon,2))]
            for city in cities:
                lat2 = city[0]
                lon2 = city[1]
                distance = vincenty((lat,lon),(lat2,lon2)).km
                if distance < closest_location:
                    closest_location = distance
                    best_location = city[2]

        return best_location 

def path_to_points(points, n=100, times=None):
    if times is None:
        times = [0] * len(points)

    if len(times) != len(points):
        if isinstance(times[0], datetime.datetime):
            times = times[0] + np.array([datetime.timedelta(0, d) for d in np.linspace(
                    0, (times[-1] - times[0]).total_seconds(), num=len(points))])
        else:
            times = np.linspace(times[0], times[-1], num=len(points))

    tuples = list(zip(points, points[1:], times, times[1:]))

    distance_between = []
    for pair in tuples:
        distance_between.append(vincenty(pair[0], pair[1]).km)

    total_distance = np.sum(distance_between)
    distance = []
    latitude = []
    longitude = []
    bearings = []
    output_time = []

    for index, pair in enumerate(tuples):
        n_pts = int(np.ceil(n * (distance_between[index] /
                                 total_distance)))
        n_pts = np.clip(n_pts, 2, n)
        p = list(map(geopy.Point, pair[0:2]))

        p_dist, p_lat, p_lon, b = points_between(p[0], p[1], n_pts)
        if len(distance) > 0:
            distance.extend(np.add(p_dist, distance[-1]))
        else:
            distance.extend(p_dist)

        latitude.extend(p_lat)
        longitude.extend(p_lon)
        bearings.extend([b] * len(p_dist))
        output_time.extend([
            pair[2] + i * (pair[3] - pair[2]) / (n_pts - 1)
            for i in range(0, n_pts)
        ])
    return distance, output_time, latitude, longitude, bearings