Geospatial Routing

Handle geographic coordinates and routing for logistics optimization.

Installation

pip install geopy networkx numpy

Quick Start

from geopy.distance import geodesic
from typing import Tuple, List

def calculate_distance(coord1: Tuple[float, float], coord2: Tuple[float, float]) -> float:
    """
    Calculate great-circle distance between two lat/lon coordinates.
    Returns distance in kilometers.
    """
    return geodesic(coord1, coord2).kilometers

def calculate_route_distance(coordinates: List[Tuple[float, float]]) -> float:
    """Calculate total distance of a route through coordinates."""
    total = 0
    for i in range(len(coordinates) - 1):
        total += calculate_distance(coordinates[i], coordinates[i + 1])
    return total

Common Patterns

Distance Matrix Construction

import numpy as np
from geopy.distance import geodesic

def build_distance_matrix(locations: dict) -> Tuple[np.ndarray, list]:
    """
    Build a distance matrix from location coordinates.
    locations: {location_id: (lat, lon)}
    Returns: (distance_matrix, location_ids)
    """
    location_ids = list(locations.keys())
    n = len(location_ids)
    matrix = np.zeros((n, n))

    for i, loc1 in enumerate(location_ids):
        for j, loc2 in enumerate(location_ids):
            if i != j:
                matrix[i][j] = geodesic(
                    locations[loc1],
                    locations[loc2]
                ).kilometers

    return matrix, location_ids


def build_time_matrix(distance_matrix: np.ndarray, speed_kmh: float = 50) -> np.ndarray:
    """Convert distance matrix to time matrix (in minutes)."""
    return (distance_matrix / speed_kmh) * 60

Geocoding and Reverse Geocoding

from geopy.geocoders import Nominatim

class GeocodingService:
    """Handle address to coordinate conversions."""

    def __init__(self, user_agent: str = "logistics_app"):
        self.geocoder = Nominatim(user_agent=user_agent)
        self.cache = {}

    def geocode(self, address: str) -> Tuple[float, float]:
        """Convert address to coordinates."""
        if address in self.cache:
            return self.cache[address]

        location = self.geocoder.geocode(address)
        if location:
            coords = (location.latitude, location.longitude)
            self.cache[address] = coords
            return coords
        return None

    def reverse_geocode(self, coords: Tuple[float, float]) -> str:
        """Convert coordinates to address."""
        location = self.geocoder.reverse(coords)
        return location.address if location else None

    def batch_geocode(self, addresses: List[str]) -> dict:
        """Geocode multiple addresses."""
        results = {}
        for address in addresses:
            coords = self.geocode(address)
            if coords:
                results[address] = coords
        return results

Routing Network

import networkx as nx

class RoutingNetwork:
    """Build and query a routing network."""

    def __init__(self):
        self.graph = nx.DiGraph()

    def add_road(self, from_node: str, to_node: str,
                 distance: float, time: float, bidirectional: bool = True):
        """Add a road segment to the network."""
        self.graph.add_edge(from_node, to_node, distance=distance, time=time)
        if bidirectional:
            self.graph.add_edge(to_node, from_node, distance=distance, time=time)

    def shortest_path(self, origin: str, destination: str,
                      weight: str = 'distance') -> dict:
        """Find shortest path between two points."""
        try:
            path = nx.shortest_path(self.graph, origin, destination, weight=weight)
            total_weight = nx.shortest_path_length(
                self.graph, origin, destination, weight=weight
            )
            return {
                'path': path,
                f'total_{weight}': total_weight,
                'segments': len(path) - 1
            }
        except nx.NetworkXNoPath:
            return None

    def find_nearest_node(self, coords: Tuple[float, float],
                          node_coords: dict) -> str:
        """Find the nearest network node to given coordinates."""
        min_dist = float('inf')
        nearest = None

        for node_id, node_coord in node_coords.items():
            if node_id in self.graph:
                dist = geodesic(coords, node_coord).kilometers
                if dist < min_dist:
                    min_dist = dist
                    nearest = node_id

        return nearest

Clustering Locations

from scipy.cluster.hierarchy import linkage, fcluster
from scipy.spatial.distance import pdist

def cluster_locations(locations: dict, max_clusters: int = None,
                      max_distance_km: float = None) -> dict:
    """
    Cluster locations geographically.
    Returns: {cluster_id: [location_ids]}
    """
    location_ids = list(locations.keys())
    coords = np.array([locations[loc] for loc in location_ids])

    # Calculate pairwise distances
    distances = pdist(coords, metric=lambda u, v: geodesic(u, v).kilometers)

    # Hierarchical clustering
    Z = linkage(distances, method='average')

    if max_clusters:
        clusters = fcluster(Z, max_clusters, criterion='maxclust')
    elif max_distance_km:
        clusters = fcluster(Z, max_distance_km, criterion='distance')
    else:
        clusters = fcluster(Z, 5, criterion='maxclust')

    # Group by cluster
    result = {}
    for loc_id, cluster_id in zip(location_ids, clusters):
        cluster_id = int(cluster_id)
        if cluster_id not in result:
            result[cluster_id] = []
        result[cluster_id].append(loc_id)

    return result

Zone-Based Routing

from shapely.geometry import Point, Polygon

class DeliveryZones:
    """Manage delivery zones for routing."""

    def __init__(self):
        self.zones = {}  # zone_id -> Polygon

    def add_zone(self, zone_id: str, polygon_coords: List[Tuple[float, float]]):
        """Add a delivery zone defined by polygon coordinates."""
        self.zones[zone_id] = Polygon(polygon_coords)

    def get_zone(self, coords: Tuple[float, float]) -> str:
        """Find which zone a coordinate belongs to."""
        point = Point(coords[1], coords[0])  # Note: shapely uses (lon, lat)
        for zone_id, polygon in self.zones.items():
            if polygon.contains(point):
                return zone_id
        return None

    def group_by_zone(self, locations: dict) -> dict:
        """Group locations by their delivery zone."""
        grouped = {}
        for loc_id, coords in locations.items():
            zone = self.get_zone(coords)
            if zone not in grouped:
                grouped[zone] = []
            grouped[zone].append(loc_id)
        return grouped

Route Optimization with Geography

def nearest_neighbor_route(depot: str, stops: List[str],
                           distance_matrix: np.ndarray,
                           location_index: dict) -> List[str]:
    """
    Build a route using nearest neighbor heuristic.
    location_index: {location_id: matrix_index}
    """
    route = [depot]
    remaining = set(stops)
    current = depot

    while remaining:
        current_idx = location_index[current]
        nearest = None
        min_dist = float('inf')

        for stop in remaining:
            stop_idx = location_index[stop]
            dist = distance_matrix[current_idx][stop_idx]
            if dist < min_dist:
                min_dist = dist
                nearest = stop

        route.append(nearest)
        remaining.remove(nearest)
        current = nearest

    route.append(depot)  # Return to depot
    return route


def two_opt_improvement(route: List[str], distance_matrix: np.ndarray,
                        location_index: dict) -> List[str]:
    """
    Improve route using 2-opt local search.
    """
    def route_distance(r):
        total = 0
        for i in range(len(r) - 1):
            total += distance_matrix[location_index[r[i]]][location_index[r[i+1]]]
        return total

    improved = True
    best_route = route[:]

    while improved:
        improved = False
        for i in range(1, len(best_route) - 2):
            for j in range(i + 1, len(best_route) - 1):
                new_route = (best_route[:i] +
                            best_route[i:j+1][::-1] +
                            best_route[j+1:])
                if route_distance(new_route) < route_distance(best_route):
                    best_route = new_route
                    improved = True

    return best_route