def nearest_nodes(
    G: nx.MultiDiGraph,
    X: float | Iterable[float],
    Y: float | Iterable[float],
    *,
    return_dist: bool = False,
) -> (
    int
    | npt.NDArray[np.int64]
    | tuple[int, float]
    | tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]
):
    """
    Find the nearest node to a point or to each of several points.

    If `X` and `Y` are single coordinate values, this function will return the
    nearest node to that point. If `X` and `Y` are iterables of coordinate
    values, it will return the nearest node to each point.

    This function is vectorized: if you have many points to search for, pass
    them in one call as numpy arrays (avoid using loops) to maximize runtime
    speed. If the graph is projected, it uses a k-d tree for Euclidean nearest
    neighbor search, which requires that scipy is installed as an optional
    dependency. If the graph is unprojected, it uses a ball tree for haversine
    nearest neighbor search, which requires that scikit-learn is installed as
    an optional dependency.

    Parameters
    ----------
    G
        Graph in which to find nearest nodes.
    X
        The points' x (longitude) coordinates, in same CRS/units as graph and
        containing no nulls.
    Y
        The points' y (latitude) coordinates, in same CRS/units as graph and
        containing no nulls.
    return_dist
        If True, optionally also return the distance(s) between point(s) and
        nearest node(s).

    Returns
    -------
    nn or (nn, dist)
        Nearest node ID(s) or optionally a tuple of ID(s) and distance(s)
        between each point and its nearest node.
    """
    # make coordinates arrays whether user passed iterable values or not
    if not (isinstance(X, Iterable) and isinstance(Y, Iterable)):
        is_scalar = True
        X_arr = np.array([X])
        Y_arr = np.array([Y])
    else:
        is_scalar = False
        X_arr = np.array(X)
        Y_arr = np.array(Y)

    if np.isnan(X_arr).any() or np.isnan(Y_arr).any():  # pragma: no cover
        msg = "`X` and `Y` cannot contain nulls."
        raise ValueError(msg)

    nodes = convert.graph_to_gdfs(G, edges=False, node_geometry=False)[["x", "y"]]
    nn_array: npt.NDArray[np.int64]
    dist_array: npt.NDArray[np.float64]

    if projection.is_projected(G.graph["crs"]):
        # if projected, use k-d tree for euclidean nearest-neighbor search
        if cKDTree is None:  # pragma: no cover
            msg = "scipy must be installed as an optional dependency to search a projected graph."
            raise ImportError(msg)
        dist_array, pos = cKDTree(nodes).query(np.array([X_arr, Y_arr]).T, k=1)
        nn_array = nodes.index[pos].to_numpy()

    else:
        # if unprojected, use ball tree for haversine nearest-neighbor search
        if BallTree is None:  # pragma: no cover
            msg = "scikit-learn must be installed as an optional dependency to search an unprojected graph."
            raise ImportError(msg)
        # haversine requires lat, lon coords in radians
        nodes_rad = np.deg2rad(nodes[["y", "x"]])
        points_rad = np.deg2rad(np.array([Y_arr, X_arr]).T)
        dist_array, pos = BallTree(nodes_rad, metric="haversine").query(points_rad, k=1)
        dist_array = dist_array[:, 0] * EARTH_RADIUS_M  # convert radians -> meters
        nn_array = nodes.index[pos[:, 0]].to_numpy()

    # convert results to correct types for return
    if is_scalar:
        nn = int(nn_array[0])
        dist = float(dist_array[0])
        if return_dist:
            return nn, dist
        # otherwise
        return nn

    # otherwise
    if return_dist:
        return nn_array, dist_array
    # otherwise
    return nn_array

    Find the nearest node to a point or to each of several points.

    If `X` and `Y` are single coordinate values, this function will return the
    nearest node to that point. If `X` and `Y` are iterables of coordinate
    values, it will return the nearest node to each point.

    This function is vectorized: if you have many points to search for, pass
    them in one call as numpy arrays (avoid using loops) to maximize runtime
    speed. If the graph is projected, it uses a k-d tree for Euclidean nearest
    neighbor search, which requires that scipy is installed as an optional
    dependency. If the graph is unprojected, it uses a ball tree for haversine
    nearest neighbor search, which requires that scikit-learn is installed as
    an optional dependency.

    Parameters
    ----------
    G
        Graph in which to find nearest nodes.
    X
        The points' x (longitude) coordinates, in same CRS/units as graph and
        containing no nulls.
    Y
        The points' y (latitude) coordinates, in same CRS/units as graph and
        containing no nulls.
    return_dist
        If True, optionally also return the distance(s) between point(s) and
        nearest node(s).

    Returns
    -------
    nn or (nn, dist)
        Nearest node ID(s) or optionally a tuple of ID(s) and distance(s)
        between each point and its nearest node.
    
osmnx/distance.py