How to use the geopy.distance.VincentyDistance function in geopy

Parameters
    ----------
    start_point: geopy.Point
    end_point: geopy.Point
    distance_meters: float

    Returns
    -------
    point: geopy.Point
        A new point between the start and the end points.
    """
    bearing = get_bearing(start_point, end_point)

    distance_km = distance_meters / 1000
    d = geo_dist.VincentyDistance(kilometers=distance_km)
    destination = d.destination(point=start_point, bearing=bearing)

    return geopy.Point(destination.latitude, destination.longitude)

def partition_edge(edge, distance_interval):
	"""
	given an edge, creates holes every x meters (distance_interval)
	:param edge: a given edge
	:param distance_interval: in meters
	:return: list of holes
	"""

	# We always return the source node of the edge, hopefully the target will be added as the source of another edge.
	holes = []
	d = geopy.distance.VincentyDistance(meters=distance_interval)
	# make sure we are using lat,lon not lon,lat as a reference.
	startpoint = geopy.Point(edge[0].get_coordinates())
	endpoint = geopy.Point(edge[1].get_coordinates())
	initial_dist = geopy.distance.distance(startpoint, endpoint).meters
	if initial_dist < distance_interval:
		# return [], distance_interval - initial_dist
		return holes
	# compute the angle=bearing at which we need to be moving.
	bearing = calculate_bearing(startpoint[0], startpoint[1], endpoint[0], endpoint[1])
	last_point = startpoint
	diff_time = edge[1].last_seen - edge[0].last_seen
	delta_time = diff_time.days*24*3600 + diff_time.seconds
	time_increment = delta_time / (int(initial_dist) / distance_interval)
	for i in range(int(initial_dist) / distance_interval):
		new_point = geopy.Point(d.destination(point=last_point, bearing=bearing))
		str_timestamp = datetime.datetime.strftime(edge[0].last_seen + datetime.timedelta(seconds=time_increment), "%Y-%m-%d %H:%M:%S+03")

var = dataset.variables[variable]

            if depth == 'bottom':
                bdata = var[time[0], :, miny:maxy, minx:maxx]
                blats = dataset.variables[latvar.name][miny:maxy, minx:maxx]
                blons = dataset.variables[lonvar.name][miny:maxy, minx:maxx]

                reshaped = bdata.reshape([bdata.shape[0], -1])
                edges = np.array(np.ma.notmasked_edges(reshaped, axis=0))

                depths = edges[1, 0, :]
                indices = edges[1, 1, :]

                from geopy.distance import VincentyDistance
                distance = VincentyDistance()

                weighted = 0
                totalweight = 0
                for i, index in enumerate(indices):
                    dist = distance.measure((lat, lon),
                                            (blats.ravel()[index],
                                             blons.ravel()[index]))
                    weight = 1 / (dist ** 2)
                    totalweight += weight
                    weighted = depths[i] * weight
                weighted /= totalweight

                depth = int(np.round(weighted))

            depthall = False
            if 'deptht' in var.dimensions or 'depth' in var.dimensions:

def new_gps_coord_after_offset_inBodyFrame(original_gps_coord, displacement, current_heading, rotation_degree_relative):
    # current_heading is in degree, North = 0, East = 90.
    # Get rotation degree in local frame.
    rotation_degree_absolute = rotation_degree_relative + current_heading
    if rotation_degree_absolute >= 360:
        rotation_degree_absolute -= 360
    vincentyDistance = geopy.distance.VincentyDistance(meters = displacement)
    original_point = geopy.Point(original_gps_coord[0], original_gps_coord[1])
    new_gps_coord = vincentyDistance.destination(point=original_point, bearing=rotation_degree_absolute)
    new_gps_lat = new_gps_coord.latitude
    new_gps_lon = new_gps_coord.longitude
    # If convert float to decimal, round will be accurate, but will take 50% more time. Not necessary.
    #new_gps_lat = decimal.Decimal(new_gps_lat)
    #new_gps_lon = decimal.Decimal(new_gps_lon)
    return (round(new_gps_lat, 7), round(new_gps_lon, 7))

: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

geometry = [Point(org.longitude, org.latitude)]
                if all(value is not None for value in (
                        horizontal_uncertainty_max, horizontal_uncertainty_min,
                        horizontal_uncertainty_max_azimuth)):
                    import geopy
                    import geopy.distance
                    start = geopy.Point(latitude=org.latitude,
                                        longitude=org.longitude)
                    lines = []
                    for distance, azimuth in (
                            (horizontal_uncertainty_max,
                             horizontal_uncertainty_max_azimuth),
                            (horizontal_uncertainty_min,
                             horizontal_uncertainty_max_azimuth + 90)):
                        azimuth = azimuth % 180
                        distance = geopy.distance.VincentyDistance(
                            kilometers=distance / 1e3)
                        end1 = distance.destination(
                            point=start, bearing=azimuth)
                        end2 = distance.destination(
                            point=start, bearing=azimuth + 180)
                        line = LineString((end1.longitude, end1.latitude),
                                          (org.longitude, org.latitude),
                                          (end2.longitude, end2.latitude))
                        lines.append(line)
                    geometry.append(MultiLineString(lines))
                else:
                    geometry.append(MultiLineString([]))

            indices.append({
                "quakeml_id": str(event.resource_id),
                "latitude": org.latitude if org else None,

def new_gps_coord_after_offset_inLocalFrame(original_gps_coord, displacement, rotation_degree):
    vincentyDistance = geopy.distance.VincentyDistance(meters = displacement)
    original_point = geopy.Point(original_gps_coord[0], original_gps_coord[1])
    new_gps_coord = vincentyDistance.destination(point=original_point, bearing=rotation_degree)
    new_gps_lat = new_gps_coord.latitude
    new_gps_lon = new_gps_coord.longitude
    # If convert float to decimal, round will be accurate, but will take 50% more time. Not necessary.
    #new_gps_lat = decimal.Decimal(new_gps_lat)
    #new_gps_lon = decimal.Decimal(new_gps_lon)
    return (round(new_gps_lat, 7), round(new_gps_lon, 7))

def _path_to_points(points, n, intimes=None):
    if intimes is None:
        intimes = [0] * len(points)

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

    d = VincentyDistance()
    dist_between_pts = []
    for pair in tuples:
        dist_between_pts.append(d.measure(pair[0], pair[1]))

    total_distance = np.sum(dist_between_pts)
    distances = []
    target_lat = []
    target_lon = []
    bearings = []
    times = []
    for idx, tup in enumerate(tuples):
        npts = int(np.ceil(n * (dist_between_pts[idx] /
                                total_distance)))
        if npts < 2:
            npts = 2
        p0 = geopy.Point(tup[0])

def get_lat_lon(self):
        """
            The only thing I know about this function is that it gets the latitudes and longitudes

            I know it deals with weird edge cases etc. but
            honestly it looks like gibberish
        """
        distance = VincentyDistance()
        height = distance.measure(
            (self.bounds[0], self.centroid[1]),
            (self.bounds[2], self.centroid[1])
        ) * 1000 * 1.25
        width = distance.measure(
            (self.centroid[0], self.bounds[1]),
            (self.centroid[0], self.bounds[3])
        ) * 1000 * 1.25

        if self.projection == 'EPSG:32661':             # north pole projection
            near_pole, covers_pole = self.pole_proximity(self.points[0])
            blat = min(self.bounds[0], self.bounds[2])
            blat = 5 * np.floor(blat / 5)

            if self.centroid[0] > 80 or near_pole or covers_pole:
                self.basemap = basemap.load_map(