diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 115f3be6e..54fda4471 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -30,11 +30,11 @@ repos:
         types_or: [markdown, yaml]
 
   - repo: https://github.com/psf/black
-    rev: "23.7.0"
+    rev: "23.9.1"
     hooks:
       - id: black
 
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.0.287"
+    rev: "v0.0.292"
     hooks:
       - id: ruff
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4140666fc..b9909d759 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,9 @@
 
 ## Unreleased
 
+- refactor the distance module and add a new routing module (#1063)
+- move shortest_path and k_shortest_path functions to new routing module, with deprecation warning (#1063)
+- rename great_circle_vec and euclidean_dist_vec functions to great_circle and euclidean, with deprecation warning (#1063)
 - better automatic UTM handling in the projection module (#1059)
 - add to_latlong parameter to the projection.project_graph function for API consistency (#1057)
 - under-the-hood code clean-up (#1047)
diff --git a/osmnx/_api.py b/osmnx/_api.py
index a694c32c8..89abdc8d9 100644
--- a/osmnx/_api.py
+++ b/osmnx/_api.py
@@ -2,10 +2,8 @@
 
 from .bearing import add_edge_bearings
 from .bearing import orientation_entropy
-from .distance import k_shortest_paths
 from .distance import nearest_edges
 from .distance import nearest_nodes
-from .distance import shortest_path
 from .elevation import add_edge_grades
 from .elevation import add_node_elevations_google
 from .elevation import add_node_elevations_raster
@@ -44,6 +42,8 @@
 from .plot import plot_orientation
 from .projection import project_gdf
 from .projection import project_graph
+from .routing import k_shortest_paths
+from .routing import shortest_path
 from .simplification import consolidate_intersections
 from .simplification import simplify_graph
 from .speed import add_edge_speeds
diff --git a/osmnx/distance.py b/osmnx/distance.py
index 8a71cba25..f44b615f0 100644
--- a/osmnx/distance.py
+++ b/osmnx/distance.py
@@ -1,7 +1,5 @@
-"""Calculate distances and shortest paths and find nearest node/edge(s) to point(s)."""
+"""Calculate distances and find nearest node/edge(s) to point(s)."""
 
-import itertools
-import multiprocessing as mp
 from warnings import warn
 
 import networkx as nx
@@ -11,6 +9,7 @@
 from shapely.strtree import STRtree
 
 from . import projection
+from . import routing
 from . import utils
 from . import utils_geo
 from . import utils_graph
@@ -30,7 +29,7 @@
 EARTH_RADIUS_M = 6_371_009
 
 
-def great_circle_vec(lat1, lng1, lat2, lng2, earth_radius=EARTH_RADIUS_M):
+def great_circle(lat1, lng1, lat2, lng2, earth_radius=EARTH_RADIUS_M):
     """
     Calculate great-circle distances between pairs of points.
 
@@ -74,7 +73,7 @@ def great_circle_vec(lat1, lng1, lat2, lng2, earth_radius=EARTH_RADIUS_M):
     return arc * earth_radius
 
 
-def euclidean_dist_vec(y1, x1, y2, x2):
+def euclidean(y1, x1, y2, x2):
     """
     Calculate Euclidean distances between pairs of points.
 
@@ -102,6 +101,72 @@ def euclidean_dist_vec(y1, x1, y2, x2):
     return ((x1 - x2) ** 2 + (y1 - y2) ** 2) ** 0.5
 
 
+def great_circle_vec(lat1, lng1, lat2, lng2, earth_radius=EARTH_RADIUS_M):
+    """
+    Do not use, deprecated.
+
+    The `great_circle_vec` function has been renamed `great_circle`. Calling
+    `great_circle_vec` will raise an error in a future release.
+
+    Parameters
+    ----------
+    lat1 : float or numpy.array of float
+        first point's latitude coordinate
+    lng1 : float or numpy.array of float
+        first point's longitude coordinate
+    lat2 : float or numpy.array of float
+        second point's latitude coordinate
+    lng2 : float or numpy.array of float
+        second point's longitude coordinate
+    earth_radius : float
+        earth's radius in units in which distance will be returned (default is
+        meters)
+
+    Returns
+    -------
+    dist : float or numpy.array of float
+        distance from each (lat1, lng1) to each (lat2, lng2) in units of
+        earth_radius
+    """
+    warn(
+        "The `great_circle_vec` function has been renamed `great_circle`. Calling "
+        "`great_circle_vec` will raise an error in a future release.",
+        stacklevel=2,
+    )
+    return great_circle(lat1, lng1, lat2, lng2, earth_radius)
+
+
+def euclidean_dist_vec(y1, x1, y2, x2):
+    """
+    Do not use, deprecated.
+
+    The `euclidean_dist_vec` function has been renamed `euclidean`. Calling
+    `euclidean_dist_vec` will raise an error in a future release.
+
+    Parameters
+    ----------
+    y1 : float or numpy.array of float
+        first point's y coordinate
+    x1 : float or numpy.array of float
+        first point's x coordinate
+    y2 : float or numpy.array of float
+        second point's y coordinate
+    x2 : float or numpy.array of float
+        second point's x coordinate
+
+    Returns
+    -------
+    dist : float or numpy.array of float
+        distance from each (x1, y1) to each (x2, y2) in coordinates' units
+    """
+    warn(
+        "The `euclidean_dist_vec` function has been renamed `euclidean`. Calling "
+        "`euclidean_dist_vec` will raise an error in a future release.",
+        stacklevel=2,
+    )
+    return euclidean(y1, x1, y2, x2)
+
+
 def add_edge_lengths(G, precision=None, edges=None):
     """
     Add `length` attribute (in meters) to each edge.
@@ -163,7 +228,7 @@ def add_edge_lengths(G, precision=None, edges=None):
         raise ValueError(msg) from e
 
     # calculate great circle distances, round, and fill nulls with zeros
-    dists = great_circle_vec(c[:, 0], c[:, 1], c[:, 2], c[:, 3]).round(precision)
+    dists = great_circle(c[:, 0], c[:, 1], c[:, 2], c[:, 3]).round(precision)
     dists[np.isnan(dists)] = 0
     nx.set_edge_attributes(G, values=dict(zip(uvk, dists)), name="length")
 
@@ -355,51 +420,12 @@ def nearest_edges(G, X, Y, interpolate=None, return_dist=False):
     return ne
 
 
-def _single_shortest_path(G, orig, dest, weight):
-    """
-    Solve the shortest path from an origin node to a destination node.
-
-    This function is a convenience wrapper around networkx.shortest_path, with
-    exception handling for unsolvable paths. It uses Dijkstra's algorithm.
-
-    Parameters
-    ----------
-    G : networkx.MultiDiGraph
-        input graph
-    orig : int
-        origin node ID
-    dest : int
-        destination node ID
-    weight : string
-        edge attribute to minimize when solving shortest path
-
-    Returns
-    -------
-    path : list
-        list of node IDs constituting the shortest path
-    """
-    try:
-        return nx.shortest_path(G, orig, dest, weight=weight, method="dijkstra")
-    except nx.exception.NetworkXNoPath:  # pragma: no cover
-        utils.log(f"Cannot solve path from {orig} to {dest}")
-        return None
-
-
 def shortest_path(G, orig, dest, weight="length", cpus=1):
     """
-    Solve shortest path from origin node(s) to destination node(s).
-
-    Uses Dijkstra's algorithm. If `orig` and `dest` are single node IDs, this
-    will return a list of the nodes constituting the shortest path between
-    them. If `orig` and `dest` are lists of node IDs, this will return a list
-    of lists of the nodes constituting the shortest path between each
-    origin-destination pair. If a path cannot be solved, this will return None
-    for that path. You can parallelize solving multiple paths with the `cpus`
-    parameter, but be careful to not exceed your available RAM.
+    Do not use, deprecated.
 
-    See also `k_shortest_paths` to solve multiple shortest paths between a
-    single origin and destination. For additional functionality or different
-    solver algorithms, use NetworkX directly.
+    The `shortest_path` function has moved to the `routing` module. Calling
+    it via the `distance` module will raise an error in a future release.
 
     Parameters
     ----------
@@ -420,49 +446,20 @@ def shortest_path(G, orig, dest, weight="length", cpus=1):
         list of node IDs constituting the shortest path, or, if orig and dest
         are lists, then a list of path lists
     """
-    _verify_edge_attribute(G, weight)
-
-    # if neither orig nor dest is iterable, just return the shortest path
-    if not (hasattr(orig, "__iter__") or hasattr(dest, "__iter__")):
-        return _single_shortest_path(G, orig, dest, weight)
-
-    # if both orig and dest are iterables, ensure they have same lengths
-    if hasattr(orig, "__iter__") and hasattr(dest, "__iter__"):
-        if len(orig) != len(dest):  # pragma: no cover
-            msg = "orig and dest must contain same number of elements"
-            raise ValueError(msg)
-
-        if cpus is None:
-            cpus = mp.cpu_count()
-        cpus = min(cpus, mp.cpu_count())
-        utils.log(f"Solving {len(orig)} paths with {cpus} CPUs...")
-
-        # if single-threading, calculate each shortest path one at a time
-        if cpus == 1:
-            paths = [_single_shortest_path(G, o, d, weight) for o, d in zip(orig, dest)]
-
-        # if multi-threading, calculate shortest paths in parallel
-        else:
-            args = ((G, o, d, weight) for o, d in zip(orig, dest))
-            pool = mp.Pool(cpus)
-            sma = pool.starmap_async(_single_shortest_path, args)
-            paths = sma.get()
-            pool.close()
-            pool.join()
-
-        return paths
-
-    # otherwise only one of orig or dest is iterable and the other is not
-    msg = "orig and dest must either both be iterable or neither must be iterable"
-    raise ValueError(msg)
+    warn(
+        "The `shortest_path` function has moved to the `routing` module. "
+        "Calling it via the `distance` module will raise an error in a future release.",
+        stacklevel=2,
+    )
+    return routing.shortest_path(G, orig, dest, weight, cpus)
 
 
 def k_shortest_paths(G, orig, dest, k, weight="length"):
     """
-    Solve `k` shortest paths from an origin node to a destination node.
+    Do not use, deprecated.
 
-    Uses Yen's algorithm. See also `shortest_path` to solve just the one
-    shortest path.
+    The `k_shortest_paths` function has moved to the `routing` module. Calling
+    it via the `distance` module will raise an error in a future release.
 
     Parameters
     ----------
@@ -484,34 +481,9 @@ def k_shortest_paths(G, orig, dest, k, weight="length"):
         a generator of `k` shortest paths ordered by total weight. each path
         is a list of node IDs.
     """
-    _verify_edge_attribute(G, weight)
-    paths_gen = nx.shortest_simple_paths(utils_graph.get_digraph(G, weight), orig, dest, weight)
-    yield from itertools.islice(paths_gen, 0, k)
-
-
-def _verify_edge_attribute(G, attr):
-    """
-    Verify attribute values are numeric and non-null across graph edges.
-
-    Raises a `ValueError` if attribute contains non-numeric values and raises
-    a warning if attribute is missing or null on any edges.
-
-    Parameters
-    ----------
-    G : networkx.MultiDiGraph
-        input graph
-    attr : string
-        edge attribute to verify
-
-    Returns
-    -------
-    None
-    """
-    try:
-        values = np.array(tuple(G.edges(data=attr)))[:, 2]
-        values_float = values.astype(float)
-        if np.isnan(values_float).any():
-            warn(f"The attribute {attr!r} is missing or null on some edges.", stacklevel=2)
-    except ValueError as e:
-        msg = f"The edge attribute {attr!r} contains non-numeric values."
-        raise ValueError(msg) from e
+    warn(
+        "The `k_shortest_paths` function has moved to the `routing` module. "
+        "Calling it via the `distance` module will raise an error in a future release.",
+        stacklevel=2,
+    )
+    return routing.k_shortest_paths(G, orig, dest, k, weight)
diff --git a/osmnx/osm_xml.py b/osmnx/osm_xml.py
index 91c78b538..b51338f82 100644
--- a/osmnx/osm_xml.py
+++ b/osmnx/osm_xml.py
@@ -153,7 +153,7 @@ def save_graph_xml(
     None
     """
     warn(
-        "The save_graph_xml has moved from the osm_xml module to the io module. "
+        "The save_graph_xml function has moved from the osm_xml module to the io module. "
         " osm_xml.save_graph_xml has been deprecated and will be removed in a "
         " future release. Access the function via the io module instead.",
         stacklevel=2,
diff --git a/osmnx/routing.py b/osmnx/routing.py
new file mode 100644
index 000000000..3b23ac974
--- /dev/null
+++ b/osmnx/routing.py
@@ -0,0 +1,174 @@
+"""Calculate weighted shortest paths between graph nodes."""
+
+import itertools
+import multiprocessing as mp
+from warnings import warn
+
+import networkx as nx
+import numpy as np
+
+from . import utils
+from . import utils_graph
+
+
+def shortest_path(G, orig, dest, weight="length", cpus=1):
+    """
+    Solve shortest path from origin node(s) to destination node(s).
+
+    Uses Dijkstra's algorithm. If `orig` and `dest` are single node IDs, this
+    will return a list of the nodes constituting the shortest path between
+    them. If `orig` and `dest` are lists of node IDs, this will return a list
+    of lists of the nodes constituting the shortest path between each
+    origin-destination pair. If a path cannot be solved, this will return None
+    for that path. You can parallelize solving multiple paths with the `cpus`
+    parameter, but be careful to not exceed your available RAM.
+
+    See also `k_shortest_paths` to solve multiple shortest paths between a
+    single origin and destination. For additional functionality or different
+    solver algorithms, use NetworkX directly.
+
+    Parameters
+    ----------
+    G : networkx.MultiDiGraph
+        input graph
+    orig : int or list
+        origin node ID, or a list of origin node IDs
+    dest : int or list
+        destination node ID, or a list of destination node IDs
+    weight : string
+        edge attribute to minimize when solving shortest path
+    cpus : int
+        how many CPU cores to use; if None, use all available
+
+    Returns
+    -------
+    path : list
+        list of node IDs constituting the shortest path, or, if orig and dest
+        are lists, then a list of path lists
+    """
+    _verify_edge_attribute(G, weight)
+
+    # if neither orig nor dest is iterable, just return the shortest path
+    if not (hasattr(orig, "__iter__") or hasattr(dest, "__iter__")):
+        return _single_shortest_path(G, orig, dest, weight)
+
+    # if only 1 of orig or dest is iterable and the other is not, raise error
+    if not (hasattr(orig, "__iter__") and hasattr(dest, "__iter__")):
+        msg = "orig and dest must either both be iterable or neither must be iterable"
+        raise ValueError(msg)
+
+    # if both orig and dest are iterable, ensure they have same lengths
+    if len(orig) != len(dest):  # pragma: no cover
+        msg = "orig and dest must be of equal length"
+        raise ValueError(msg)
+
+    # determine how many cpu cores to use
+    if cpus is None:
+        cpus = mp.cpu_count()
+    cpus = min(cpus, mp.cpu_count())
+    utils.log(f"Solving {len(orig)} paths with {cpus} CPUs...")
+
+    # if single-threading, calculate each shortest path one at a time
+    if cpus == 1:
+        paths = [_single_shortest_path(G, o, d, weight) for o, d in zip(orig, dest)]
+
+    # if multi-threading, calculate shortest paths in parallel
+    else:
+        args = ((G, o, d, weight) for o, d in zip(orig, dest))
+        pool = mp.Pool(cpus)
+        sma = pool.starmap_async(_single_shortest_path, args)
+        paths = sma.get()
+        pool.close()
+        pool.join()
+
+    return paths
+
+
+def k_shortest_paths(G, orig, dest, k, weight="length"):
+    """
+    Solve `k` shortest paths from an origin node to a destination node.
+
+    Uses Yen's algorithm. See also `shortest_path` to solve just the one
+    shortest path.
+
+    Parameters
+    ----------
+    G : networkx.MultiDiGraph
+        input graph
+    orig : int
+        origin node ID
+    dest : int
+        destination node ID
+    k : int
+        number of shortest paths to solve
+    weight : string
+        edge attribute to minimize when solving shortest paths. default is
+        edge length in meters.
+
+    Yields
+    ------
+    path : list
+        a generator of `k` shortest paths ordered by total weight. each path
+        is a list of node IDs.
+    """
+    _verify_edge_attribute(G, weight)
+    paths_gen = nx.shortest_simple_paths(utils_graph.get_digraph(G, weight), orig, dest, weight)
+    yield from itertools.islice(paths_gen, 0, k)
+
+
+def _single_shortest_path(G, orig, dest, weight):
+    """
+    Solve the shortest path from an origin node to a destination node.
+
+    This function is a convenience wrapper around networkx.shortest_path, with
+    exception handling for unsolvable paths. It uses Dijkstra's algorithm.
+
+    Parameters
+    ----------
+    G : networkx.MultiDiGraph
+        input graph
+    orig : int
+        origin node ID
+    dest : int
+        destination node ID
+    weight : string
+        edge attribute to minimize when solving shortest path
+
+    Returns
+    -------
+    path : list
+        list of node IDs constituting the shortest path
+    """
+    try:
+        return nx.shortest_path(G, orig, dest, weight=weight, method="dijkstra")
+    except nx.exception.NetworkXNoPath:  # pragma: no cover
+        utils.log(f"Cannot solve path from {orig} to {dest}")
+        return None
+
+
+def _verify_edge_attribute(G, attr):
+    """
+    Verify attribute values are numeric and non-null across graph edges.
+
+    Raises a `ValueError` if attribute contains non-numeric values and raises
+    a warning if attribute is missing or null on any edges.
+
+    Parameters
+    ----------
+    G : networkx.MultiDiGraph
+        input graph
+    attr : string
+        edge attribute to verify
+
+    Returns
+    -------
+    None
+    """
+    try:
+        values = np.array(tuple(G.edges(data=attr)))[:, 2]
+        values_float = values.astype(float)
+        if np.isnan(values_float).any():
+            warn(f"The attribute {attr!r} is missing or null on some edges.", stacklevel=2)
+    except ValueError as e:
+        msg = f"The edge attribute {attr!r} contains non-numeric values."
+        raise ValueError(msg) from e
diff --git a/osmnx/stats.py b/osmnx/stats.py
index 7fa05efdf..0dd3433f8 100644
--- a/osmnx/stats.py
+++ b/osmnx/stats.py
@@ -244,9 +244,9 @@ def circuity_avg(Gu):
     # calculate straight-line distances as euclidean distances if projected or
     # great-circle distances if unprojected
     if projection.is_projected(Gu.graph["crs"]):
-        sl_dists = distance.euclidean_dist_vec(y1=y1, x1=x1, y2=y2, x2=x2)
+        sl_dists = distance.euclidean(y1=y1, x1=x1, y2=y2, x2=x2)
     else:
-        sl_dists = distance.great_circle_vec(lat1=y1, lng1=x1, lat2=y2, lng2=x2)
+        sl_dists = distance.great_circle(lat1=y1, lng1=x1, lat2=y2, lng2=x2)
 
     # return the ratio, handling possible division by zero
     sl_dists_total = sl_dists[~np.isnan(sl_dists)].sum()
diff --git a/osmnx/utils.py b/osmnx/utils.py
index 11bf5251a..c1de56492 100644
--- a/osmnx/utils.py
+++ b/osmnx/utils.py
@@ -4,7 +4,7 @@
 import logging as lg
 import os
 import sys
-import unicodedata
+import unicodedata as ud
 from contextlib import redirect_stdout
 from pathlib import Path
 from warnings import warn
@@ -300,10 +300,8 @@ def log(message, level=None, name=None, filename=None):
         # prepend timestamp
         message = f"{ts()} {message}"
 
-        # convert to ascii so it doesn't break windows terminals
-        message = (
-            unicodedata.normalize("NFKD", str(message)).encode("ascii", errors="replace").decode()
-        )
+        # convert to ascii so it works in windows command prompts
+        message = ud.normalize("NFKD", message).encode("ascii", errors="replace").decode()
 
         # print explicitly to terminal in case jupyter notebook is the stdout
         if getattr(sys.stdout, "_original_stdstream_copy", None) is not None:
diff --git a/tests/test_osmnx.py b/tests/test_osmnx.py
index d24247d62..6b30a7384 100755
--- a/tests/test_osmnx.py
+++ b/tests/test_osmnx.py
@@ -1,6 +1,6 @@
 """Test suite for the package."""
 
-# use agg backend so you don't need a display on ci
+# use agg backend so you don't need a display on CI
 # do this first before pyplot is imported by anything
 import matplotlib as mpl
 
@@ -280,6 +280,7 @@ def test_routing():
     route = ox.shortest_path(G, orig_node, dest_node, weight="time")
     # test good weight
     route = ox.shortest_path(G, orig_node, dest_node, weight="travel_time")
+    route = ox.distance.shortest_path(G, orig_node, dest_node, weight="travel_time")
 
     route_edges = ox.utils_graph.route_to_gdf(G, route, "travel_time")
     attributes = ox.utils_graph.get_route_edge_attributes(G, route)
@@ -299,8 +300,13 @@ def test_routing():
 
     # test k shortest paths
     routes = ox.k_shortest_paths(G, orig_node, dest_node, k=2, weight="travel_time")
+    routes = ox.distance.k_shortest_paths(G, orig_node, dest_node, k=2, weight="travel_time")
     fig, ax = ox.plot_graph_routes(G, list(routes))
 
+    # test great circle and euclidean distance calculators
+    assert ox.distance.great_circle_vec(0, 0, 1, 1) == pytest.approx(157249.6034105)
+    assert ox.distance.euclidean_dist_vec(0, 0, 1, 1) == pytest.approx(1.4142135)
+
     # test folium with keyword arguments to pass to folium.PolyLine
     gm = ox.plot_graph_folium(G, popup_attribute="name", color="#333333", weight=5, opacity=0.7)
     rm = ox.plot_route_folium(G, route, color="#cc0000", weight=5, opacity=0.7)