A* Pathfinding for Maritime Weather Routing

Adaptation of A* graph search to minimum-fuel ocean routing with physics-based cost evaluation and time-varying weather

Abstract

The A* algorithm, originally developed for shortest-path search in weighted graphs, is adapted here for minimum-fuel ocean routing subject to time-varying meteorological and oceanographic conditions. WindMar discretises the ocean surface into a regular latitude-longitude grid at configurable resolution (default 0.5°, approximately 30 nautical miles at the equator) with eight-connected cell adjacency. Land cells are excluded using a high-resolution global land mask, and the search corridor is bounded by a configurable margin around the great-circle path between origin and destination. The cost function at each edge evaluates fuel consumption via the full hydrodynamic vessel model, incorporating wind resistance, wave-added resistance, and ocean current effects on speed over ground. An admissible heuristic based on great-circle distance, maximum attainable speed, and a conservative fuel rate ensures that the search terminates with a provably optimal path. Post-processing applies Douglas-Peucker simplification with land-crossing checks to reduce waypoint density, followed by optional variable-speed optimisation across the smoothed route. The complete pipeline — grid construction, A* search, smoothing, and route statistics — executes in 0.5–2 seconds for typical Mediterranean routes, enabling interactive use in the WindMar web application.

1. Introduction

The problem of finding an optimal route for a vessel across an ocean subject to environmental forcing has attracted sustained attention since Dijkstra (1959) established the theoretical foundations for shortest-path search in weighted graphs. Hart, Nilsson, and Raphael (1968) extended Dijkstra's algorithm with a heuristic estimate of the remaining cost to the goal, producing the A* algorithm that retains optimality guarantees while substantially reducing the number of nodes explored. Hagiwara (1989) was among the first to apply graph-search methods to maritime weather routing, constructing isochrone-based routing charts that account for wave-induced speed loss. More recently, Szlapczynska (2015) demonstrated multi-objective weather routing using evolutionary algorithms, and Vettor and Guedes Soares (2016) developed a complete ship weather routing system that integrates vessel performance models with numerical weather prediction data.

WindMar adopts the A* algorithm rather than dynamic programming or evolutionary approaches for several reasons. First, A* with an admissible heuristic guarantees that the returned path is optimal with respect to the cost function, a property that evolutionary algorithms cannot provide without exhaustive evaluation. Second, the heuristic guides the search toward the destination, dramatically reducing the number of cells explored compared to Dijkstra's algorithm, which expands uniformly in all directions. Third, the algorithm's node-expansion structure maps naturally onto a regular geographic grid, where each cell represents a navigable ocean area and each edge represents a candidate leg with computable fuel cost. The physics-based cost evaluation at each edge draws on the hydrodynamic vessel model for resistance and fuel calculations and on the weather data pipeline for time-varying environmental conditions at each grid point.

This article is organised as follows. Section 2 formulates the optimisation problem. Section 3 describes the construction of the navigable ocean grid. Section 4 presents the A* search algorithm and its adaptation to fuel-minimisation routing. Section 5 details the per-leg voyage calculation that underpins the edge cost function. Section 6 describes the Douglas-Peucker route smoothing procedure. Section 7 discusses temporal weather integration. Section 8 analyses the computational performance of the pipeline.

2. Problem Formulation

The route optimisation problem solved by WindMar can be stated as follows. Given an origin waypoint, a destination waypoint, a departure time, a vessel model with known hydrodynamic characteristics, and a time-varying weather field covering the transit area, find the sequence of intermediate waypoints and (optionally) the speed at each leg that minimises total fuel consumption while satisfying all safety and geographic constraints. The primary objective is fuel minimisation; however, the cost function can incorporate weighted time penalties to produce fuel-time trade-off solutions.

The decision variables are the geographic positions of intermediate waypoints (constrained to lie on ocean grid cells) and, when variable-speed optimisation is enabled, the vessel speed at each leg. The fuel consumption at each leg is not a simple function of distance; it depends on the wind, wave, and current conditions encountered, the vessel's heading relative to the environmental forcing, and the chosen speed. This coupling between route geometry and environmental conditions makes the problem inherently non-linear and motivates the use of graph search with physics-based edge cost evaluation rather than analytic optimisation.

Safety constraints are enforced through a multiplicative cost factor. The safety_constraints.get_safety_cost_factor() method evaluates seakeeping criteria (roll angle, pitch angle, vertical acceleration) at each candidate leg and returns a factor of 1.0 for safe conditions, a value greater than 1.0 for marginal conditions (penalising but not prohibiting the leg), and infinity for dangerous conditions (effectively blocking the leg from the search). Geographic constraints are enforced by excluding land cells from the grid and applying zone penalties via zone_checker.get_path_penalty(), which returns a multiplicative factor for legs that traverse restricted or undesirable areas.

3. Graph Construction

The search graph is constructed as a regular latitude-longitude grid covering a corridor around the great-circle path between origin and destination. The RouteOptimizer._build_grid() method creates this grid with the following procedure. First, the latitude and longitude extents of the origin and destination are computed, and a configurable margin (default margin_deg=5.0) is added on all sides to allow the search to explore routes that deviate from the direct path. The grid is then populated with cells at the configured resolution, where each cell is represented by a GridCell dataclass containing its latitude, longitude, row index, and column index. The GridCell is hashable, enabling efficient use in Python sets and dictionaries for the closed set and g-score lookup table.

Land filtering is performed by querying the global_land_mask library at each grid point. Any cell classified as land is excluded from the grid dictionary, producing an ocean-only navigable graph. This filtering is critical for preventing the A* search from generating routes that cross landmasses. The land mask operates at approximately 0.05° resolution, which is finer than the default routing grid resolution of 0.5°, ensuring that narrow straits and coastal passages are correctly classified. Further details on the ocean masking procedure are provided in the Weather Fields article.

3.1 Grid Connectivity

Each cell is connected to its eight immediate neighbours in a regular grid topology: north, northeast, east, southeast, south, southwest, west, and northwest. The eight direction vectors are defined as the constant DIRECTIONS = [(-1,0), (-1,1), (0,1), (1,1), (1,0), (1,-1), (0,-1), (-1,-1)], where tuple elements represent (row_delta, col_delta). Cardinal moves (N, E, S, W) traverse one grid cell spacing, while diagonal moves (NE, SE, SW, NW) traverse a distance of √2 times the cell spacing. The actual geographic distance for each move is computed via the haversine formula rather than Euclidean approximation, ensuring correct distance calculation at all latitudes where the longitudinal spacing of a 0.5° grid varies from approximately 30 nm at the equator to approximately 15 nm at 60° latitude.

3.2 Grid Parameters

3.3 Haversine Distance

All distance calculations between grid cells use the haversine formula with the Earth radius set to 3440.065 nautical miles. This value is the standard mean Earth radius (6371 km) converted to nautical miles and provides sub-0.3% accuracy for the distance ranges encountered in ocean routing (tens to thousands of nautical miles).

The initial bearing from one cell to another is computed using the standard two-argument arctangent formula:

4. A* Search Algorithm

The A* search implemented in RouteOptimizer._astar_search() operates on the ocean grid constructed in the previous section. The algorithm maintains three principal data structures: (1) an open set implemented as a min-heap priority queue of SearchNode objects ordered by f-score, (2) a dictionary g_scores mapping grid cell coordinates (row, col) to the best-known accumulated cost from the origin, and (3) a closed set of fully expanded cells. Each SearchNode encapsulates the cell's f-score, the cell itself, the accumulated g-score, the estimated arrival time at that cell, and a reference to the parent node for path reconstruction.

4.1 Cost Function Components

The evaluation function for each node follows the standard A* formulation:

The g-score g(n) represents the actual accumulated fuel consumption from the origin to node n, computed by summing the fuel costs of all legs along the path. Each leg's fuel cost is evaluated by _calculate_move_cost(), which queries the weather provider at the midpoint of the leg for the estimated arrival time, computes fuel consumption via the vessel hydrodynamic model, adjusts for ocean current effects on speed over ground, and applies safety and zone penalty multipliers. The cost returned is the product of fuel consumption in metric tonnes, the safety cost factor, and any zone penalty factor.

4.2 Admissible Heuristic

The heuristic h(n) must never overestimate the true remaining cost for A* to guarantee optimality. WindMar computes the heuristic as the great-circle distance from the current cell to the destination, divided by the maximum attainable service speed, multiplied by the fuel consumption rate, and then scaled by a factor of 0.8. The 0.8 multiplier provides a deliberate underestimate: even under the most favourable conditions (following currents, calm seas), the actual fuel consumption per nautical mile will exceed this conservative estimate. This ensures admissibility across all possible weather conditions while still providing sufficient guidance to focus the search toward the destination.

4.3 Node Expansion

At each iteration, the algorithm pops the node with the lowest f-score from the open set. If this node is the destination cell, the search terminates and the optimal path is reconstructed by following parent pointers from the destination back to the origin. Otherwise, the node is added to the closed set and its eight neighbours are examined. For each neighbour that is present in the ocean grid (i.e., not land and not outside the corridor) and not already in the closed set, the algorithm computes the tentative g-score as the current node's g-score plus the move cost to the neighbour. If this tentative g-score is lower than any previously recorded g-score for that neighbour, the neighbour's record is updated and it is pushed onto the open set with the new f-score. The search terminates either when the destination is reached or when the number of explored cells exceeds max_cells (default 50,000), a safeguard against excessive computation in pathological cases.

4.4 Pseudocode

OPEN ← priority queue with start node (f=h(start), g=0)
CLOSED ← empty set
g_scores ← {start: 0}

while OPEN is not empty:
    node ← pop minimum f from OPEN
    if node is destination: return reconstruct_path(node)
    add node to CLOSED

    for each neighbor in 8-connected(node):
        if neighbor in CLOSED or neighbor is land: continue

        cost, time ← calculate_move_cost(node, neighbor)
        tentative_g ← node.g + cost

        if tentative_g < best_g[neighbor]:
            best_g[neighbor] ← tentative_g
            neighbor.parent ← node
            neighbor.arrival_time ← node.arrival_time + time
            push neighbor to OPEN with f = tentative_g + h(neighbor)

return failure (no path found within max_cells)

4.5 Edge Cost Evaluation

The _calculate_move_cost() method is the most computationally intensive component of each node expansion. For a candidate move from cell A to cell B, the method performs the following steps. First, it computes the great-circle distance and initial bearing between the two cells using the haversine formula. Second, it queries the weather provider at the geographic midpoint of the leg for the estimated arrival time at that midpoint, obtaining wind, wave, and current conditions. Third, it passes the environmental conditions and vessel heading to the vessel model's calculate_fuel_consumption() method, which returns the fuel consumed in metric tonnes for the leg. Fourth, it computes the effect of ocean currents on speed over ground: the current velocity component along the vessel's heading is projected using effect = current_kts × cos(relative_angle), and the speed over ground is computed as sog = stw + current_component. Fifth, the safety cost factor and zone penalty are retrieved and multiplied into the final cost. The method returns a tuple of (cost, travel_time_hours), where the cost is fuel_mt × safety_factor × zone_penalty.

5. Per-Leg Voyage Calculation

Once the A* search has produced an optimal waypoint sequence, the VoyageCalculator class performs a detailed per-leg analysis of the route. This class is also invoked during the search itself (via the move cost function) to evaluate candidate legs, but its primary role is to produce the comprehensive VoyageResult that includes per-leg breakdowns of distance, bearing, weather, speed loss, fuel consumption, and power requirements. The calculator processes each leg sequentially, accumulating departure and arrival times to ensure that weather queries reflect the vessel's actual temporal position along the route.

5.1 Leg Performance Calculation

For each leg, the _calculate_leg_performance() method determines the vessel's achievable speed through water (STW) and the corresponding fuel consumption. The method first evaluates the total resistance at the commanded calm-water speed using the hydrodynamic model, which sums calm-water hull resistance, wind resistance (Blendermann model), and wave-added resistance (STA-wave method). The required propulsive power is then compared against 90% of the vessel's maximum continuous rating (MCR). If the required power is within limits, the vessel maintains the commanded speed. If the required power exceeds 90% MCR, the vessel must reduce speed involuntarily, and the reduction factor is computed using the cube law relationship between power and speed.

After determining the actual STW (with or without involuntary speed reduction), the method recalculates resistance and fuel consumption at the reduced speed. The speed loss percentage is recorded in the LegResult for diagnostic purposes, enabling the user to identify legs where heavy weather forced a significant departure from the planned speed.

5.2 Current Effect on Speed Over Ground

The speed over ground (SOG) is computed by vector addition of the vessel's STW along its heading and the ocean current vector. The VoyageCalculator._calculate_sog() method decomposes both the vessel velocity and the current velocity into north and east components, sums the components, and computes the magnitude of the resultant vector. This approach correctly handles currents from any direction relative to the vessel heading, including beam currents that alter the vessel's course-made-good without changing the scalar SOG along the planned track.

In the simplified current model used within the A* move cost evaluation (_calculate_current_effect()), the current effect is projected along the vessel heading as a scalar component: the relative angle between the vessel's heading and the current direction is computed, and the along-track current component is current_kts × cos(relative_angle). A positive component (following current) increases SOG and reduces passage time; a negative component (opposing current) decreases SOG and increases both passage time and fuel consumption per nautical mile of ground track.

5.3 Per-Leg Output Fields

5.4 Weather Data at Leg Midpoint

Weather conditions for each leg are queried at the geographic midpoint of the leg and at the estimated time of arrival at that midpoint. The LegWeather dataclass encapsulates the full set of environmental parameters: wind speed and direction, significant wave height and period, wave direction, current speed and direction, and (when available) decomposed wind-wave and swell components. The decomposition flag has_decomposition determines whether the seakeeping model evaluates vessel motions using the total sea state or the more accurate decomposed approach described in the Hydrodynamics & RAO article.

6. Route Smoothing

The raw path produced by the A* search consists of waypoints aligned to the regular grid, resulting in a staircase-like trajectory with unnecessary intermediate points. The RouteOptimizer._smooth_path() method applies the Douglas-Peucker line simplification algorithm (Douglas and Peucker, 1973) to reduce the waypoint count while preserving the essential shape of the route. This is a standard cartographic technique that recursively identifies and removes waypoints whose perpendicular distance from the line connecting their neighbours falls below a specified tolerance.

WindMar configures the simplification tolerance at tolerance_nm=5.0 nautical miles, meaning that any waypoint whose removal would cause the route to deviate by less than 5 nm from its original shape is discarded. This value represents a balance between route fidelity (preserving the optimiser's chosen deviations around weather systems) and waypoint economy (reducing the number of legs that must be individually evaluated for fuel consumption and safety).

6.1 Douglas-Peucker Algorithm

The algorithm operates recursively on a sequence of waypoints. Given the endpoints of a segment, it finds the intermediate waypoint with the greatest perpendicular distance from the line connecting the endpoints. If this distance exceeds the tolerance, the waypoint is retained and the algorithm recurses on the two sub-segments. If no intermediate waypoint exceeds the tolerance, all intermediate waypoints are removed and the segment is replaced by a direct line. The recursion produces a minimal set of waypoints that approximates the original route within the specified tolerance.

6.2 Land-Crossing Check

A critical extension to the standard Douglas-Peucker algorithm in the maritime context is the land-crossing check. Before simplifying any segment by removing intermediate waypoints, the algorithm verifies that the resulting direct line between the retained endpoints does not cross land. If the simplified segment would cross a landmass, the intermediate waypoint with the greatest distance is retained regardless of whether it exceeds the tolerance, and the recursion proceeds on the sub-segments. This ensures that the smoothed route remains entirely over water, even when the original grid-aligned path detours around peninsulas, islands, or narrow straits.

6.3 Re-evaluation

After smoothing, the route statistics are recalculated using _calculate_route_stats(), which iterates over the smoothed waypoint sequence and computes per-leg distance, bearing, weather, fuel consumption, and safety status. When variable-speed optimisation is enabled, each leg is also evaluated at multiple speeds (7 speed steps from 6.0 to 18.0 knots, as defined by SPEED_STEPS and SPEED_RANGE_KTS) to identify the speed that minimises fuel consumption per nautical mile while respecting safety constraints. The smoothed route is then compared against a direct great-circle route to compute fuel and time savings percentages, which are reported in the OptimizedRoute result.

7. Temporal Weather Integration

A distinguishing feature of WindMar's route optimisation is its use of time-varying weather data throughout the search and route evaluation process. Rather than evaluating all legs against a single weather snapshot at departure time, the system queries weather conditions at the estimated time of arrival at each grid cell during the A* search and at each leg midpoint during the detailed voyage calculation. The TemporalGridWeatherProvider supplies weather data as a function of both geographic position and time, drawing on the multi-timestep forecast grids described in the Data Pipeline article.

During the A* search, the SearchNode dataclass carries an arrival_time field that tracks the cumulative transit time from the origin. When a node is expanded and its neighbours are evaluated, the move cost function uses the node's arrival time plus half the estimated leg transit time as the query time for the weather provider. This means that the weather conditions used to evaluate a leg in the middle of a five-day voyage reflect the forecast for that future time, not the conditions at departure. The arrival time at the neighbour node is then updated as the parent's arrival time plus the full leg transit time, propagating the temporal position forward through the graph.

7.1 Impact on Route Selection

Temporal weather integration can produce routes that differ significantly from those computed using a static weather snapshot. Consider a vessel departing into calm conditions with a storm system forecast to arrive in 48 hours along the direct path. A static-weather optimiser would see no reason to deviate from the great-circle route, as the departure-time conditions are benign everywhere. The temporal optimiser, by contrast, evaluates the later legs against the forecast storm conditions and may route the vessel on a longer but calmer southerly arc that avoids the worst of the weather. This capability is particularly important for multi-day ocean crossings where weather systems can develop, intensify, and propagate during the voyage.

7.2 Iterative Time Accumulation

The temporal integration introduces a coupling between the route and the weather: the weather encountered depends on when the vessel arrives at each point, which depends on the speed achieved under the weather encountered at earlier points. In the A* search, this coupling is handled implicitly through the arrival time field in each SearchNode. In the post-search voyage calculation, the coupling is handled by sequential iteration over legs, where the departure time of each leg equals the arrival time of the previous leg. This sequential accumulation ensures that the weather query times reflect the vessel's actual temporal trajectory along the route, including any speed reductions caused by adverse conditions on earlier legs. The trilinear interpolation scheme used by the weather provider (spatial bilinear interpolation combined with temporal linear interpolation between forecast time steps) is described in the Data Pipeline article.

8. Performance Characteristics

The A* search with an admissible heuristic has a worst-case time complexity of O(N log N), where N is the number of ocean grid cells in the search corridor. The log N factor arises from the heap operations (push and pop) on the priority queue. In practice, the heuristic substantially reduces the number of cells explored: for a Mediterranean route at 0.5° resolution, the search typically explores 2,000–8,000 cells out of a total grid of 15,000–30,000 ocean cells, achieving an effective pruning ratio of 70–85%. The number of cells explored is recorded in the OptimizedRoute.cells_explored field and the total optimisation time in optimization_time_ms.

8.1 Computational Cost Breakdown

8.2 Variable-Speed Optimisation

When variable-speed optimisation is enabled, the _find_optimal_speed() method evaluates 7 discrete speed steps uniformly distributed across the range 6.0–18.0 knots (SPEED_RANGE_KTS). For each candidate speed, the method computes the fuel consumption, the effective SOG accounting for currents, the transit time, and the safety cost factor. The scoring criterion depends on the optimisation objective: for fuel optimisation, the method selects the speed that minimises fuel per nautical mile; for time optimisation, it selects the speed that minimises transit time. Safety penalties are applied multiplicatively to the score, ensuring that speeds which produce dangerous vessel motions are heavily penalised even if they are fuel-efficient. The optimal speed, fuel consumption, and transit time are returned as a tuple for incorporation into the leg results.

8.3 Memory Usage

Memory consumption is proportional to the number of grid cells. Each GridCell occupies approximately 64 bytes (four floats plus Python object overhead), and each entry in the g-scores dictionary adds approximately 100 bytes (tuple key plus float value plus hash table overhead). For the maximum grid size of 50,000 cells, the total memory footprint of the search data structures is approximately 8 MB, well within the constraints of a web application backend. The SearchNode objects in the priority queue are garbage-collected as nodes are popped and processed, so the queue's peak memory usage is bounded by the number of cells on the frontier at any given time, typically 500–2,000 nodes.

8.4 Typical Performance Metrics

References

1. Hart, P.E., Nilsson, N.J. and Raphael, B. (1968). “A Formal Basis for the Heuristic Determination of Minimum Cost Paths.” IEEE Transactions on Systems Science and Cybernetics, vol. 4, no. 2, pp. 100–107.
2. Dijkstra, E.W. (1959). “A Note on Two Problems in Connexion with Graphs.” Numerische Mathematik, vol. 1, pp. 269–271.
3. Hagiwara, H. (1989). “Weather Routing of (Sail-Assisted) Motor Vessels.” PhD Thesis, Delft University of Technology.
4. Szlapczynska, J. (2015). “Multi-objective Weather Routing with Customised Criteria and Constraints.” Journal of Navigation, vol. 68, no. 2, pp. 338–354.
5. Douglas, D.H. and Peucker, T.K. (1973). “Algorithms for the Reduction of the Number of Points Required to Represent a Digitized Line or its Caricature.” Cartographica, vol. 10, no. 2, pp. 112–122.
6. IMO (2021). “Guidance on Treatment of Innovative Energy Efficiency Technologies for Calculation and Verification of the Attained EEDI and EEXI.” MEPC.1/Circ.896.
7. Vettor, R. and Guedes Soares, C. (2016). “Development of a Ship Weather Routing System.” Ocean Engineering, vol. 123, pp. 1–14.