contrailopt.dag

Utilities for horizontal DAG construction.

Functions

validate_flight_profile(ds, n_nodes)

Validate and format a (waypoint, altitude_ft) flight profile for the track solver.

Classes

AirportCoords(*, icao_code, longitude, ...)	Coordinates and elevation of an airport, identified by ICAO code.
EdgeInterpolation(*, air_temperature, ...)	Met fields interpolated at sample points.
EdgeMetLookup(*, ds, edge_ptr, edge_idx, ...)	Pre-interpolated met data on edge sample points.
HorizontalDAG(*, lon, lat, adj_ptr, adj, ...)	Directed graph on (lon, lat) nodes with CSR adjacency.
Track(*, lon, lat, node_time)	An ordered sequence of timed waypoints along a flown path.

class contrailopt.dag.AirportCoords(*, icao_code: str, longitude: float, latitude: float, elevation_ft: float)

Bases: object

Coordinates and elevation of an airport, identified by ICAO code.

classmethod from_icao(icao_code: str) → Self

Look up airport coordinates by ICAO code.

property coords: tuple[float, float]

Return (longitude, latitude) coordinates as a tuple.

class contrailopt.dag.HorizontalDAG(*, lon: NDArray[floating], lat: NDArray[floating], adj_ptr: NDArray[int64], adj: NDArray[int64], edge_dist: NDArray[floating], h_origin: int, h_dest: int)

Bases: object

Directed graph on (lon, lat) nodes with CSR adjacency.

All geometry (distances, azimuths, interpolation, polygon exclusion) is computed on a sphere, not on a planar lon/lat grid.

lon: NDArray[floating]

Longitude of each node in degrees (n,). Assumed to be in the range [-180, 180). (This assumption is used in crosses_antimeridian()).

lat: NDArray[floating]

Latitude of each node in degrees (n,). Assumed to be in the range [-90, 90].

adj_ptr: NDArray[int64]

CSR row pointers (n + 1,). Neighbors of node i are adj[adj_ptr[i]: adj_ptr[i+1]].

adj: NDArray[int64]

Neighbor (destination) indices for each directed edge (m,).

edge_dist: NDArray[floating]

Great-circle distance in meters for each directed edge (m,).

h_origin: int

Index of the distinguished origin node.

h_dest: int

Index of the distinguished destination node.

property n_nodes: int

The number of nodes in the graph.

property n_edges: int

The number of directed edges in the graph.

property out_degree: NDArray[int64]

The out-degree of each node.

property edge_src: NDArray[int64]

The source endpoint index of each directed edge.

property edges: NDArray[int64]

The directed edges of the graph as (src, dest) index pairs in an (m, 2) array.

property crosses_antimeridian: bool

Determine if the great circle between origin and destination crosses the antimeridian.

neighbors(i: int) → NDArray[int64]

Return the neighbors of a specified node.

edge_index(src: int, dst: int) → int

Return the CSR index of the directed edge (src, dst).

Performs a linear scan over the neighbors of src. This could be replaced with np.searchsorted if needed, but that would require enforcing sorted neighbors during construction.

Raises:

ValueError – If no edge from src to dst exists.

neighbors_batch(nodes: NDArray[int64]) → NDArray[int64]

Return neighbors (duplicates included with multiplicity) for a batch of nodes.

expand_neighbors(nodes: NDArray[int64]) → tuple[NDArray[int64], NDArray[floating], NDArray[int64], NDArray[int64]]

Expand CSR adjacency for a batch of nodes into flat edge arrays.

Returns:

• flat_nbr (numpy.ndarray) – (e,) neighbor indices for all edges leaving nodes.

• flat_dist (numpy.ndarray) – (e,) edge distances in meters.

• src_idx (numpy.ndarray) – (e,) index into nodes for each flat entry, so nodes[src_idx[k]] is the source node of flat edge k.

• flat_edge_idx (numpy.ndarray) – (e,) index of each edge in the CSR arrays (adj, edge_dist).

• Here ``e = out_degree[nodes].sum()``, the total number of outgoing

• edges from all ``nodes`.`

adjacency_matrix() → NDArray[bool]

Return dense boolean adjacency matrix A where A[i, j] is True for i->j.

distance_matrix(missing: float = inf) → NDArray[floating]

Return dense distance matrix A where A[i, j] is the distance aong i->j.

Distances are set to infinity by default where edges are missing.

reverse() → Self

Return a new DAG with all edge directions flipped and origin/dest swapped.

prune_unreachable() → Self

Return a new DAG with only nodes reachable from origin that also reach dest.

prune_edges(degree: int) → Self

Keep the best degree outgoing and incoming edges per node.

Each edge is scored by the sum of its azimuth deviations: how far the edge direction deviates from the azimuth toward the destination (at the tail) plus how far the reverse deviates from the azimuth toward the origin (at the head). An edge is kept if it ranks among the best degree outgoing edges of its source or among the best degree incoming edges of its destination.

Parameters:

degree (int) – Number of outgoing and incoming edges to keep per node.

Returns:

A new DAG with at most degree outgoing and incoming edges per node.

Return type:

Self

exclude_polygons(polygons: list[list[tuple[float, float]]]) → Self

Return a new DAG with edges crossing any polygon removed.

Uses spherely for geodesic intersection tests on the sphere.

Parameters:

polygons (list[list[tuple[float, float]]]) – List of polygons, where each polygon is a list of (lon, lat) vertices.

Returns:

A new DAG with offending edges removed and then pruned.

Return type:

HorizontalDAG

sample_edges(spacing_m: float) → tuple[NDArray[floating], NDArray[floating], NDArray[int64], NDArray[int64]]

Sample points along every edge at most spacing_m meters apart.

Points are uniformly spaced along each edge, and both edge endpoints (source and destination nodes) are included as samples.

Returns:

• sample_lon (numpy.ndarray) – (s,) longitude of each sample point.

• sample_lat (numpy.ndarray) – (s,) latitude of each sample point.

• edge_idx (numpy.ndarray) – (s,) edge index for each sample point.

• edge_ptr (numpy.ndarray) – (m + 1,) CSR-style pointer so edge i’s samples are at sample_lon[edge_ptr[i]: edge_ptr[i+1]].

• Here ``s = edge_ptr[-1]``, the total number of sample points across all edges.

plot(ax: GeoAxes | None = None, linewidth: float = 2.0, show_edges: bool = True) → GeoAxes

Plot the DAG on a cartopy map.

classmethod from_network(lon: NDArray[floating], lat: NDArray[floating], tail: NDArray[int64], head: NDArray[int64], origin_idx: int = 0, dest_idx: int = -1, max_angle_deg: float = 40.0) → Self

Build a DAG from a static network graph using the dual azimuth constraint.

classmethod from_points(lon: NDArray[floating], lat: NDArray[floating], origin_idx: int = 0, dest_idx: int = -1, max_angle_deg: float = 40.0, max_dist_m: float = 500000.0) → Self

Build a DAG from lon/lat arrays using the dual azimuth constraint.

classmethod from_poisson(origin_lon: float, origin_lat: float, dest_lon: float, dest_lat: float, poisson_spacing_m: float = 80000.0, max_cross_track: float | None = None, max_angle_deg: float = 40.0, max_dist_m: float = 500000.0, dtype: type[floating] = <class 'numpy.float64'>, rng: Generator | None = None) → Self

Build a DAG from Poisson-disk sampled points along the OD great circle.

topo_wavefronts() → Generator[NDArray[int64], None, None]

Yield topological wavefronts reachable from origin.

Only nodes reachable from h_origin are emitted. Unreachable nodes (those with incoming edges from outside the reachable subgraph) are excluded.

The first wavefront contains only the origin. Wavefront k contains nodes whose reachable in-degree drops to zero after removing wavefronts 0 … k-1.

class contrailopt.dag.Track(*, lon: NDArray[floating], lat: NDArray[floating], node_time: NDArray[datetime64])

Bases: object

An ordered sequence of timed waypoints along a flown path.

The flight-profile optimizer (contrailopt.optimize.solve_track()) follows a fixed lateral path, so it needs only the ordered nodes, their times, and along-track distances.

This interface supplies lon, lat, node_time, cum_dist, segment_dist, n_nodes, h_origin, h_dest computed straight from the coordinates, so it can often be used in place of HorizontalDAG.

property n_nodes: int

The number of waypoints.

property h_origin: int

Index of the origin node (always the first waypoint).

property h_dest: int

Index of the destination node (always the last waypoint).

property segment_dist: NDArray[floating]

Great-circle distance between consecutive waypoints (n - 1,).

property cum_dist: NDArray[floating]

Cumulative along-track distance at each waypoint (n,), starting at zero.

property crosses_antimeridian: bool

Whether the origin-to-destination span wraps the antimeridian.

plot(ax: GeoAxes | None = None, linewidth: float = 2.0) → GeoAxes

Plot the track on a cartopy map.

contrailopt.dag.validate_flight_profile(ds: Dataset, n_nodes: int) → Dataset

Validate and format a (waypoint, altitude_ft) flight profile for the track solver.

Returns a dataset carrying the met variables under their original names (air_temperature, u_wind, v_wind, and optional eef_per_m), each cast to float32 and oriented (waypoint, altitude_ft).

NaN in the core weather variables raises, while NaN in eef_per_m is zero-filled.

class contrailopt.dag.EdgeInterpolation(*, air_temperature: NDArray[floating], eastward_wind: NDArray[floating], northward_wind: NDArray[floating], eef_per_m: NDArray[floating] | None)

Bases: object

Met fields interpolated at sample points.

class contrailopt.dag.EdgeMetLookup(*, ds: Dataset, edge_ptr: NDArray[int64], edge_idx: NDArray[int64], sample_lon: NDArray[floating], sample_lat: NDArray[floating], cum_dist: NDArray[floating], delta_dist: NDArray[floating], sample_azimuth: NDArray[floating])

Bases: object

Pre-interpolated met data on edge sample points.

ds: Dataset

xr.Dataset with dims (sample, altitude_ft, time) containing weather variables interpolated onto edge sample coordinates.

edge_ptr: NDArray[int64]

CSR-style pointer array (n_edges + 1,). Samples for edge i are at indices edge_ptr[i]:edge_ptr[i+1].

edge_idx: NDArray[int64]

Edge index for each sample point (n_samples,).

sample_lon: NDArray[floating]

Longitude of each sample point (n_samples,).

sample_lat: NDArray[floating]

Latitude of each sample point (n_samples,).

cum_dist: NDArray[floating]

Cumulative distance from edge source to each sample point in meters (n_samples,).

delta_dist: NDArray[floating]

Distance in meters from this sample to the next (n_samples,). The last sample of each edge has delta_dist = 0. Equal to diff(cum_dist) within each edge.

sample_azimuth: NDArray[floating]

Azimuth in radians from each sample to the next (n_samples,). The last sample of each edge copies the previous sample’s azimuth.

classmethod from_met(met: MetDataset | Dataset, dag: HorizontalDAG, altitude_ft: NDArray[floating], takeoff_time: Timestamp, flight_hours: int, spacing_m: float, eef: DataArray | MetDataArray | None = None) → Self

Interpolate met data onto dag edge sample points.

Parameters:

• met (MetDataset | xr.Dataset) – Gridded met dataset with “air_temperature”, “eastward_wind”, and “northward_wind”. If “eef_per_m” is present, it will also be included in the output with NaN values filled to 0.0 (no EEF forecast is treated as zero forcing). NaN values in weather variables are not allowed and will raise an error. Either a pycontrails MetDataset or a raw xr.Dataset with similar structure can be passed.

• dag (HorizontalDAG) – Horizontal DAG whose edges will be sampled.

• altitude_ft (numpy.ndarray) – An array of altitudes in feet to interpolate onto.

• takeoff_time (pandas.Timestamp) – Departure time for the flight, used to select met time steps.

• flight_hours (int) – Number of hourly time steps to retain starting from takeoff_time.

• spacing_m (float) – Spacing in meters between sample points along edges. Passed to dag.sample_edges.

• eef (xr.DataArray | MetDataArray | None, default None) – Optional “eef_per_m” DataArray on its own lon/lat grid. If provided, EEF is interpolated onto sample points independently from the weather grid, avoiding the need to pre-merge onto a common grid. Takes precedence over “eef_per_m” in met if both are present. Assumed to adhere to pycontrails MetDataArray conventions.

Returns:

EdgeMetLookup with weather interpolated onto (sample, altitude_ft, time) dims.

Return type:

EdgeMetLookup