Forecast

The unified forecast reader. A single Forecast class reads the v2 canonical netCDF schema (see Canonical NetCDF Forecast Schema v2) regardless of whether the file originated from GFS or ERA5; the source is resolved from the file’s institution global attribute and exposed as Forecast.source.

Forecast is a thin orchestrator: it owns the public API — getNewCoord() (scalar) and getNewCoords() (batch) — and the trajectory advection, while delegating the wind lookup to pluggable backends (EarthSHAB.utils.forecast_backends), the altitude/time interpolation to EarthSHAB.utils.wind_interp, and the advection step to EarthSHAB.utils.advection. The wind field is selected by three orthogonal config fields — wind_interpolation, advection, and backend — all documented on Wind Interpolation Methods. These utils modules are shared byte-for-byte with the HAB-COM project so interpolation/physics edits stay portable between the two codebases.

Forecast reader for the v2 canonical netcdf forecast schema.

Single class for both GFS- and ERA5-sourced netcdf files. This module replaces both ERA5.py and GFS.py.

Standard v2 file format:

dims:       valid_time, pressure_level, latitude, longitude
data vars:  u, v, z, t  (z is geopotential m**2 s**-2; t may be absent)
coords:     latitude descending, longitude in [-180, 180) ascending,
            pressure_level descending hPa,
            valid_time 'seconds since 1970-01-01' / proleptic_gregorian
global:     Conventions = 'CF-1.7'
            institution = 'NOAA/NCEP (GFS)' or 'ECMWF (ERA5)'

Every netcdf file is a bounded forecast subset (no full-world arrays with mask-based subsetting), although full world forecasts can still be downloaded by applying full world lat/lon bounds. The reader slices the entire stored array.

Hot-swappable wind-lookup design

Three orthogonal axes select the wind field, each a constructor kwarg with a config.forecast[...] fallback:

• backend – ‘xarray’ (slow/original) | ‘numpy’ | ‘numba’ (default)

• wind_interpolation – ‘bilinear’ | ‘linear_neighbors’ | ‘linear_full’ | ‘spline_full’

• advection – ‘geodesic’ (WGS84) | ‘tangent_plane’ (fixed-anchor spherical)

The skeleton (this class) owns the public API – getNewCoord (scalar) and getNewCoords (batch) – plus advection. The tier swaps only the wind lookup, implemented by the backend classes in utils.forecast_backends; the vertical altitude/time math is shared (utils.wind_interp). getNewCoord is a thin length-1 wrapper over getNewCoords so the two share one code path.

The shape mirrors EarthSHAB’s Forecast.py (same Forecast(start_coord) construction and getNewCoord return) so this is upstreamable.

class Forecast.Forecast(start_coord=None, *, forecast_file=None, start_time=None, dt=None, sim_time=None, wind_interpolation=None, advection=None, backend=None)

Reader for v2 canonical forecast files (GFS- or ERA5-sourced).

__init__(start_coord=None, *, forecast_file=None, start_time=None, dt=None, sim_time=None, wind_interpolation=None, advection=None, backend=None)

close()

Close the underlying netCDF dataset (and any backend dataset).

CPython’s __del__ on netCDF4.Dataset is unreliable enough that long-running batch evaluators accumulate HDF5 chunk-cache state across launches and segfault. Callers should close() explicitly. Idempotent.

closest(arr, k)

Given an ordered array and a value, return the index of the closest item.

closestIdx(arr, k)

Given an ordered array and a value, return the index of the closest item.

datetime_epoch(timedate_obj)

Convert a datetime (or ISO string) to GMT epoch seconds.

fill_missing_data(data)

Linearly interpolate over missing entries in a 1-D profile column.

getNearestAlt(hour_index, lat, lon, alt)

GFS-style signature: (hour_index, lat_deg, lon_deg, alt_m).

getNearestAltbyIndex(int_hr_idx, lat_i, lon_i, alt_m)

Nearest altitude index at a specific (time, lat, lon) cell.

getNewCoord(coord, dt)

Advance a single balloon coord by dt seconds (EarthSHAB-shaped).

Thin length-1 wrapper over getNewCoords so scalar and batch share one code path (and are bit-identical at M=1).

Parameters:

• coord – dict with keys lat, lon, alt, timestamp

• dt – integration interval in seconds

Returns:

[lat_new, lon_new, x_wind_vel, y_wind_vel, x_wind_vel_diag, y_wind_vel_diag, bearing, nearest_lat, nearest_lon, nearest_alt]

getNewCoords(lats, lons, alts, timestamp, dt, *, diagnostic=True)

Advance a batch of coords by dt seconds under the local wind field.

Parameters:

• lats/lons/alts – 1-D arrays (length M) of current positions

• timestamp – shared datetime for all members

• dt – integration interval in seconds

• diagnostic – also compute the linear_full diagnostic winds and the nearest-cell diagnostic fields (needed for the full scalar return / EarthSHAB contract). Set False in hot batch loops to skip the extra wind sample + per-member nearest-index lookups.

Returns:

list of 10 length-M arrays mirroring getNewCoord’s 10-tuple: lat_new, lon_new, x_wind, y_wind, x_wind_diag, y_wind_diag, bearing, nearest_lat, nearest_lon, nearest_alt.

wind_alt_Interpolate2(alt_m, diff_time, lat_idx, lon_idx, lat=None, lon=None)

Scalar wind interpolation -> [u, v, u_diag, v_diag].

Kept for API parity with EarthSHAB’s Forecast.py. Routes through the selected backend; (u_diag, v_diag) is always the linear_full path. For ‘bilinear’ the actual lat/lon are required (nearest-cell methods fall back to the grid coordinates of the supplied indices).