diff --git a/docs/internal/pyneats_units_audit.md b/docs/internal/pyneats_units_audit.md new file mode 100644 index 0000000..9896016 --- /dev/null +++ b/docs/internal/pyneats_units_audit.md @@ -0,0 +1,779 @@ +# PyNEATS unit audit against the AO upload schema + +Reviewed revision: `f7502fe6344e5e085df1042ac817f02a05bd100d` (`origin/main`, 2026-10-13), +on branch `docs/units-audit`. Every `file:line` reference below points to that revision, +before the comment patch was applied. The patch moves lines in the files it touches. + +## 1. Outcome + +- **Two known interface differences.** `tas` (knots in the AO JSON, m/s in PyNEATS) and + `calorific_value` (MJ/kg in the AO JSON, J/kg in PyNEATS) are handled by the NEATS + application before PyNEATS is called. The task owner confirmed this. The application source + was not available, so the conversions were not seen in code. PyNEATS consumes both values + in the units it expects: m/s and J/kg. No PyNEATS code converts either value again, and none + contradicts the confirmed conventions. +- **Most AO fields map consistently.** Latitude and longitude, flight level, take-off mass, + per-point aircraft mass and hydrogen content carry the same unit and convention at the + PyNEATS boundary as in the AO schema. Flight level becomes metres of ISA pressure altitude. + Hydrogen content is a weight percentage. +- **Additional discrepancies were found in the inspected scope.** The ones on the AO input + path are listed in Section 7: + - Non-UTC offsets in trajectory timestamps are silently discarded rather than applied. + - Some schema-valid timestamp forms, and a bare single-flight object, abort the run. + - Fuel flow derived from `am` is shifted by one sample. + - Partially null optional columns are not filled. + + Inside PyNEATS, away from the AO fields, the audit found unit errors with large numerical + effects (Section 7): + - The default local aCCF H2O formula receives potential vorticity in SI units, but its + coefficients are defined for PVU. The H2O term loses all sensitivity to the weather, and + H2O ATR is underestimated 1.6–3.3 times per flight. + - The small-emitter output `contrails_ef_J` contains an AGWP in W m-2 yr, not joules. + - OpenAirClim is told the NOx species is NO, while the NOx PyNEATS supplies is + NO2-equivalent. Method D O3 and CH4 are then about 53 % too high. + - The shipped OpenAirClim configuration pins the metric start year to 2025. +- **Most conventions the AO schema leaves undefined are defined in EU law** (Section 6.0). + Implementing Regulation (EU) 2024/2493 added Annex IIIa to the MRR (2018/2066). Its + definitions settle the following, and PyNEATS already follows each of them: + - `ff` is whole-aircraft fuel flow, kg/s (§1(33)); + - `ee` is the overall propulsion efficiency (§1(34) and Commission guidance v2.1); + - `load_factor` is payload mass ÷ maximum payload mass (§1(32)); + - `hydrogen_per_carbon_ratio` is atomic (§1(36)); + - `takeoff_mass` is the mass at the beginning of the take-off run (§1(29)); + - `calorific_value` is the net (lower) calorific value (§1(38)); + - `hydrogen_content` is a mass percentage (Commission guidance, Table 3). +- **Still open**, and no public source settles them (Section 6.0): + - the numeric scale of `aromatic_content`, `sulphur` and `naphthalene`. Their basis is set + by law (aromatics and naphthalene by volume, sulphur by mass), but not whether to send + 25 or 0.25; + - `ee` is a "percentage" in the law but a 0–1 fraction in the schema; + - EOBT versus AOBT: ICAO Item 13 is the EOBT, while PyNEATS labels the value AOBT; + - who chooses the NWP run for `nwp_date_time`; + - the MRR defines `tas` in m/s, while the AO schema uses knots. The upstream conversion + is right for uploads that follow the schema. +- **Computational behaviour is unchanged.** The patch adds and corrects comments and + docstrings only. An AST comparison confirms that no code changed (Section 11). + +## 2. Scope, material and method + +**Reference schema.** `ao_data_json_schema_supporting_arrays (4).json`, supplied as +`/data/rdcodina/pyneats/ao_data_json_schema_supporting_arrays (4) (1).json`. At top level it +accepts either one flight object or an array of flight objects. Paths below are relative to +`flight_information`. + +**Unavailable material.** +- The companion schema `ao_data_json_schema.json` was not supplied and is not in the + repository. It is not published either: the Commission guidance (v2.1, §3.3) places the + data-reporting templates on the NEATS User Group SharePoint, which requires a login. The two + schemas have **not** been compared. +- `scripts/schemas/flight_schema.json` is PyNEATS's own input schema. It is not the companion + schema. +- The NEATS application source was not available on this machine or in the + `eurocontrol-asu` GitHub organisation. Every statement about application preprocessing + therefore rests on the task owner's confirmation and is not verified in code. +- The `climaccf` package is not a declared dependency and is not installed. CLIMaCCF + behaviour was checked against a read-only clone of `github.com/dlr-pa/climaccf` (HEAD + `0184719`, 2025-06-20) and against the Dietmüller et al. 2023 and Matthes et al. 2023 GMD + papers. +- No BADA data files were read. pyBADA was inspected as source code only. + +**Dependency versions (from `uv.lock`).** pycontrails 0.60.2, pyBADA 0.1.12, +openairclim 0.13.0, xarray 2024.10.0. Findings that rest on a library interface quote paths +relative to `site-packages`. + +**Method.** +- Five parallel tracers followed every AO field from the JSON boundary through parsing, + interpolation, performance, emissions, climate functions and metrics to the outputs: + - the input boundary, time and position; + - fuel properties; + - performance; + - emissions, climate and outputs; + - runners and weather. +- They established units from arithmetic, constants, downstream library interfaces, tests + and golden data. A variable name or a plausible magnitude was not accepted as evidence. +- Tracer findings were merged. +- Two independent skeptics checked every finding. One asked whether the behaviour is real. + The other asked whether something else already accounts for it. +- A completeness critic then commissioned six gap investigations, and their findings were + verified in turn. +- The highest-impact claims were checked again by hand: PV units, `contrails_ef_J`, the + mass-derived fuel flow, the OpenAirClim `t_0`, BADA4 `VStall` and BADA3 `MPL`. +- No pipeline or BADA calculation was run. Numerical consequences quoted below come from + arithmetic on the committed golden outputs (`tests/data/golden/*_output.json`). + +**Evidence tags.** *Demonstrated* means shown by computation, a library interface, tests or +golden data. *Stated* means written in a comment, docstring or document only. + +**I/O variants covered.** +- JSON file input (`FleetRunnerParams.trajectory_json_filepath`, `scripts/run_json.py`). +- DataFrame input (`FleetRunnerParams.trajectory_dataframe`, `FlightRunner(source=df)`). +- Large emitter (Method C: CoCiP plus local aCCF or CLIMaCCF) and small emitter (Method D: + OpenAirClim). +- Performance-only runner. +- Flights with and without a trajectory, and every combination of provided `ff`, `ee` and + `am`. + +## 3. Audit table: AO input fields + +| JSON field | JSON unit or convention | PyNEATS boundary unit | Internal or output unit, if different | Existing conversion and layer | Status | Evidence | +| --- | --- | --- | --- | --- | --- | --- | +| *(top level)* object or array | `anyOf`: one `NeatsData` object or an array of them | list of flight objects only | — | none. A bare object is iterated by its keys and raises `AttributeError` | **New discrepancy** (D1) | `steps/parsing/neats_io.py:45-46`; `runners/fleet.py:561-566`; `scripts/run_json.py:238-240` *(demonstrated, snippet)* | +| `departure_date_time` | UTC, ISO 8601; "as indicated in Item 13 of the flight plan" | attr `aobt`: string copied verbatim | — | none. Never parsed or normalised; used only in flight keys and outputs | **Ambiguous / unverified** (A7): Item 13 is the EOBT (ICAO Doc 4444); the MRR does not choose; the PyNEATS label "AOBT" is unsupported | `neats_io.py:106`; `views.py:48-55`; golden output `aobt` equals the input string *(demonstrated)*; `flight_schema.json:38` and `output-format.md:50` say "AOBT" *(stated)* | +| `arrival_date_time` | UTC, ISO 8601 | mapped, then dropped | — | removed by the `allowed_attrs` filter, because it is not in the Flight4D attr lists | **Not consumed** | `neats_io.py:98,107,122`; `steps/parsing/views.py:35-55` *(demonstrated)* | +| `nwp_date_time` | UTC, ISO 8601, nullable; the NWP instance to use | not read | — | none. The NWP instance is whatever weather store or factory the caller supplies | **Not consumed** (open question, A8) | `grep -rn nwp` finds no code; `runners/fleet.py:578-584` *(demonstrated)* | +| `aircraft_properties.takeoff_mass` | kg, 0.1 kg resolution | attr `takeoff_weight`, kg | kg; pyBADA mass | rename only | **Consistent**: MRR Annex IIIa §1(29) defines it as kg "at beginning of the take-off run" (usage note in A6) | `neats_io.py:112`; `steps/performance/bada_model.py:786-808,849` *(demonstrated)* | +| `aircraft_properties.load_factor` | "unitless", [0, 1]; basis not stated | attr `payload_factor`, fraction | ZFW = OEW + pf × MPL in kg, or pf × (MTOW − OEW) when MPL is unavailable, which is always the case for BADA3 | rename only. Default 1.0. Not clipped | **Consistent**: defined by MRR Annex IIIa §1(32) as a fraction of maximum payload mass; default 1 per §5(2)(b) (A3) | `neats_io.py:113`; `bada_model.py:963-980`; `neats_default_parameters.py:141-143` *(demonstrated)*; `flight_schema.json:78-81` "payload / max-payload ratio" *(stated)* | +| `fuel_properties.hydrogen_content` | "percentage of hydrogen mass content", 0.01 resolution | wt.% (for example 13.79) | wt.%. Scales `ei_h2o` = 1.23 × (wt.% / 13.8) kg/kg | none | **Consistent** | `neats_io.py:115`; `core/neats_fuel.py:81-82,106`; pycontrails `core/fuel.py:53` (JetA 13.8); golden inputs 13.79 *(demonstrated)*. PyNEATS docs that say "fraction" are in Section 9 | +| `fuel_properties.hydrogen_per_carbon_ratio` | "unitless"; atomic or mass basis not stated | attr `h_c_ratio`, treated as an **atomic (molar)** H/C ratio | converted to wt.% H = 100 · r·1.008 / (12.011 + r·1.008); used only when `hydrogen_content` is absent | conversion in `NEATSFuel.__init__` | **Consistent**: MRR Annex IIIa §1(36) defines it as the number of H atoms per C atom (A4) | `neats_fuel.py:83-87`; range check: atomic [1.9061, 2.1857] maps exactly to [13.79, 15.50] wt.% (`tests/test_cases_pipeline.py:1492-1516`, mapping recomputed by the auditor) | +| `fuel_properties.aromatic_content` | "unitless"; scale and basis not stated | attr `aromatic_content`, stored | not used in any calculation | none | **Not consumed**. Volume basis per MRR §5(4) default; JSON scale undefined (A5) | `neats_io.py:117`; `neats_fuel.py:94-96`; default 0.25 commented as 25 % (`neats_default_parameters.py:300-302`) *(stated)* | +| `fuel_properties.calorific_value` | MJ/kg of fuel | attr `q_fuel`, **J/kg** (lower heating value) | J/kg in pycontrails `Fuel.q_fuel` and in the BADA fuel-flow correction ratio | × 1e6 in the **NEATS application** (owner-confirmed, not inspected). PyNEATS performs no conversion | **Known difference handled upstream** | `neats_io.py:118`; `neats_fuel.py:92`; `neats_default_parameters.py:315-320` (43.13e6, 42.8e6); `bada_model.py:765-771,841,875`; golden inputs 4.1944e7 *(demonstrated)* | +| `fuel_properties.sulphur` | "unitless"; scale and basis not stated | attr `sulphur_content`, stored | not used. `ei_so2` is always the JetA value of 1.2e-3 kg/kg | none | **Not consumed**. Mass basis per MRR §5(4) default; JSON scale undefined (A5) | `neats_io.py:119`; `neats_fuel.py:97-98,107-108`; pycontrails `emissions.py:770-771` *(demonstrated)* | +| `fuel_properties.naphthalene` | "unitless"; scale and basis not stated | dropped | — | key `naphthalene` in `neats_io` against `naphtalene` in `Flight4D.ATTRS_OPTIONAL`, so the value never reaches the fuel object (I2) | **Not consumed** | `neats_io.py:120,122`; `views.py:45`; `neats_fuel.py:195`; snippet *(demonstrated)* | +| `trajectory.trajectory_data[].ts` | UTC, ISO 8601 date-time | column `time`: naive string `%Y-%m-%d %H:%M:%S`, read as UTC | tz-naive UTC `datetime64[ns]`; durations in s (forward differences); resampled to 60 s by default | `strptime(...%z)` then `strftime` **discards the offset without applying it** | **New discrepancy** (D2), conditional on non-UTC offsets; UTC input is correct | `neats_io.py:61-68`; `neats_parser.py:233-242`; pycontrails `core/vector.py:292-304,2162-2167`; snippet *(demonstrated)* | +| `trajectory.trajectory_data[].lat` | decimal degrees, WGS 84 | column `latitude`, degrees | degrees | none | **Consistent** | `neats_io.py:73`; golden output equals input *(demonstrated)* | +| `trajectory.trajectory_data[].lon` | decimal degrees, WGS 84, [-180, 180] | column `longitude`, degrees | degrees; no wrapping | none | **Consistent** | `neats_io.py:74` *(demonstrated)* | +| `trajectory.trajectory_data[].fl` | flight level in hundreds of feet | column `altitude` in FL until the parser | m of **ISA pressure altitude** (FL × 30.48). pycontrails derives `level` (hPa) and `air_pressure` (Pa) from it | × 100 then `ft_to_m` (× 0.3048) in `NeatsTrajectoryParser` | **Consistent** | `neats_io.py:77`; `neats_parser.py:217-225`; pycontrails `core/vector.py:1503-1505,1544-1545`; golden `altitude[0]` = 121.92 m for FL 4 *(demonstrated)* | +| `trajectory.trajectory_data[].tas` | knots | column `true_airspeed`, **m/s** | m/s (BADA adapters, `overall_propulsion_efficiency`). Derived from ground speed and wind (m/s) when absent | × 1852/3600 in the **NEATS application** (owner-confirmed, not inspected). PyNEATS performs no conversion | **Known difference handled upstream**. The MRR (§1(4)) defines TAS in m/s; see §6.0 | `neats_io.py:82`; `bada_model.py:623-637,868,1050`; pycontrails `physics/jet.py:155-178` *(demonstrated)* | +| `trajectory.trajectory_data[].ff` | kg/s; whole-aircraft or per-engine not stated | column `fuel_flow`, kg/s, consumed as **whole aircraft** | `fuel_burn` = ff × segment duration, in kg per forward segment | none. Never scaled by engine count | **Consistent**: MRR Annex IIIa §1(33) defines whole-aircraft flow; some RSTS tables say "each engine" (A1) | `neats_io.py:79`; `bada_model.py:693,1071-1076`; pycontrails `emissions.py:183-195` (`fuel_flow / n_engine`), `cocip.py:929-931` (`fuel_flow / true_airspeed`) *(demonstrated)* | +| `trajectory.trajectory_data[].ee` | "unitless", [0, 1]; definition not stated | column `engine_efficiency`, fraction | consumed as **overall propulsion efficiency** F·TAS / (ff·q_fuel) in the SAC slope, proportional to 1/(1 − η). Computed values are clipped to [0, 0.5]; provided values are not clipped | none | **Consistent** in definition: MRR §1(34) and Commission guidance give overall propulsion efficiency; scale is "percentage" in law but [0, 1] in the schema (A2) | `neats_io.py:80`; `bada_model.py:1046-1063`; pycontrails `physics/jet.py:178,189`, `models/sac.py:206` *(demonstrated)* | +| `trajectory.trajectory_data[].am` | kg, 0.1 kg resolution | column `aircraft_mass`, kg | kg. When `ff` is absent, `ff` = −Δam/Δt in kg/s (D3) | none | **Consistent** (derived `ff` in D3) | `neats_io.py:81`; `bada_model.py:705-749,859-860` *(demonstrated)* | +| `flight_identification`, `departure_airport`, `arrival_airport` | identifiers | attrs `flight_id`, `departure_airport`, `arrival_airport` | — | none (upper-cased for the airport-fuel lookup) | **Not applicable** (no unit) | `neats_io.py:102-104` | +| `aircraft_properties.aircraft_type`, `aircraft_version`, `engine_uid` | identifiers; `engine_uid` may be a list of up to 4 | attrs `aircraft_type`, `aircraft_series` (renamed), `engine_uid` | FleetRunner splits a list into one copy per engine id and averages numeric results; FlightRunner rejects a list | rename only | **Not applicable** (no unit) | `neats_io.py:109-111`; `runners/fleet.py:612-640,385-390`; `neats_parser.py:298-303` | +| `confidential`, `ets_exemption` | flag, letter | not read | — | — | **Not consumed** (no unit) | `neats_io.py:99-123` | +| `trajectory.trj_data_source` | `Actual` or `Planned` | not read. `model_type` is the function argument, default `"NM"` | — | — | **Not consumed** (no unit) | `neats_io.py:22,105`; `runners/fleet.py:566` | + +### Relevant mappings beyond the AO JSON path + +| Mapping | Convention | PyNEATS boundary unit | Internal or output unit | Existing conversion and layer | Status | Evidence | +| --- | --- | --- | --- | --- | --- | --- | +| DataFrame entry (`trajectory_dataframe`, `FlightRunner(source=df)`): `altitude`, `time` | Not defined by AO upload schema | `altitude` in **FL**, not metres. `time` as `datetime64` (naive means UTC; tz-aware is converted) or a string in `date_format` (default `%Y-%m-%d %H:%M:%S`) | as for the JSON path | the same parser conversions | **Ambiguous / unverified** (A9) | `neats_parser.py:141,217-258`; example comment `examples/fleet_computation_from_dataframe.py:172` *(stated)*; snippet: ISO strings become NaT and are dropped *(demonstrated)* | +| Airport fuel-property file (`airport_fuel_path`) | Not defined by AO upload schema | `float()` only. `q_fuel` assumed J/kg; `h_c_ratio` used as atomic | J/kg, wt.% after `NEATSFuel` | **none**. PyNEATS reads the file directly, so the application's × 1e6 does not apply to it | **Ambiguous / unverified** (A10) | `neats_parser.py:35-41,44-120,284-296`; `tests/test_airport_qfuel.py` (J/kg) | +| Default branch: no `calorific_value` | Not defined by AO upload schema | `DEFAULT_Q_FUEL` = 42.8e6 J/kg | the BADA fuel flow is still multiplied by 43.13/42.8 = 1.0077 | correction ratio in `BADAPerformanceModel` | **Ambiguous / unverified** (A11) | `neats_fuel.py:92`; `bada_model.py:565-571,764-772,841,875`; `docs/explanation/mrv-specification.md:220` *(stated)* | +| `REFERENCE_Q_FUEL` against BADA4 LHV | Not defined by AO upload schema | 43.13e6 J/kg (pycontrails JetA) | — | not read from BADA files, although BADA4 carries a per-aircraft LHV | **Ambiguous / unverified** (A12) | `neats_default_parameters.py:315-317`; pyBADA `bada4.py:1855-1874`; pyBADA dummy BADA4 files: 43.217e6 and 43.71e6 J/kg | +| Per-point nulls in `tas`, `ff`, `ee`, `am` | schema allows `null` per point | a column survives unless every value is null | NaN spreads through interpolation and is skipped by sums | none | **New discrepancy** (D4) | `neats_io.py:93-94`; `steps/interpolation/pycontrails_interpolation.py:108-112`; `bada_model.py:1071-1076`; snippet *(demonstrated)* | + +## 4. Audit table: public outputs + +The AO upload schema defines inputs only. It is not an output contract, so every row below +reads "Not defined by AO upload schema" and records the convention the code produces. +`runner.results` from `FleetRunner` carries only flight-level entries. Per-waypoint columns +appear in `FlightView.to_dict()` dumps (the golden files) and in in-memory `Flight` objects. + +| Output field | JSON unit or convention | PyNEATS unit | Basis | Status | Evidence | +| --- | --- | --- | --- | --- | --- | +| `flight_information.fuel_burn_kg` | Not defined by AO upload schema | kg | flight total: nansum of per-segment `fuel_burn` | Consistent internally. Undercounts with partial nulls (D4) | `climate_metrics/gwp.py:359-379,434`; golden 3773.20 kg | +| `flight_information.co2_baseline_kg` | Not defined by AO upload schema | kg CO2 | flight total (`attrs['total_co2']`) | Consistent | `gwp.py:353,433`; golden 11919.55 = 3.159 × 3773.20 | +| `flight_information.contrails_ef_J` | Not defined by AO upload schema | **J** for Method C (sum of CoCiP `ef`). **W m-2 yr** for Method D (`AGWP_20_CONT`, no efficacy) | flight total | **New discrepancy** for Method D (O2) | `gwp.py:411,426`; openairclim `write_output.py:48-51` | +| `climate_metrics[species][H].EAGWP_Wm2yr` | Not defined by AO upload schema | W m-2 yr | flight total, efficacy-weighted; H = 20, 50, 100 yr | Consistent with code docstrings. The MRV spec label differs (Section 9) | `gwp.py:108-172,385-399`; golden recomputation ratio 1.0 | +| `climate_metrics[species][H].CO2eq_kg` | Not defined by AO upload schema | kg CO2-equivalent | flight total = EAGWP / C(H); for CO2 it equals `total_co2` | Consistent | `gwp.py:146-151,384`; `core/physics.py:62-67` | +| Low-altitude zero records | Not defined by AO upload schema | `fuel_burn_kg` and `co2_baseline_kg` are NaN; every species, CO2 included, is 0.0 | flight | Ambiguous (O9) | `runners/fleet.py:295-319`; `tests/test_zero_climate_result.py` | +| `fuel_flow` (column) | Not defined by AO upload schema | kg/s, whole aircraft | instantaneous rate per waypoint (AO value, −Δam/Δt, or BADA × q_ref/q) | Consistent (D3 for the derived branch) | `bada_model.py:693,741,875,1071-1074`; pycontrails `emissions.py:183-195` | +| `fuel_burn`, `fuel` (columns) | Not defined by AO upload schema | kg | per forward segment [tᵢ, tᵢ₊₁]; last value NaN | Consistent | `bada_model.py:693,1076-1079`; pycontrails `core/flight.py:1966-1969` | +| `aircraft_mass`, `thrust` (columns) | Not defined by AO upload schema | kg; N (whole aircraft) | instantaneous | Consistent | `bada_model.py:1066-1067`; `bada_adapters.py:279` | +| `true_airspeed`, `ground_speed`, `segment_duration`, `rocd`, `acceleration`, `delta_tau` (columns) | Not defined by AO upload schema | m/s, m/s, s, **ft/min**, m/s², K | per segment (forward) or per point | Consistent. `rocd` is the only non-SI kinematic column | `bada_model.py:622-674`; pycontrails `core/flight.py:617-717` | +| `engine_efficiency` (column) | Not defined by AO upload schema | dimensionless fraction | instantaneous; computed values clipped to [0, 0.5] | Consistent (definition question A2) | `bada_model.py:1056-1063` | +| `co2`, `h2o`, `so2`, `sulphates`, `oc`, `nox`, `co`, `hc`, `nvpm_mass` (columns) | Not defined by AO upload schema | kg | per segment = `fuel_burn` × EI | Consistent | pycontrails `emissions.py:760-777` | +| `nvpm_number` (column) | Not defined by AO upload schema | count (particles) | per segment | Consistent | pycontrails `emissions.py:777` | +| `nox_ei`, `co_ei`, `hc_ei`, `nvpm_ei_m`; `nvpm_ei_n` (columns) | Not defined by AO upload schema | kg/kg fuel; particles/kg fuel | emission index, not mass. NOx is NO2-equivalent by ICAO EDB convention | Consistent | pycontrails `emissions.py:263-295,324,476-479`; `models/apcemm/inputs.py:70` | +| `total_co2`, `total_h2o`, `total_so2`, `total_sulphates`, `total_oc` (attrs) | Not defined by AO upload schema | kg | flight total (nansum) | Consistent | pycontrails `emissions.py:780-788` | +| `total_nox`, `total_co`, `total_hc`, `total_nvpm_mass`; `total_nvpm_number` (attrs) | Not defined by AO upload schema | kg; count | flight sum **over waypoints inside the met domain only** (NaN EI elsewhere) | Ambiguous (O13) | golden EAF2143: 29 of 93 segments, 31.8 % of fuel, excluded | +| `ATR_20_O3`, `ATR_20_CH4`, `ATR_20_H2O` (columns) | Not defined by AO upload schema | K | per segment: pulse ATR over 20 yr, no efficacy; CH4 includes PMO (× 1.29) | `ATR_20_H2O`: **New discrepancy** (O1) | `climate_functions/local_accf.py:152-253,300-306` | +| `ef` (column, CoCiP) | Not defined by AO upload schema | J | contrail energy forcing per waypoint, summed over the contrail lifetime (energy, not forcing) | Consistent | pycontrails `contrail_properties.py:1485-1488`, `cocip.py:1434,1448` | +| `rf_sw_mean`, `rf_lw_mean`, `rf_net_mean`, `sdr_mean`, `rsr_mean`, `olr_mean` (columns) | Not defined by AO upload schema | W m-2 | local contrail RF and fluxes averaged over the contrail lifetime (forcing, not energy) | Consistent | pycontrails `cocip.py:1436-1446,2022-2093` | +| `AGWP_{20,50,100}_{CH4,O3,H2O,CONT}` (attrs, Method D) | Not defined by AO upload schema | W m-2 yr | flight total, one-year pulse, no efficacy; CH4 includes PMO | Consistent units. Window anchored to 2025 (O4) | `open_airclim.py:251-268`; openairclim `write_output.py:48-51` | +| `time` (`FlightView.to_dict`) | Not defined by AO upload schema | integer Unix seconds, truncated, from tz-naive UTC | instantaneous | Consistent; was undocumented (O11) | `core/views.py:229-232`; golden 1752085020 = first input `ts` | +| `contrail_age` (`FlightView.to_dict`) | Not defined by AO upload schema | integer **nanoseconds** | per waypoint | Consistent; was undocumented (O11) | pycontrails `utils/json.py:60` (NumpyEncoder) | +| Performance-only runner: `timestamps`, `true_airspeed`, `fuel_flow`, `aircraft_mass`, `engine_efficiency` | Not defined by AO upload schema | ISO 8601 **without offset** (naive UTC), m/s, kg/s, kg, fraction | per waypoint; one record per engine copy | Ambiguous (O12) | `runners/performance_runner.py:32,80-84,146-158` | + +## 5. Known differences handled upstream + +| Quantity | AO JSON unit | PyNEATS input unit | Conversion | Layer | Provenance | PyNEATS consumers checked | +| --- | --- | --- | --- | --- | --- | --- | +| `tas` → `true_airspeed` | knots | m/s | × 1852/3600 | NEATS application, before PyNEATS is called | **confirmed by the task owner**; application source not inspected | BADA adapters (`bada_model.py:868`), `overall_propulsion_efficiency` (`bada_model.py:1050`; pycontrails `jet.py:155-178`), the pycontrails Emissions/CoCiP interfaces. All treat the value as m/s. No knots handling exists in `src/` | +| `calorific_value` → `q_fuel` | MJ/kg of fuel | J/kg of fuel (LHV) | × 1e6 | NEATS application, before PyNEATS is called | **confirmed by the task owner**; application source not inspected | `NEATSFuel` (`neats_fuel.py:92`), the BADA correction ratio (`bada_model.py:841,875`), SAC/CoCiP through `Fuel.q_fuel` and the fleet column `q_fuel`. All treat it as J/kg. No × 1e6 exists in `src/` | + +Direct callers of PyNEATS must supply m/s and J/kg themselves. This covers the FleetRunner +JSON path, `scripts/run_json.py` and DataFrame users. Two places get no upstream protection: + +- **Airport fuel file.** PyNEATS reads `airport_fuel_path` directly, so a `q_fuel` column in + MJ/kg would enter unconverted (A10). +- **No magnitude check.** PyNEATS does not check the size of `q_fuel`. On the BADA path, + 42.8 (MJ/kg) fails loudly through the 100 kg/s fuel-flow outlier check. When `ff`, `ee` and + `am` are all provided, BADA is skipped and the value reaches SAC/CoCiP silently. + +`scripts/schemas/flight_schema.json:106-109` labels `calorific_value` as MJ/kg. The golden +inputs it validates contain J/kg values, and the code expects J/kg (Section 9). + +## 6. Ambiguities and unverified conventions on the AO input path + +### 6.0 What the regulation and official documents settle + +Investigated after the code audit. Each field the AO schema leaves undefined was traced to +its authoritative definition, and two independent checks then re-opened every cited source. +Sources were ranked in this order: +1. EU law: Implementing Regulation (EU) 2024/2493, which inserts Articles 56a and 56b and + Annex IIIa into the MRR, Implementing Regulation (EU) 2018/2066. The English Annex IIIa is + unchanged in the consolidated MRR of 27 May 2025. The later act 2025/842 and the + corrigenda affect other language versions only. +2. Official Commission and NEATS documents: + - the Commission DG CLIMA "MRV non-CO2 data collection Guidance", v2.1, February 2026; + - the DLR/DWD "Reference set of technical specifications for the MRV" (RSTS, + EC-CLIMA/2024/NP/0014, March 2025); + - the pre-NEATS "System Architecture Description" (SAD, 10 January 2025). This is a draft + input specification that the AO schema has since superseded. +3. Standards: ICAO Doc 4444, ICAO Annex 16 Vol. II, DEF STAN 91-091 / ASTM D1655. +4. Repository artefacts. + +"§" below refers to MRR Annex IIIa. + +| AO field | Authoritative definition | Source | PyNEATS | Remaining gap | +| --- | --- | --- | --- | --- | +| `ff` | "the mass of fuel in kilograms that passes through the aircraft fuel system and into the aircraft's engines per second": **whole aircraft, kg/s** | §1(33); SAD "total fuel mass flow rate"; RSTS Table 2.13 (CoCiP input) "total" | matches | RSTS Tables 2.5, 2.7, 2.8 and 2.10 say "fuel mass flow rate of each engine". The law never says "total" or "all engines" in so many words; the whole-aircraft reading is its plain meaning | +| `ee` | "the percentage of useful thrust generated by an aircraft engine relative to the energy input from fuel". The guidance calls it "aircraft overall propulsion efficiency"; the SAD and RSTS give "Overall propulsion efficiency", unitless, [0, 1] | §1(34), §1(35) ("by all-time stamps"); guidance §3.3; SAD II.11.2.7 | matches (F·TAS / (ff·LHV), pycontrails SAC) | The law says "percentage" while the schema, SAD and RSTS give [0, 1]; the schema's maximum of 1 rejects percentage input. Net or gross thrust is not stated. The guidance wording could also be read as one value per flight | +| `load_factor` | "the weight of passengers, cargo and baggage, including mail and hand luggage, expressed as fraction of the maximum payload mass"; maximum payload "can be retrieved by the applied fuel burn module"; default 1 | §1(31), §1(32), §5(2)(b); SAD and RSTS "fraction of the maximum payload mass, e.g., as defined in BADA" | matches (ZFW = OEW + pf·MPL) | none for the definition. §1(3) also uses "load factor" in the aerodynamic (g-load) sense, so the template should name the payload meaning explicitly. For BADA3, PyNEATS uses MTOW − OEW instead of MPL (I4) | +| `hydrogen_per_carbon_ratio` | "the number of hydrogen atoms (H) per carbon atom (C) per molecule of the fuel": **atomic (molar)** | §1(36); ReFuelEU Union Database field "Molar ratio of hydrogen to carbon" | matches | The SAD draft said "ratio of hydrogen mass to carbon mass"; the law prevails. `mrv-specification.md` still says "molecular mass ratio" (Section 9) | +| `hydrogen_content` | not in the law. NEATS default "Hydrogen Content 13.79 % mass"; RSTS and SAD "Percentage of hydrogen mass content" | guidance v2.1 §3.3 Table 3; RSTS Table 2.12 | matches (wt.%) | none. PyNEATS's own `input-format.md` and `flight_schema.json` say "fraction (0–1)" (Section 9) | +| `calorific_value` | "net calorific value of the fuel on board", which is the lower heating value (MRR Art. 3 NCV definition); MJ/kg in the schema, SAD and ReFuelEU; default 42.8 MJ/kg | §1(38); MRR Art. 3; guidance Table 3 | matches (J/kg after the application's × 1e6) | none | +| `takeoff_mass` | "the aircraft mass at beginning of the take-off run, including everything and everyone carried at that moment, in kilograms" | §1(29) | unit matches. PyNEATS applies it at the first retained trajectory point, which may be a taxi point (A6) | The SAD says "actual", while the guidance names the operational flight plan as the source | +| `am` | "the mass in kilogrammes of the aircraft along the trajectory, which equals to subtracting from the take-off mass the fuel burn" | §1(28) | matches | none | +| `aromatic_content`, `sulphur`, `naphthalene` | Basis: defaults "25 % volume", "0,3 % mass", "3,0 % volume" (the ASTM Jet A-1 limits); ReFuelEU reports "aromatics and naphthalenes by percentage volume and of sulphur by percentage mass". Sulphur and naphthalene are not defined in §1 and appear only as defaults. Not used by NEATS yet | §5(4); §1(37); ReFuelEU Reg. (EU) 2023/2405 Art. 10(d); guidance Table 3 footnote | not used. The defaults 0.25 / 0.003 / 0.03 are fractions on these bases | **Scale (25 or 0.25) is not defined** by any public source. The volume basis of the AO aromatics field is inferred from the defaults, and the SAD draft said "mass". The ReFuelEU Union Database uses a percent scale ("20.5 % is reported as 20.5") | +| `departure_date_time` | The law says only "day and time of departure … UTC" (§1(1)). The schema ties the field to flight-plan Item 13, which holds "the estimated off-block time (EOBT)", never an actual time | ICAO Doc 4444 App. 2, Item 13 | label mismatch: PyNEATS stores it as `aobt` and calls it "actual off-block time" (no numbers change) | EOBT or AOBT is not fixed in law. Item 13 carries only HHMM; the date comes from Item 18 DOF/ | +| `arrival_date_time` | "day and time of … arrival, UTC"; the event (landing, in-block, EOBT + EET) is undefined everywhere | §1(1) | dropped by PyNEATS | the AIBT label in `flight_schema.json` has no basis | +| `ts` | time stamp in seconds; interval "not exceeding 60 seconds"; UTC. The RSTS asks for "YYYY-MM-DDTHH:MM:SS" | §1(6), §1(7); RSTS Table 2.1 | matches for UTC input (D2 for offsets and fractional seconds) | none | +| `lat`, `lon`, `fl` | latitude and longitude "in decimal degree"; altitude "in pressure altitude"; the SAD gives the FL number in hundreds of feet | §1(5); SAD | matches | none | +| `tas` | "the speed of the aircraft relative to the air mass …, in meters per second (m/s)". The AO schema and RSTS §2.2.3 (operator-supplied data) say knots | §1(4); RSTS §2.2.3 | PyNEATS uses m/s (the legal unit); the application converts the schema's knots | Not a defect for uploads that follow the schema. An operator who follows the regulation's m/s instead would have TAS scaled by 0.514 without any error; a plausibility check upstream would catch it | +| `nwp_date_time` | The law requires a "common reference NWP model, provided by the Commission through NEATS" for enhanced weather data. Method C uses "current weather information" usable "for individual flight planning", i.e. forecast data. The guidance says "Weather data is always provided by NEATS as Secondary data". The RSTS says "flights are evaluated based on the latest weather forecast prior to departure" (ICON runs at 00/06/12/18 UTC, +61 h) | §1(25), §1(26), §4; Art. 56a(7)(b), 56b(7); guidance §2.3, §3.3; RSTS §2.4.5, §4.2 | not read; weather is chosen by the caller | No public source defines the field (run initialisation time or valid time; behaviour when null). The run rule rests on a single non-binding RSTS sentence. Art. 56a(4) also allows AO or third-party weather for basic weather data | +| NOx basis (internal, for O3) | ICAO Annex 16 Vol. II: NOx "calculated as if the nitric oxide were in the form of nitrogen dioxide", "EI(NOx) as NO2". OpenAirClim v0.13 response surfaces are per kg NO, and `nox = "NO2"` applies 30/46; v0.19.0 (2026-09-28) made NO2 the default | ICAO Annex 16 Vol. II Ch. 1, App. 3 §7.1.2; EASA EEDB introduction; openairclim `calc_response.py` | **mismatch** (O3): `nox = "NO"` | No MRV text prescribes the OpenAirClim setting. The method owner should confirm in writing | + +Consequences for this audit: +- A1, A2, A3 and A4 are no longer open questions about the template. The regulation defines + each one, and PyNEATS follows the definition. The remaining action is to align the AO + template wording with the law: "whole aircraft", "overall propulsion efficiency as a + fraction", "payload ÷ maximum payload, not seat load factor", "atomic H/C". +- A5 narrows to the numeric scale. A7 becomes a labelling correction (EOBT). A8 stays open + for the NEATS application team. + +Each entry gives: +1. the field and code path; +2. the expected and observed convention, with confidence; +3. the consequence or the open question; +4. whether the application might already handle it; +5. a recommendation. Recommendations are not implemented. + +**A1. `ff`: whole aircraft or per engine.** *Resolved (§6.0): whole aircraft, MRR Annex IIIa §1(33). PyNEATS matches. The analysis below explains the consequences if an operator nevertheless sends per-engine values.* +1. `trajectory_data[].ff` → `neats_io.py:79` → `bada_model.py:693,1071-1076` → pycontrails + Emissions and CoCiP. +2. The AO schema says only "kg·s⁻¹ at the data point". Every PyNEATS consumer treats the + value as whole-aircraft (high confidence): + - pycontrails divides by `n_engine` to get per-engine flow (`emissions.py:183-195`); + - CoCiP uses `fuel_flow / true_airspeed` as fuel per metre (`cocip.py:929-931`); + - the fuel-burn guardrail compares Σ ff·Δt with MTOW − OEW; + - every internally produced `ff` is whole-aircraft: BADA, and −Δam/Δt. +3. A per-engine `ff` would make fuel burn, CO2, every emission mass and the CoCiP fuel per + metre n_engine times too small (× 0.5 for a twin, × 0.25 for a quad). Thrust-setting EIs + would be looked up at 1/n_engine of the real flow. No guard catches this: + - the guardrail rejects only overestimates; + - the 100 kg/s outlier check never inspects a provided `ff`; + - when `ff`, `ee` and `am` are all present, every check is skipped (`bada_model.py:471-477`). +4. Possibly. The application might normalise `ff`, but nothing confirms it. +5. **Clarify the AO template**: "total fuel flow of the aircraft, all engines, kg·s⁻¹". No + conversion is needed if AOs send the total. Optionally, apply the existing fuel-burn + guardrail on the early-exit path as well, so a kg/h value fails loudly. + +**A2. `ee`: which efficiency.** *Resolved (§6.0): overall propulsion efficiency, MRR §1(34) and Commission guidance §3.3. PyNEATS matches. Open: "percentage" in law against [0, 1] in the schema.* +1. `trajectory_data[].ee` → `neats_io.py:80` → `bada_model.py:1046-1048` (kept) → pycontrails + `sac.slope_mixing_line` (`sac.py:206`). +2. The AO schema says "engine efficiency, unitless, [0, 1]". The consumer needs **overall + propulsion efficiency** η = F·TAS / (ff·LHV), which is what PyNEATS computes when `ee` is + absent (`jet.py:178`), clipped to [0, 0.5]. Confidence is high for the consumer side. + `examples/README.md:140` calls the field "engine thermal efficiency". +3. The SAC slope scales with 1/(1 − η). A thermal efficiency of about 0.5 against an overall + efficiency of about 0.35 raises G by a factor of 1.3. `ee` = 1.0 is schema-valid and makes + G infinite. Provided values are never clipped. +4. Unknown. +5. **Clarify the AO template**: "overall propulsion efficiency = thrust × TAS / (fuel flow × + LHV), fraction, typically ≤ 0.5". Optionally clip provided values in PyNEATS to the + [0, 0.5] range already used for computed values. + +**A3. `load_factor`: basis.** *Resolved (§6.0): payload mass ÷ maximum payload mass, MRR §1(32). PyNEATS matches.* +1. `load_factor` → `payload_factor` (`neats_io.py:113`) → `_estimate_initial_mass_iterative` + (`bada_model.py:963-980`). +2. The AO schema says "the load factor, unitless". PyNEATS uses it as a payload-mass fraction + of the BADA maximum payload, ZFW = OEW + pf·MPL. The first-pass guess uses pf·(MTOW − OEW). + `flight_schema.json:78-81` states "payload / max-payload ratio" (medium confidence that AOs + will use the same reading). It is used only when neither `am` nor `takeoff_mass` is given. + The default is 1.0, meaning full payload. Values above 1 are not clipped. +3. An AO sending a passenger load factor (for example 0.85) would get 0.85 × MPL of payload. + That differs from the real payload of passengers, baggage and cargo, so initial mass and + fuel change. For BADA3 aircraft the basis is MTOW − OEW rather than MPL (I4). +4. Possibly. The application might convert a seat load factor. +5. **Clarify the AO template**: "payload mass ÷ maximum payload, fraction 0–1; not a seat load + factor". + +**A4. `hydrogen_per_carbon_ratio`: atomic or mass basis.** *Resolved (§6.0): atomic, MRR §1(36). PyNEATS matches.* +1. `neats_io.py:116` → `neats_fuel.py:83-87`. +2. The AO schema says "unitless". PyNEATS assumes an atomic (molar) ratio. The DLR range + [1.9061, 2.1857] equals [13.79, 15.50] wt.% exactly under that formula. Confidence that + PyNEATS intends atomic is high; confidence that AOs will send atomic is unknown. + `mrv-specification.md:233-238` calls it a "molecular mass ratio" and leaves out the × 100. +3. Only matters when `hydrogen_content` is absent. A mass ratio (about 0.16 for Jet A) would + give about 1.3 wt.% H, so `ei_h2o` would be about 10 times too low. That affects H2O + emissions, CoCiP SAC, the OpenAirClim H2O inventory and the nvPM SAF correction. There is + no range check (test case `TC_HC_RATIO_OUT_OF_RANGE` is a strict xfail). +4. Possibly. +5. **Clarify the AO template**: "atomic (molar) H/C ratio, about 1.9–2.2 for Jet A". Correct + the MRV spec formula (Section 9). + +**A5. `aromatic_content`, `sulphur`, `naphthalene`: scale and basis.** *Partly resolved (§6.0): basis from the MRR §5(4) defaults and ReFuelEU (aromatics and naphthalene by volume, sulphur by mass). The JSON scale is still undefined.* +1. `neats_io.py:117-120` → `neats_fuel.py:94-102`. +2. The AO schema says "unitless" with no scale and no mass/volume basis. The PyNEATS defaults + are fractions (0.25, 0.003 commented as "0.3 %", 0.03). The airport test fixture looks like + percentages (19.1, 0.07, 0.36). `flight_schema.json` says "mass fraction". None of these + values enters a calculation. Naphthalene never reaches the fuel object (I2). +3. No numerical effect today. Any future use (for example sulphur-scaled `ei_so2`) would face + an undefined scale. The SO2 and sulphate outputs use the fixed Jet A-1 EI (600 ppm S basis) + whatever sulphur is supplied. +4. Possibly. +5. **Clarify the AO template** (fraction, % or ppm; mass or volume) before any use. Document + that SO2 and sulphate outputs are fixed Jet A-1 values. + +**A6. `takeoff_mass`: point of application and precedence.** *Definition settled (§6.0): mass at the beginning of the take-off run, MRR §1(29). PyNEATS's point of application does not guarantee that event (see below).* +1. `takeoff_weight` → `_choose_mass_strategy` / `_single_pass_performance` + (`bada_model.py:786-808,849`). +2. The unit is kg and consistent. It is applied at the first **retained** trajectory point: + - after resampling to 1 min; + - for BADA3 aircraft, after the speed filter; + - on the altitude-fallback path, at the first point at or above FL20. + It is ignored, with only a debug log, when an `am` column exists, and nothing checks it + against `am[0]`. It is not clipped to MTOW. +3. Ground or taxi points before rotation are charged against it. On the FL20 fallback, the + climb fuel below FL20 is never subtracted. Conflicting `am` and `takeoff_mass` pass without + notice. +4. Possibly, if the application trims the trajectory to the take-off point. +5. **Clarify the AO template**: mass at brake release or at the first trajectory point. + Optionally log at warning level when `am` overrides `takeoff_mass`. + +**A7. `departure_date_time`: EOBT or AOBT.** *Partly resolved (§6.0): Item 13 holds the EOBT (ICAO Doc 4444). The MRR does not choose. The PyNEATS "AOBT" label is unsupported.* +1. `neats_io.py:106` → attr `aobt` → `report.py:114-135`. +2. The AO schema says "Item 13 of the flight plan in UTC", which is the EOBT. PyNEATS docs + (`flight_schema.json:38`, `output-format.md:50`) label it "actual off-block time (AOBT)". + The string is copied verbatim and never parsed. +3. No numerical effect. The output label may describe the wrong quantity. Flight keys compare + raw strings, so the same instant written as `Z` and as `+0000` gives two different keys. +4. Unknown. +5. **Clarify the AO template** and PyNEATS output docs. If the field stays, normalise it to a + UTC ISO string at the boundary. + +**A8. `nwp_date_time`: who selects the NWP instance.** *Still open (§6.0): weather comes from the common NEATS NWP model. The RSTS rule is "latest weather forecast prior to departure". The field itself is not defined publicly.* +1. Not read anywhere in PyNEATS. Weather comes from `zarr_paths` (`fleet.py:578-584`) or + `DWDFactory(...)(asofdate, hour)` (`weather_factory.py:259-296`). +2. The AO schema requires each flight to be computed with the named NWP instance. +3. A fleet that mixes flights requesting different NWP instances runs on one store. Flights + outside the store's time range get NaN meteorology, and NaN wind is set to 0 + (`weather_provider.py:382-383`). `DWDFactory` builds its window with `strftime` and no time + zone conversion, so `asofdate` must be naive UTC or UTC-aware. +4. Likely. The application may group flights by `nwp_date_time`. +5. **Confirm with the application owners.** At minimum, the `FleetRunnerParams` docstring now + states that the field is not read. + +**A9. DataFrame entry contract.** +1. `split_df_into_flights` / `FlightRunner(source=df)` → `NeatsTrajectoryParser`. +2. `altitude` must be a flight level and `time` must match `date_format`. Until this patch the + user-facing entry points did not say so. +3. Altitude supplied in metres becomes 30.48 times too high. Time strings in ISO form (with + `T`, an offset or fractions) become NaT. Those rows are dropped silently, and the flight + fails only if every row is dropped. +4. Not applicable (direct library use). +5. Documented in the patch. Optionally rename the column or check its range. + +**A10. Airport fuel file units.** +1. `_load_airport_fuel_props` (`neats_parser.py:44-120`) → fallback (`neats_parser.py:284-296`). +2. Values are passed through `float()` only. The tests use `q_fuel` in J/kg. `h_c_ratio` is + treated as atomic. There is no `hydrogen_content` column. The file's own units are + unverified, and no production file is in the repository. +3. A `q_fuel` column in MJ/kg would enter 1e6 times too small, because the application's + conversion does not apply to this file. A mass-based `h_c_ratio` would make `ei_h2o` about + 10 times too low. +4. No. PyNEATS reads the file itself. +5. **Confirm units with the file owner.** The patch documents the assumed units in the loader + docstring. Optionally add a magnitude guard for `q_fuel`. + +**A11. `q_fuel` correction on the default branch.** +1. `bada_model.py:565-571,764-772,841,875`. +2. `mrv-specification.md:220` says no correction is applied when the AO gives no `q_fuel`. The + code always applies the ratio 43.13/42.8, because `NEATSFuel` always carries a `q_fuel` + (default 42.8e6 J/kg). The `else reference_q_fuel` branch cannot be reached. +3. BADA fuel flow, fuel burn, CO2 and every fuel-proportional quantity are +0.77 % for flights + without `calorific_value`. The iterative take-off mass shifts as well. +4. No. +5. **Decide in PyNEATS**: apply the correction only to a supplied `q_fuel`, or amend the spec. + +**A12. `REFERENCE_Q_FUEL` against BADA4 LHV.** +1. `neats_default_parameters.py:315-317` → `bada_model.py:841`. +2. A constant 43.13e6 J/kg (the pycontrails JetA value) stands in for the LHV that BADA fuel + flow corresponds to. BADA4 files carry their own LHV; pyBADA's dummy files have 43.217e6 + and 43.71e6 J/kg. BADA3 has none. +3. The correction is off by LHV_file / 43.13e6, which is +0.2 % to +1.35 % on the dummy files. + Licensed files were not checked. +4. No. +5. **Check against BADA4 files**. Consider reading the adapter's LHV for BADA4. + +**A13. Mixed fuel-flow provenance (documented in part).** +- When `ff` is provided without `am`: + - output `fuel_flow` and `fuel_burn` come from the AO; + - `aircraft_mass` and the iterative take-off mass are integrated from BADA fuel flow + (`bada_model.py:927,990-994,1018-1022,1067`); + - so mass(t₀) − mass(t_end) ≠ Σ fuel_burn. +- A computed `ee` always uses BADA fuel flow. That is documented design + (`mrv-specification.md:190`; `input-prioritization.md:44`). +- Units agree throughout. The issue is provenance, not scale. +- **Recommendation**: document this, or integrate mass from the provided `ff`. The decision + belongs to the method owner. + +## 7. New discrepancies + +### On the AO input path + +**D1. Bare single-flight object rejected.** +1. Top-level JSON → `neats_json_to_flights` (`neats_io.py:45-46`) → `FleetRunner._load_trajectories` + (`fleet.py:561-566`) and `scripts/run_json.py:238-240`. +2. Expected: the AO schema's `anyOf` allows one object. Observed: iterating a dict yields its + keys, so `'str'.get` raises `AttributeError` and aborts the whole run (high confidence, + snippet). The list-only contract is stated in the function's type and docstring and in + `input-format.md:3-11`. +3. A single-object upload fails completely. There is no numerical effect. +4. Likely, if the application always wraps uploads in a list. +5. **Convert in PyNEATS** (wrap a `Mapping` in a list at the file-reading entry points), or + document that the array form is required. + +**D2. Trajectory timestamp offsets discarded; some valid forms rejected.** +1. `ts` → `neats_io.py:61-68` → `neats_parser.py:233-242`. +2. Expected: every RFC 3339 `date-time` maps to the same instant in UTC. Observed: + - `strptime('%Y-%m-%dT%H:%M:%S%z')` followed by `strftime('%Y-%m-%d %H:%M:%S')` **drops** + the offset without applying it. `2025-11-02T12:02:00+02:00` becomes 12:02 UTC instead of + 10:02 UTC. + - Fractional seconds and a lowercase `z` are schema-valid but raise `ValueError`. + - Nothing catches the error, so it aborts the batch before per-flight error handling. + - `Z`, `+0000` and `+00:00` are handled correctly. Every `ts` in the tests uses `+0000`. + High confidence (snippets). +3. A non-zero offset shifts every waypoint by the offset. Meteorology and radiation are then + sampled at the wrong time, and CoCiP day or night can flip. A schema-valid fractional + second aborts every flight in the batch. Schema-conformant UTC input is unaffected. +4. Possibly, if the application normalises `ts` to whole seconds with `Z` or `+0000`. +5. **Convert in PyNEATS**: parse with `pd.to_datetime(..., utc=True, format='ISO8601')` or + `datetime.fromisoformat(...).astimezone(timezone.utc)`. Also **clarify the AO template**: + UTC only, whole seconds. + +**D3. Fuel flow derived from `am` is misaligned in time.** +1. `am` → `_derive_fuel_flow_from_mass_if_needed` (`bada_model.py:705-749`). It runs when `am` + is present and `ff` is absent. +2. Expected (docstring and comment): ff[i] = −(m[i+1] − m[i]) / Δt[i]. Observed: + `np.diff(mass, prepend=mass[0])` gives m[i] − m[i−1], divided by the **forward** + Δt[i] = t[i+1] − t[i]. So ff[0] = 0, the last value is NaN, and every rate sits one sample + late. High confidence. Units are consistent (kg/s, whole aircraft). +3. On a uniform 60 s grid, rates lag by one sample, and Σ fuel_burn = m[0] − m[n−2], so the + last segment's burn is lost. Where spacing is non-uniform (rows removed by the speed + filter, or a non-default `interpolation_time`), each rate is off by Δt[i−1]/Δt[i]. This + affects `fuel_flow`, `fuel_burn`, emissions and the guardrail total. It does not affect + `engine_efficiency`, which uses BADA fuel flow. +4. No. +5. **Convert in PyNEATS**: use `-np.diff(mass, append=np.nan) / segment_duration`, computed + before the speed filter or with `segment_duration` recomputed after it. Regenerate golden + files that use `am` without `ff`. + +**D4. Partially null optional columns are not filled.** +1. `tas`, `ff`, `ee`, `am` → `neats_io.py:93-94` → `pycontrails_interpolation.py:108-112` → + `bada_model.py:1046,1071-1076`. +2. Expected: the schema allows `null` per point, and such points should be filled with values + in the same convention, or the flight rejected. Observed: a column survives unless every + value is null. `np.interp` spreads NaN to resampled samples between neighbours. A provided + column is used wholesale, with no per-point BADA fill. NaN `ee` gives NaN SAC slope at + those points. High confidence. +3. Fuel, CO2 and emission totals (nansum) undercount the affected segments, and the guardrail + weakens by the same amount. On a 60 s input grid, only the null point itself is affected. + With 5 min input spacing, about 9 resampled samples are affected per null. +4. Possibly, if the application never sends partial nulls. +5. **Convert in PyNEATS**: fill NaN points of provided columns from the BADA result, or reject + partial columns. + +### Elsewhere in the inspected PyNEATS path (not AO fields) + +**O1. Local aCCF H2O: potential vorticity in SI, coefficients in PVU.** This is the most +consequential finding. +1. DWD `pv` → pycontrails `PotentialVorticity` (`weather_factory.py:242`) → `LocalACCFModel.compute_h2o` + (`local_accf.py:239-240`) → `h2o_raw_formula` (`local_accf.py:228-231`): + 4.05e-16 + 1.48e-16·|pv|. This is the **default** large-emitter non-CO2 model + (`neats_default_parameters.py:115`). +2. Expected: the CLIMaCCF reference implementation applies this regression to PV in PVU + (`climaccf/accf.py`: `self.pvu = ds['pv'].values * 1e6`, present since 2022). Matthes et + al. 2023, Annex A says "PV given in standard PV units [10⁻⁶ K kg⁻¹ m² s⁻¹]". Observed: the + column is in SI, K m² kg⁻¹ s⁻¹ (pycontrails units; golden |pv| about 1e-7 to 1.3e-5), and is + used unconverted. High confidence: golden ATR_20_H2O / fuel_burn × 3 equals 4.05e-16 on + every valid sample, which means the weather has no effect. +3. The missing factor is 1e6 on `pv`. The per-sample H2O aCCF is too small by + 1 + 0.365·|PV_PVU|: about 1.9 times at the golden median and up to about 5.9 times. + Fuel-weighted per flight, ATR_20_H2O, H2O EAGWP and CO2eq are underestimated 1.59–3.26 + times. EAF2143 H2O CO2eq_20 would rise from about 764 kg to about 2,500 kg. The ClimAccf + backend does not have this error, because CLIMaCCF converts internally, so the two + backends disagree. +4. No. PV comes from PyNEATS's own weather path. +5. **Convert in PyNEATS**: multiply by 1e6 inside `compute_h2o`. Do **not** convert in the + weather layer, because the ClimAccf backend needs SI. Regenerate golden outputs, and add a + LocalACCF-against-CLIMaCCF regression test. + +**O2. Method D `contrails_ef_J` holds an AGWP, not joules.** +1. `gwp.py:426`: `total_ef_J = float(flight.attrs.get("AGWP_20_CONT", 0.0))`, small emitter only. +2. Expected: contrail energy forcing in J, as in Method C (`gwp.py:411`). Observed: the raw + OpenAirClim 20-yr contrail AGWP in W m-2 yr, with no efficacy and always horizon 20. High + confidence. +3. A different physical quantity sits under a J-named key. It cannot be converted to energy + forcing by a fixed factor. Summing or comparing `contrails_ef_J` across Methods C and D + mixes units. CO2eq is not affected. +4. No. +5. **Change in PyNEATS**: emit null for Method D, or use a separately named field. + +**O3. OpenAirClim told the NOx species is NO.** +1. pycontrails NOx mass (`emissions.py:263-273,773`) → inventory `NOx` in kg + (`open_airclim.py:131,181`) → `src/pyneats/resources/openairclim.toml:10` + `nox = "NO"` → openairclim `calc_response.py:106-133` (`corr_nox` = 1.0). +2. Expected: the inventory mass is NO2-equivalent. This is the ICAO EDB convention, the + pycontrails APCEMM interface documents `nox_ei` as kg(NO2)/kg (`models/apcemm/inputs.py:70`), + and the CLIMaCCF literature uses it. The TOML should then read `"NO2"`, so that openairclim + applies 30/46 to its per-kg-NO response surfaces (`resp_RF_O3.nc` title "NO Emission → O3 + RF"). Observed: `"NO"`, carried over verbatim from the upstream openairclim demo + `example.toml`. Medium-high confidence: the NO2 basis is a documented convention, not + demonstrated by computation. +3. Method D O3 RF and AGWP are exactly 46/30 = 1.533 times too high. CH4 and PMO are about + 1.533 times too high. H2O, CO2 and contrails are unaffected. The local aCCF path is + consistent: K/kg(NO2) times an NO2-basis EI. +4. Unlikely. The application could pass its own `base_toml_path`. +5. **Change in PyNEATS**: set `nox = "NO2"` (or force it in `_create_worker_toml`) after + confirming the EDB basis with the method owner. Do not convert `nox_ei`, because that would + break the aCCF path. + +**O4. OpenAirClim metric window pinned to 2025.** +1. Year of the first waypoint → `_create_worker_toml` (`open_airclim.py:94`) sets + `time.range = [Y, Y+101]`. The bundled TOML keeps `metrics.t_0 = [2025]` + (`openairclim.toml:98`), and the harvest reads `AGWP_{h}_2025` (`open_airclim.py:251`). +2. Expected: AGWP integrated from the emission year. Observed: + - Y > 2025: openairclim `_check_metrics` raises `ValueError` (`read_config.py:681-690`), + and the flight becomes an error record. This covers every 2026 flight. + - Y ≠ 2025 with the default base inventory: `calc_cont.check_cont_input` also raises. + - Y < 2025 without a base inventory: silently wrong windows. + High confidence for the failure. All golden flights are from 2025. +3. With the shipped configuration, Method D works only for flights dated 2025, unless the + deployment overrides the TOML. +4. Possibly, if the application overrides the TOML. Not visible. +5. **Decide in PyNEATS** (method owner): anchor every flight to the 2025 pulse convention + consistently (inventory year, base inventory, CO2 coefficients), or set `t_0` and the + inventory year per flight. + +**O5. DWD `sdr` is TOA, and is scaled by 3600 twice on the ClimAccf path.** +1. `weather_factory.py:247` maps DWD `sdr` to `surface_solar_downward_radiation`. `build_cache` + multiplies by `sdr_accumulate_dt_s` (default 3600) and labels the result J m-2 + (`weather_factory.py:405-411`). `make_accf_surface_view` relabels it W m-2 without + converting (`climaccf.py:88-93`), and pycontrails ACCF multiplies by 3600 again + (`models/accf.py:405-423`). +2. Expected: surface downward solar radiation, accumulated once in J m-2. Observed: + - The store history (`swflx_dn`, "solar incoming radiation at TOA") shows it is **TOA** + incoming solar. + - With the default cache settings it is scaled by 3600 twice. + - The live path and the example cache (dt = None) are scaled once. + High confidence in the inconsistency. +3. Only the non-default `accf` backend is affected. Whether CLIMaCCF's consumed aCCFs (O3, CH4, + H2O) read `ssrd` is unverified. CoCiP is not affected: it computes its own TOA insolation. +4. No. +5. **Change in PyNEATS**: keep one stored convention with a matching units attribute, make the + adapter convert rather than relabel, and stop calling TOA flux "surface". + +**O6. DWD `rhi` relabelled from % to 1 without dividing by 100.** +1. `weather_factory.py:522-523`; `_standardize_vars` at `weather_factory.py:588-603`. +2. Expected: values in % are divided by 100 when relabelled to a fraction. Observed: only the + label changes. `clc`, two lines away, is divided by 100. The source units of the DWD files + are unverified. +3. If the source is in %, RH is 100 times too large for its only consumer, the non-default + pycontrails ACCF. +4. No. +5. **Change in PyNEATS**: divide by 100 when the source attribute is `%`. + +**O7. OLR renamed to ECMWF `ttr` without a sign change.** +1. `climaccf.py:74-87`. +2. Expected: ECMWF defines ttr = −OLR (pycontrails `datalib/ecmwf/variables.py:174-188`). + Observed: a rename only, so values stay positive. +3. Non-default backend only. The effect inside CLIMaCCF is unverified. +4. No. +5. Negate on rename, unless the CLIMaCCF source shows it expects positive OLR. + +**O8. Weather factory relabels units without converting.** Ambiguity rather than a +demonstrated error. +- `_standardize_vars` overwrites `attrs['units']` with the pycontrails unit after a warning, + without arithmetic (`weather_factory.py:572-608`). +- The module docstring claimed "automatic unit conversions" (corrected in the patch). The + only real conversions are: + - `level` from Pa to hPa (a heuristic on the first element); + - `clc` divided by 100; + - SDR multiplied by dt in `build_cache`. +- Any source variable in non-pycontrails units would be mislabelled silently. Source units of + the DWD files were not available. +- **Recommendation**: raise on mismatch, or convert. + +## 8. Further ambiguities in outputs and methodology + +- **O9. Low-altitude zero results.** A flight rejected for never reaching 550 hPa reports + NaN fuel and CO2 baselines but 0.0 EAGWP and CO2eq for every species, **CO2 included** + (`fleet.py:295-319`). This is intentional and pinned by tests + (`tests/test_zero_climate_result.py`). Summing `CO2eq_kg` across flights undercounts + CO2. **Document** the semantics in `output-format.md`, or reconsider zeroing CO2. +- **O10. Method C silent species drop.** A conversion exception or a missing ATR column drops + a species without an error record (`gwp.py:216-217,239-240`). Partial overrides of the + k-factor, efficacy or `rf_backward` tables fall back to a scale of 1.0 (`gwp.py:155-165`). + An H2O entry missing from `k_atr` would be off by a factor of 31. This is reachable only with + non-default parameters. +- **O11. `FlightView.to_dict` encodings.** `time` is integer Unix seconds, truncated rather + than rounded as the old comment said. `contrail_age` is integer **nanoseconds**. Both are + now documented in the docstring. They apply to golden-file dumps, not to `runner.results`. +- **O12. Performance-only output.** + - `true_airspeed` is in m/s, not the AO's knots. + - `timestamps` carry no offset (naive UTC); the docstring example showed `+00:00` and was + corrected. + - Multi-engine copies are not merged. + - **Document**, or emit `Z`. +- **O13. Totals limited to the met domain.** `total_nox`, `total_co`, `total_hc` and + `total_nvpm_*` sum only waypoints inside the met domain (default 550–140 hPa). EIs are NaN + elsewhere. `total_co2` and `total_h2o` cover every segment. For EAF2143, 31.8 % of the + flight's fuel lies outside the domain. The OpenAirClim NOx inventory fills those NaN with 0, + while CO2 and H2O stay complete, so one inventory mixes two bases. **Document** the basis, + or decide on fallback EIs. +- **O14. PMO merged into CH4.** PMO takes CH4 efficacy and conversion factors in both methods + (`open_airclim.py:264-266`; `local_accf.py:221-222`). `EFFICACY['PMO']` and the PMO + conversion-factor columns are unused. The effect is about +0.9 % on the combined CH4+PMO + EAGWP in Method D. **Clarify** with the specification owner. +- **O15. ClimAccf `unit_K_per_kg_fuel=True`.** CLIMaCCF would substitute its own fleet-mean + NOx EI (0.013), or apply no conversion when `merged=False`. The default is False. **Reject** + True in `ACCFParams`, or document it (done in the patch). +- **O16. RF-backward coupling.** `RF_BACKWARD_FACTOR` is approximately + 1/`ACCF_SCALE_*['V1.0']`. That is documented intent (`physics.py:92-97`; + `mrv-specification.md:310`) with at most 0.16 % rounding mismatch. It is not a discrepancy, + and the coupling is now stated in a comment. + +## 9. Documentation and schema inconsistencies (reported, not edited) + +The patch corrected wrong unit statements in **code comments and docstrings**: the items +marked "corrected" below and those listed in Section 11. Markdown documents and JSON schemas +were left unchanged on purpose, because the scope excludes schema changes and fixes to +discrepancies: + +| Location | Statement | Code or demonstrated convention | Action | +| --- | --- | --- | --- | +| `docs/reference/input-format.md:43`; `scripts/schemas/flight_schema.json:94-96` | `hydrogen_content` is a "mass fraction (0–1)" | wt.% (13.79). A caller sending 0.138 gets `ei_h2o` 100 times too low | fix docs | +| `scripts/schemas/flight_schema.json:106-109` | `calorific_value` in MJ/kg | J/kg (golden inputs validated by this schema hold 4.1944e7) | fix schema text: PyNEATS input is J/kg after the application's × 1e6 | +| `docs/explanation/mrv-specification.md:233-238` | H/C is a "molecular mass ratio"; formula without × 100 | atomic ratio; × 100 to give wt.% | fix spec | +| `docs/explanation/mrv-specification.md:220` | no BADA correction when `q_fuel` is not provided | correction always applied (A11) | decide, then align | +| `docs/explanation/mrv-specification.md:264,277,302` | EAGWP and AGWP in W m-2 | W m-2 yr | fix spec | +| `docs/explanation/mrv-specification.md:272` | s_yr = 31,536,000 s | 31,556,952 s (`physics.py:58`); contrail EAGWP differs by 0.066 % | decide which year length, then align | +| `docs/explanation/mrv-specification.md:284-286` | C(H) = 24.16 / 47.57 / 74.36 e-15 | 24.12 / 47.50 / 74.24 e-15 (commit `bf3aa7b`, newer); CO2eq differs by about 0.15–0.17 % | update spec | +| `docs/reference/output-format.md` | per-segment `flight_results`; `climaccf_version`; error-record keys; "fleet totals" | `runner.results` holds flight-level entries; the keys differ; `total_*` are per-flight; no units are given | rewrite with a units table | +| `docs/reference/input-format.md:85-113` | minimal example at FL 4 and FL 10 | rejected by the 550 hPa guard | use cruise levels | +| `docs/howto/usage.md:85-90,112,168-169`; `docs/reference/example-scripts.md:11` | `ZarrPaths(met_path=…)`, `DWDFactory(...).build()`, "fuel flow & thrust" | `met_store`/`rad_store`, `__call__(asofdate, hour)`, no thrust output | fix docs | +| `examples/README.md:140` | `engine_efficiency` = "engine thermal efficiency" | overall propulsion efficiency | fix docs | +| `scripts/schemas/flight_schema.json:102-116`; `docs/reference/input-format.md` | aromatics and naphthalene as "mass fraction" | volume basis per MRR Annex IIIa §5(4) and ReFuelEU (the pre-NEATS SAD draft said mass); scale undefined | align with the law once the scale is decided | +| `scripts/schemas/flight_schema.json:38,43`; `docs/reference/output-format.md:50`; `report.py` | departure = "AOBT", arrival = "AIBT" | AO schema: Item 13 = EOBT; arrival event undefined | relabel neutrally ("departure time as supplied") | +| `docs/internal/performance_use_cases_analysis.md` | TAS always recomputed; early exit on FF+EE; FF never derived from mass | TAS kept when provided; early exit needs FF+EE+AM; FF derived from AM | mark superseded | +| `scripts/parquet_to_neats_json.py` | treats `altitude` as feet and stamps `+0000` without converting | basis (barometric or geometric) and time zone unverified | document expected input | +| `tests/test_performance_fallback.py:4,74,80`; `tests/test_guardrails.py:133,142,151` | FL15; 500 hPa | FL20; 550 hPa | update test docstrings | + +## 10. Incidental observations outside the unit scope + +These affect how the audited fields are consumed but are not unit issues. They are reported +for a separate decision and were not fixed. + +- **I1. BADA4 speed filter never runs.** `BADA4Adapter.v_stall_cas` calls + `VStall(mass=, config=)`, but the pyBADA 0.1.12 BADA4 signature is + `VStall(mass, HLid, LG, nz=1.0, **kwargs)` (`bada_adapters.py:175-180`; pyBADA + `bada4.py:2349`). The `TypeError` is swallowed (`bada_model.py:463-469`), so taxi and + low-speed points stay in for BADA4 aircraft. The unit tests mock `VStall`. +- **I2. Naphthalene key mismatch.** `neats_io` emits `naphthalene`, `Flight4D` allows + `naphtalene`, and `NEATSFuel.from_attrs` reads `naphthalene`. The AO value and the airport + value both fail to reach the fuel object. +- **I3. Zero values replaced.** `NEATSFuel.from_attrs` uses `or` between alternative keys, so + an explicit `sulphur_content` of 0.0 becomes the default 0.003 on the fuel object. +- **I4. BADA3 MPL hidden.** `BADA3Adapter.MPL` returns `None` with the comment "BADA3 has no + MPL attribute". pyBADA 0.1.12 `Bada3Aircraft` does set MPL (OPF maximum payload × 1000 kg, + `bada3.py:396-405`). BADA3 flights therefore estimate ZFW with pf × (MTOW − OEW), which is + heavier than pf × MPL. The comment was corrected; the behaviour was not changed. +- **I5. Parser `timezone` parameter has no effect.** pycontrails stores time as naive UTC + (`vector.py:2162-2167`). The docstrings were corrected. +- **I6. Outputs not strict JSON.** `scripts/run_json.py` writes no output. Golden JSON + contains bare `NaN` and `Infinity` tokens (28,945 and 318 in `fleet_5_flights_output.json`). +- **I7. `fleet_to_flights` docstring.** It said fuel is rebuilt from the Fleet columns; it is + rebuilt from attrs (corrected). +- **I8. Humidity-scaling adapter fallbacks.** The fallback branches rebuild `Flight` without + its fuel (`weather_provider.py:139-147`). They are unreachable with stock pycontrails. + +## 11. Patch summary and checks performed + +**Files modified.** Comments and docstrings were added or corrected in 27 source files (+826/−178 lines). New +file: `docs/internal/pyneats_units_audit.md` (this report). + +- Input boundary and parsing: `steps/parsing/neats_io.py` (a field-by-field unit mapping in + the `neats_json_to_flights` docstring, the two upstream-conversion comments, the `ts` + offset behaviour), `neats_parser.py` (FL → m pressure altitude; UTC parsing; corrected the + no-op `timezone` docstring; airport-file units), `views.py` (units of every Flight4D + column and attr), `params.py` (corrected 500 hPa → 550 hPa ≈ FL160). +- Fuel and constants: `core/neats_fuel.py` (wt.%, atomic H/C, J/kg template, fixed JetA + EIs), `core/neats_default_parameters.py` (units of every constant; corrected FL15 → FL20 + and 1.1/10 % → 1.2/20 %; `REFERENCE_Q_FUEL` is the pycontrails JetA value, not a BADA + value), `core/physics.py` (year length and C(H) against the spec; conversion-factor + direction and unverified absolute units), `core/fleet_utils.py` (fuel rebuilt from attrs), + `core/views.py` (`to_dict` time and timedelta encodings), `core/compute_parameters.py`. +- Performance: `steps/performance/bada_model.py` (parameter units and meaning; TAS template; + q_fuel correction; the mass-derived fuel-flow comment now describes the actual backward + difference; fuel-flow provenance), `bada_adapters.py` (adapter units; the BADA4 `VStall` + signature; the corrected BADA3 MPL comment), `speed_filter.py`, `altitude_filter.py`, + `steps/interpolation/pycontrails_interpolation.py`. +- Climate and outputs: `steps/climate_functions/local_accf.py` (input-column units; the PV + SI-versus-PVU note; NOx as NO2), `climaccf.py`, `open_airclim.py` (inventory units, NOx + species, `t_0`), `climate_functions/views.py`, `steps/climate_metrics/gwp.py` (output + units; the Method D `contrails_ef_J` note), `climate_metrics/report.py`, + `climate_metrics/views.py`, `steps/emissions/views.py` (per-segment masses, EIs, totals + basis). +- Runners and weather: `runners/fleet.py` (`FleetRunnerParams` input units, the + `runner.results` shape with units), `runners/flight.py` (the `source` DataFrame + contract), `runners/performance_runner.py` (output units; corrected the timestamp + example), `steps/weather/weather_factory.py` (corrected the "automatic unit conversions" + claim; the `rhi`, `sdr`, `pv` and time-window notes). + +The two upstream conversions are documented with the requested template at the JSON +boundary (`neats_io.py`), at the `q_fuel` definition (`NEATSFuel`) and where TAS enters +the performance model (`bada_model.py`). The runner entry points reference them in their +docstrings. + +**Checks performed.** + +- **AST equivalence.** For every modified `.py` file, the AST of the working copy with + docstrings removed was compared with that of `HEAD`: 27 of 27 files are identical. The + checker was first tested on purpose-made edits: it passes a docstring or comment change + and flags a one-token code change. This is the evidence that computational behaviour is + unchanged. +- **`uv run ruff check src/ tests/`**: all checks passed (baseline: passed). +- **`uv run ruff format --check src/ tests/`**: 84 files already formatted (baseline: same). +- **`uv run mypy src/pyneats/`**: no issues in 65 source files (baseline: same). +- **`make test-unit`** (tests without BADA or weather data): 150 passed, 1 skipped, + 4 xfailed, 103 deselected. The deselected tests need BADA or weather data, which is not + available here, so they were not run. +- **`git diff`** was read in full. Repeated template blocks were trimmed, and one unrelated + docstring fix was reverted to keep the patch to units. + +**Regulatory research (Section 6.0).** One researcher per question, each followed by two +independent verifiers. One verifier checked source fidelity by re-opening every cited source; +the other checked applicability and looked for contradicting sources. No finding was refuted. +Their corrections (citation pages, a missed guidance table, a superseded SAD draft, the +OpenAirClim v0.19.0 change) are incorporated above. After the research, the code comments +for `ff`, `ee`, `load_factor`, `hydrogen_per_carbon_ratio`, `takeoff_mass`, the +aromatics/sulphur/naphthalene basis and the EOBT label were updated to cite MRR Annex IIIa, +and the checks above were re-run. Not searched: PyNeats GitHub issues and PRs (`gh` is not +installed here), and the NEATS User Group SharePoint (login required). + +**Not done.** No pipeline, BADA, CoCiP or OpenAirClim calculation was run. No new tests were +added. Numerical equivalence is established by the AST comparison (no code changed). It is not +established by running tests. diff --git a/src/pyneats/core/compute_parameters.py b/src/pyneats/core/compute_parameters.py index 2b173e0..93bc68f 100644 --- a/src/pyneats/core/compute_parameters.py +++ b/src/pyneats/core/compute_parameters.py @@ -13,11 +13,15 @@ DEFAULT_JOBLIB_PREFERENCE : str Default joblib backend preference. DEFAULT_WEATHER_OFFSET_H : int - Default weather offset in hours. + Default weather time-window length in hours: ERA5 and the DWD Zarr path + load ``[asofdate, asofdate + offset]``. DEFAULT_HORIZONTAL_RES_DEG : float - Default horizontal resolution in degrees. + Default horizontal grid resolution in degrees (ERA5 ``grid``). DEFAULT_SDR_ACCUMULATE_DT_S : int - Default SDR accumulation time in seconds. + Default SDR accumulation step in seconds (3600 s = 1 h). When the DWD Zarr + cache is built, the ``sdr`` values (assumed W m-2) are multiplied by it and + labelled J m**-2; the live NetCDF path does not apply it (see + docs/internal/pyneats_units_audit.md). DEFAULT_MET_CHUNKS : Mapping[str, Any] Default chunk sizes for meteorological data. DEFAULT_RAD_CHUNKS : Mapping[str, Any] @@ -40,13 +44,14 @@ DEFAULT_BATCH_SIZE: Final[int] = 16 DEFAULT_JOBLIB_PREFERENCE: Final[str] = "processes" -DEFAULT_WEATHER_OFFSET_H: Final[int] = 36 +DEFAULT_WEATHER_OFFSET_H: Final[int] = 36 # h -DEFAULT_HORIZONTAL_RES_DEG: Final[float] = 0.25 -DEFAULT_SDR_ACCUMULATE_DT_S: Final[int] = 3600 # 1 hour +DEFAULT_HORIZONTAL_RES_DEG: Final[float] = 0.25 # degrees +DEFAULT_SDR_ACCUMULATE_DT_S: Final[int] = 3600 # s (1 hour) # Default chunk dimensions for xarray operations using zarr backend +# (number of elements per chunk along each dimension) DEFAULT_MET_CHUNKS: Final[Mapping[str, Any]] = { "time": 1, diff --git a/src/pyneats/core/fleet_utils.py b/src/pyneats/core/fleet_utils.py index e5597e5..12cf5f7 100644 --- a/src/pyneats/core/fleet_utils.py +++ b/src/pyneats/core/fleet_utils.py @@ -2,7 +2,7 @@ Utilities for Fleet ↔ List[Flight] conversions. This module provides utility functions to convert between lists of Flight objects and Fleet objects, -preserving and restoring fuel information and original columns. +carrying fuel information across the conversion and restoring original columns. """ from __future__ import annotations @@ -34,8 +34,10 @@ def flights_to_fleet(flights: list[Flight]) -> Fleet: """ Convert a list of Flight objects to a Fleet, preserving fuel information. - Fuel properties (q_fuel, ei_h2o) are extracted from flight.fuel and stored - as columns in the Fleet dataframe. The original fuel object is set to None. + Fuel properties are copied from flight.fuel into per-waypoint columns + ``q_fuel`` [J/kg, lower calorific value] and ``ei_h2o`` [kg H2O / kg fuel], + which pycontrails SAC/CoCiP read on the Fleet. The original fuel object is + set to None. Flights with different columns are harmonized by padding missing columns with NaN values. Original column sets are stored in attrs to enable restoration during fleet_to_flights(). @@ -101,9 +103,11 @@ def fleet_to_flights( """ Convert a Fleet object back to a list of Flight objects, restoring fuel information. - Fuel properties (q_fuel, ei_h2o) are extracted from Fleet columns and used - to reconstruct flight.fuel objects. Columns that were added during Fleet - conversion (NaN padding) are removed to restore original column sets. + flight.fuel is rebuilt with ``NEATSFuel.from_attrs(flight.attrs)``, i.e. + from the flight attributes, not from the Fleet columns. Columns that were + added during Fleet conversion (NaN padding and the ``q_fuel``/``ei_h2o`` + columns) are removed without being read, restoring original column sets + (see docs/internal/pyneats_units_audit.md). Any flights present in the original manifest but missing from the Fleet are returned as error records. @@ -111,7 +115,8 @@ def fleet_to_flights( Parameters ---------- fleet : Fleet - Fleet object with fuel information as columns. + Fleet object produced by flights_to_fleet(). Its fuel columns are not + read; fuel is rebuilt from each flight's attrs. step_name : str Name of the processing step, used in error messages for dropped flights. @@ -142,7 +147,7 @@ def fleet_to_flights( for col in fleet_columns - original_columns: flight.data.pop(col, None) - # Restore fuel object + # Rebuild the fuel object from attrs; the q_fuel/ei_h2o columns are not read flight.fuel = NEATSFuel.from_attrs(flight.attrs) # Create error records for dropped flights diff --git a/src/pyneats/core/neats_default_parameters.py b/src/pyneats/core/neats_default_parameters.py index 5a1882e..e5051fe 100644 --- a/src/pyneats/core/neats_default_parameters.py +++ b/src/pyneats/core/neats_default_parameters.py @@ -24,27 +24,27 @@ DEFAULT_CLIMATE_IMPACT : str Default climate impact metric. DEFAULT_INTERPOLATION_TIME : str - Default interpolation time step. + Default interpolation time step, as a pandas frequency string ("1min" = 60 s). DEFAULT_BADA4_VERSION : str Default BADA4 version. DEFAULT_BADA3_VERSION : str Default BADA3 version. DEFAULT_PAYLOAD_FACTOR : float - Default payload factor. + Default payload factor: dimensionless fraction of maximum payload mass. DEFAULT_FUEL_RESERVE_FRACTION : float - Default fuel reserve fraction. + Default fuel reserve fraction: dimensionless fraction of trip fuel mass. DEFAULT_MAX_MASS_ESTIMATION_ITER : int Default max mass estimation iterations. DEFAULT_MAX_REL_MASS_DIFF : float - Default max relative mass difference. + Default max relative mass difference (dimensionless fraction). DEFAULT_TRUE_AIR_SPEED_SMOOTHING_WINDOW : int - Default smoothing window for true air speed. + Default smoothing window for true air speed, in number of points. DEFAULT_DELTA_TAU_COMPUTE_METHOD : Literal Default delta tau compute method. DEFAULT_DELTA_TAU_FILL_METHOD : Literal Default delta tau fill method. DEFAULT_ROCD_PHASE_THRESHOLD : float - Default rate of climb/descent phase threshold. + Default rate of climb/descent phase threshold, in ft/min. DEFAULT_EMISSIONS_KWARGS : Mapping[str, Any] Default emissions model keyword arguments. DEFAULT_COCIP_KWARGS : Mapping[str, Any] @@ -52,23 +52,23 @@ DEFAULT_CLIMACCF_KWARGS : Mapping[str, Any] Default ClimAccf model keyword arguments. DEFAULT_AROMATICS_CONTENT : float - Default aromatics content for fuel. + Default aromatics content for fuel (fraction scale; unused). DEFAULT_SULPHUR_CONTENT : float - Default sulphur content for fuel. + Default sulphur content for fuel (fraction scale; unused). DEFAULT_NAPHTHALEN_CONTENT : float - Default naphthalene content for fuel. + Default naphthalene content for fuel (fraction scale; unused). DEFAULT_HYDROGEN_CONTENT : float - Default hydrogen content for fuel. + Default hydrogen content for fuel, in wt.% (percent by mass). REFERENCE_Q_FUEL : float - Reference q_fuel value. + Reference q_fuel value, in J/kg (pycontrails JetA lower calorific value). DEFAULT_Q_FUEL : float - Default q_fuel value. + Default q_fuel value, in J/kg (lower calorific value). OPEN_AIRCLIM_LEVELS : tuple of int - Default pressure levels for Open Air Clim. + Default pressure levels for Open Air Clim, in hPa. OPEN_AIRCLIM_GRID_RES : float - Default grid resolution for Open Air Clim. + Default grid resolution for Open Air Clim, in degrees. OPEN_AIRCLIM_COMPUTATION_HORIZON : int - Default computation horizon for Open Air Clim. + Default computation horizon for Open Air Clim, in years. """ from collections.abc import Mapping @@ -121,6 +121,7 @@ # Interpolation defaults # Specifies the temporal resolution for trajectory +# pandas frequency string passed to Flight.resample_and_fill ("1min" = 60 s) DEFAULT_INTERPOLATION_TIME: Final[str] = "1min" # Performance model defaults @@ -134,10 +135,13 @@ # Very conservative threshold on the fuel flow. An error in BADA will be thown if higher # Further studies need to be performed to better handle out of envelop data points +# Unit: kg/s, whole aircraft; compared only with the BADA-computed fuel flow +# (after the q_fuel correction). DEFAULT_FF_OUTLIER_THRESHOLD = 100 # Mass estimation parameters # Used in iterative aircraft mass computation from fuel consumption +# Fraction of the maximum payload mass MPL (MTOW - OEW when MPL is unavailable) DEFAULT_PAYLOAD_FACTOR: Final[float] = ( 1.0 # Multiplier for expected payload (1.0 = full capacity) ) @@ -176,21 +180,24 @@ ) # Fuel burn guardrail: reject if total fuel > (MTOW - OEW) * threshold -# A ratio of 1.1 allows 10% margin above the useful payload capacity +# (masses in kg). A ratio of 1.2 allows a 20% margin above MTOW - OEW. DEFAULT_FUEL_BURN_THRESHOLD: Final[float] = 1.2 # Altitude filter for performance fallback # Minimum Flight Level (FL) threshold: points below this FL are filtered # when retrying after a performance step failure. -# FL15 = 2000 ft ≈ 610 m — filters out ground-level / taxi data. +# FL20 = 2000 ft = 609.6 m — filters out ground-level / taxi data. +# Converted to m (FL x 100 ft) and compared with the altitude column [m]. DEFAULT_MIN_ALTITUDE_FL: Final[float] = 20 -# Maximum ratio of filtered points before giving up on altitude fallback. +# Maximum ratio of filtered points before giving up on altitude fallback +# (dimensionless fraction [0, 1] of trajectory points). # If more than 80% of trajectory points are below min FL, the flight is # likely a low-altitude operation — retrying is unlikely to help. DEFAULT_MAX_ALTITUDE_FILTER_RATIO: Final[float] = 0.8 # Guardrail: reject non-CO2 species producing > MAX_NONCO2_CO2_RATIO × CO2 baseline +# (dimensionless ratio of CO2-equivalent mass to emitted CO2 mass, both in kg) MAX_NONCO2_CO2_RATIO: Final[float] = 1000.0 # Emissions model parameters @@ -262,6 +269,7 @@ "horizontal_resolution": None, # Use native data resolution "unit_K_per_kg_fuel": False, # Output in ATR units (not K per kg fuel) } +# Compared with the pycontrails air_pressure column [Pa]; 40000 Pa = 400 hPa DEFAULT_ACCF_VALIDITY_PRESSURE = 40000 # Pressure level (Pa) for ACCF validity check # Weather field interpolation configuration @@ -297,6 +305,9 @@ # Jet fuel composition and thermodynamic properties # Used in emissions and energy calculations. Default values from consortium (DLR/TO70) +# Aromatics, sulphur and naphthalene defaults are fraction scale (0.25 = 25%); the AO +# schema does not define the scale of the corresponding fields. None of the three is +# used in calculations. Hydrogen content is in wt.% (percent by mass), not a fraction. DEFAULT_AROMATICS_CONTENT: Final[float] = ( 0.25 # 25% - Default value (not used in current calculations) ) @@ -307,14 +318,19 @@ 0.03 # 3% - Default value (not used in current calculations) ) DEFAULT_HYDROGEN_CONTENT: Final[float] = ( - 13.79 # 13.79% - Default value provided by consortium (DLR/TO70) for fuel composition + 13.79 # wt.% - Default value from consortium (DLR/TO70) for fuel composition ) -# Fuel energy content (specific energy) +# Fuel energy content (specific energy), lower calorific value in J/kg of fuel # Critical for converting fuel consumption to energy and emissions estimates +# REFERENCE_Q_FUEL is the pycontrails JetA q_fuel (pycontrails/core/fuel.py); it is +# not read from BADA files. The BADA fuel flow is divided by q_fuel / REFERENCE_Q_FUEL. +# Because NEATSFuel always carries a q_fuel (DEFAULT_Q_FUEL when none is supplied), +# this correction always applies (see docs/internal/pyneats_units_audit.md). REFERENCE_Q_FUEL: Final[float] = ( - 43_130_000.0 # J/kg - Baseline value from PyContrails/BADA + 43_130_000.0 # J/kg - Baseline value from PyContrails JetA ) +# DEFAULT_Q_FUEL is used when neither the AO nor the airport fuel file supplies q_fuel DEFAULT_Q_FUEL: Final[float] = ( 42_800_000.0 # J/kg - Default value provided by consortium (DLR/TO70) ) diff --git a/src/pyneats/core/neats_fuel.py b/src/pyneats/core/neats_fuel.py index a99560b..247a22c 100644 --- a/src/pyneats/core/neats_fuel.py +++ b/src/pyneats/core/neats_fuel.py @@ -34,23 +34,51 @@ class NEATSFuel(SAFBlend): Parameters ---------- hydrogen_content : float, optional - Hydrogen content in weight percent (wt.% H). + Hydrogen content in wt.% (percent by mass, e.g. 13.8 for Jet-A), the + same scale as the AO field "percentage of hydrogen mass content". Not + a fraction. Defaults to ``DEFAULT_HYDROGEN_CONTENT`` (13.79 wt.%). h_c_ratio : float, optional - Atomic H/C ratio. + H/C ratio, dimensionless, atomic (molar): hydrogen atoms per carbon + atom (Reg. (EU) 2018/2066 Annex IIIa §1(36)). Converted to wt.% with + ``r * 1.008 / (12.011 + r * 1.008) * 100``; used only when + ``hydrogen_content`` is None. q_fuel : float, optional - Lower heating value (J/kg). + PyNEATS input unit: J/kg of fuel (lower calorific value). + The AO JSON field ``calorific_value`` is in MJ/kg of fuel. + The NEATS application multiplies by 1e6 before calling PyNEATS. + Values received here must already be in J/kg; do not convert again. + Defaults to ``DEFAULT_Q_FUEL`` (42.8e6 J/kg). pct_blend_gate : float, optional Blend gate for SAFBlend compatibility. sulphur_content : float, optional - Sulphur content (not used in calculations). + Sulphur content, % by mass basis (Reg. (EU) 2018/2066 Annex IIIa §5(4) default); + fraction vs percentage scale not defined by the AO schema. Stored + only, not used in calculations. aromatics_content : float, optional - Aromatics content (not used in calculations). + Aromatics content, % by volume basis (Annex IIIa §5(4) default); + fraction vs percentage scale not defined by the AO schema. Stored + only, not used in calculations. naphthalene : float, optional - Naphthalene content (not used in calculations). + Naphthalene content, % by volume basis (Annex IIIa §5(4) default); + fraction vs percentage scale not defined by the AO schema. Stored + only, not used in calculations. name : str, optional Name of the fuel. + + Notes + ----- + Derived emission indices, all in kg per kg of fuel: + + - ``ei_co2`` is the pycontrails JetA value, independent of hydrogen content. + - ``ei_h2o`` is the JetA value scaled by the resolved hydrogen content + (wt.%) divided by 13.8 (JetA hydrogen content, wt.%). + - ``ei_so2`` (and ``ei_sulphates`` derived from it) is the JetA value + whatever ``sulphur_content`` is, so SO2 and sulphate outputs do not + depend on the supplied sulphur (see docs/internal/pyneats_units_audit.md). """ + # Class defaults are fraction scale (see neats_default_parameters); values + # supplied by the AO keep the AO scale, which the schema does not define. aromatics_content: float = ( DEFAULT_AROMATICS_CONTENT # default value - Not used in calculations ) @@ -82,13 +110,14 @@ def __init__( h_wt_pct = float(hydrogen_content) elif h_c_ratio is not None: r = float(h_c_ratio) - # Convert H/C atomic ratio to hydrogen *mass percent* - # (×100 to keep PyContrails' convention) + # Convert H/C atomic (molar) ratio to hydrogen *mass percent* using + # atomic masses H = 1.008, C = 12.011 g/mol + # (×100 to keep PyContrails' wt.% convention) h_wt_pct = (r * 1.008) / (12.011 + r * 1.008) * 100.0 else: h_wt_pct = DEFAULT_HYDROGEN_CONTENT - # --- Resolve q_fuel + # --- Resolve q_fuel [J/kg, LHV]; already in J/kg on entry (see docstring) qf = float(q_fuel) if q_fuel is not None else DEFAULT_Q_FUEL # --- Storing Optional Fuel Attributes from AOs, not used at this stage @@ -101,7 +130,8 @@ def __init__( if naphthalene is not None: object.__setattr__(self, "naphthalene", naphthalene) - # --- Derive dependent indices consistently + # --- Derive dependent indices consistently [kg/kg fuel] + # ei_co2 and ei_so2 are the JetA values; only ei_h2o follows h_wt_pct ei_co2 = base.ei_co2 ei_h2o = base.ei_h2o * (h_wt_pct / base.hydrogen_content) ei_so2 = base.ei_so2 @@ -172,6 +202,14 @@ def from_attrs(cls, attrs: Mapping[str, Any]) -> NEATSFuel: Handles the case where fuel properties are stored as arrays (from column serialization) by extracting the first element. + Values are passed through unconverted, so they must already be in the + units of the constructor (``q_fuel`` in J/kg, ``hydrogen_content`` in + wt.%). Only the key ``naphthalene`` is read; the airport-file path stores + the value under ``naphtalene`` and the AO JSON path drops it, so a parsed + naphthalene value is not picked up. Alternative keys are joined with + ``or``, so an explicit 0.0 under the first key falls back to the second + key, then to the default (see docs/internal/pyneats_units_audit.md). + Parameters ---------- attrs : Mapping[str, Any] diff --git a/src/pyneats/core/physics.py b/src/pyneats/core/physics.py index 927b147..4d7d863 100644 --- a/src/pyneats/core/physics.py +++ b/src/pyneats/core/physics.py @@ -7,29 +7,36 @@ Attributes ---------- METRICS_HORIZONS : tuple of int - Supported time horizons for metrics. + Supported time horizons for metrics, in years. SURFACE_EARTH : float Surface area of the Earth in m². SECONDS_PER_YEAR : int - Number of seconds in a year. + Number of seconds in a year: 31,556,952 s = 365.2425 d (Gregorian mean + year). docs/explanation/mrv-specification.md uses 31,536,000 s (365 d) + (see docs/internal/pyneats_units_audit.md). SOLAR_CONSTANT : float Solar constant in W m^-2. CO2_AGWP_COEFF_WM2YR_PER_KG : Mapping[int, float] - CO2 GWP coefficients in W·m⁻²·yr·kg⁻¹. + CO2 AGWP coefficients C(H) in W m-2 yr per kg CO2, keyed by horizon in + years. CONVERSION_FACTORS_AGWP_TO_RF : dict - Conversion factors from AGWP to RF for various species and horizons. + Conversion factors from RF to AGWP(H) (Dahlmann et al., 2025; the name + reads the other way), keyed by horizon in years and species. Absolute + units not verified. CONVERSION_FACTORS_ATR_TO_RF : dict - Conversion factors from ATR to RF for various species and horizons. + Conversion factors from RF to ATR(H) (Dahlmann et al., 2025; the name + reads the other way), keyed by horizon in years and species. Absolute + units not verified. EFFICACY : dict - Efficacy values for different species. + Efficacy values for different species (dimensionless). RF_BACKWARD_FACTOR : dict - Backward calculation factors for radiative forcing. + Backward calculation factors for radiative forcing (divisors). ACCF_SCALE_03 : Mapping[str, float] - Scaling factors for O3 ACCFs. + Scaling factors for O3 ACCFs (divisors), keyed by aCCF version. ACCF_SCALE_CH4 : Mapping[str, float] - Scaling factors for CH4 ACCFs. + Scaling factors for CH4 ACCFs (divisors), keyed by aCCF version. ACCF_SCALE_H2O : Mapping[str, float] - Scaling factors for H2O ACCFs. + Scaling factors for H2O ACCFs (divisors), keyed by aCCF version. DEFAULT_PRESSURE_LEVELS_HPA : tuple of float Default pressure levels in hPa. """ @@ -55,11 +62,15 @@ METRICS_HORIZONS: Final[tuple[int, ...]] = (20, 50, 100) SURFACE_EARTH: Final[float] = 5.101e14 # m² +# Gregorian mean year (365.2425 d); used as s/yr to turn contrail energy forcing +# [J] into EAGWP [W m-2 yr] SECONDS_PER_YEAR: Final[int] = 31_556_952 # s SOLAR_CONSTANT: Final[float] = 1360.0 # [W m^-2] # CO2 GWP coefficients from AirClim (Dahlmann et al., 2025) -# C(H) in W·m⁻²·yr·kg⁻¹ (per kgCO2) +# C(H) in W·m⁻²·yr·kg⁻¹ (per kgCO2); keys are horizons in years. +# Values from commit bf3aa7b; docs/explanation/mrv-specification.md still lists +# the older 24.16 / 47.57 / 74.36e-15 (see docs/internal/pyneats_units_audit.md). CO2_AGWP_COEFF_WM2YR_PER_KG: Mapping[int, float] = { 20: 24.12e-15, 50: 47.50e-15, @@ -68,6 +79,9 @@ # New generic conversion factors to convert RF to AGWP or ATR # (Dahlmann et al., 2025) for Pulse 2025 scenario +# Absolute units are not stated here and not verified. gwp.py uses only the ratio +# k_AGWP(H) / k_ATR(H0), combined with efficacy**2 / RF_BACKWARD_FACTOR, to map +# ATR [K] to EAGWP [W m-2 yr], so that combined scale must carry W m-2 yr per K. CONVERSION_FACTORS_AGWP_TO_RF: Final[dict[int, dict[str, float]]] = { 20: {"CH4": 10.6318, "O3": 1.0192, "H2O": 1.0192, "PMO": 10.6318}, 50: {"CH4": 13.0968, "O3": 1.0192, "H2O": 1.0192, "PMO": 13.0968}, @@ -95,17 +109,21 @@ # Dahlmann et al., 2025 # Therefore used here to discount accf output value before applying # conversion factors from Dahlmann et al., 2025 +# Used as a divisor in gwp.py. The values are approximately 1 / ACCF_SCALE_*["V1.0"] +# (1/1.97 = 0.508, 1/2.03 = 0.493, 1/1.92 = 0.521), so the two tables are coupled. RF_BACKWARD_FACTOR: Final[dict[str, float]] = { "CH4": 0.492, "O3": 0.508, "H2O": 0.52, } -# ACCFs scaling factors +# ACCFs scaling factors, keyed by aCCF version; local_accf divides the aCCF output +# by them ACCF_SCALE_03: Final[Mapping[str, float]] = {"V1.0": 1.97, "V1.0A": 11.0} ACCF_SCALE_CH4: Final[Mapping[str, float]] = {"V1.0": 2.03, "V1.0A": 35.0} ACCF_SCALE_H2O: Final[Mapping[str, float]] = {"V1.0": 1.92, "V1.0A": 3.0} +# Pressure levels in hPa DEFAULT_PRESSURE_LEVELS_HPA: Final[tuple[float, ...]] = ( 550, 500, diff --git a/src/pyneats/core/views.py b/src/pyneats/core/views.py index 6634033..3a544f0 100644 --- a/src/pyneats/core/views.py +++ b/src/pyneats/core/views.py @@ -218,6 +218,21 @@ def matches(cls, flight: Flight) -> bool: # Overload method def to_dict(self) -> dict[str, Any]: + """ + Convert the FlightView to a dictionary of JSON-compatible values. + + Data columns and attrs are merged; data keys overwrite duplicate attrs. + + Returns + ------- + dict[str, Any] + ``time`` is serialised as integer seconds since the Unix epoch, + truncated (``astype("datetime64[s]")``), from the tz-naive UTC times + that pycontrails stores. ``timedelta64[ns]`` columns (e.g. + ``contrail_age``) are serialised by the pycontrails ``NumpyEncoder`` + as integer nanoseconds. Other values keep the units of the step that + produced them. + """ np_encoder = json_utils.NumpyEncoder() def encode(key: str, obj: Any) -> Any: @@ -227,11 +242,12 @@ def encode(key: str, obj: Any) -> Any: # Convert numpy objects to python objects if isinstance(obj, np.ndarray | np.generic): - # round time to unix seconds + # truncate time to integer unix seconds (no rounding) if key == "time": return np_encoder.default(obj.astype("datetime64[s]").astype(int)) - # round specific keys in precision + # other arrays: converted to lists, values unchanged + # (timedelta64[ns] becomes integer nanoseconds) return np_encoder.default(obj) # Pass through everything else diff --git a/src/pyneats/runners/fleet.py b/src/pyneats/runners/fleet.py index c64ee39..66e5845 100644 --- a/src/pyneats/runners/fleet.py +++ b/src/pyneats/runners/fleet.py @@ -143,7 +143,29 @@ class StepFailure: # - kw_only=True avoids dataclass field ordering issues when subclassing. @dataclass(kw_only=True) class FleetRunnerParams(RunnerConfig): - """Fleet runner configuration, extending the FlightRunner-style RunnerConfig.""" + """Fleet runner configuration, extending the FlightRunner-style RunnerConfig. + + Attributes: + zarr_paths: Weather Zarr stores (met, rad, optional wind). Every flight + in the run is computed with these stores; the AO JSON field + ``nwp_date_time`` is not read by PyNEATS. + trajectory_json_filepath: JSON file holding a list (array) of AO-format + flight objects, already converted by the NEATS application (``tas`` + in m/s, ``calorific_value`` in J/kg). A bare single object is not + accepted (see ``neats_json_to_flights``). + trajectory_dataframe: Concatenated DataFrame of all flights. Trajectory + columns use the Flight4D names, in the units the default ``neats`` + parser expects: ``latitude`` and ``longitude`` (degrees), + ``altitude`` (flight level, hundreds of ft), ``time`` (datetime64, + naive = UTC, or a string in the parser ``date_format``), + ``fuel_flow`` (kg/s, whole aircraft), ``aircraft_mass`` (kg), + ``engine_efficiency`` (dimensionless fraction [0, 1]) and + ``true_airspeed`` (m/s). All other columns + (``flight_id``, ``departure_airport``, ``arrival_airport``, + ``aobt``, ``aircraft_type``, ..., ``q_fuel`` in J/kg, + ``hydrogen_content`` in wt.% (percent by mass)) are flight-level + and are moved to ``attrs`` by ``split_df_into_flights``. + """ # Required fleet resource zarr_paths: ZarrPaths @@ -298,7 +320,12 @@ def _create_error_record( def _build_zero_climate_result(error_record: dict[str, Any]) -> dict[str, Any]: - """Build a climate result with all-zero metrics for a low-altitude flight.""" + """Build a climate result with all-zero metrics for a low-altitude flight. + + ``co2_baseline_kg`` and ``fuel_burn_kg`` are NaN (not computed), + ``contrails_ef_J`` is 0.0 and every ``EAGWP_Wm2yr`` (W m-2 yr) and + ``CO2eq_kg`` (kg CO2e) entry is 0.0. + """ fi = dict(error_record.get("flight_information", {})) fi["low_altitude_flight"] = True fi["co2_baseline_kg"] = float("nan") @@ -367,7 +394,12 @@ def _process_parallel_results( def _merge_engine_group( group: list[FlightWithEmissions], parent_key: str ) -> FlightWithEmissions: - """Average numeric DataFrame columns across engine copies; restore original flight identity.""" + """Average numeric columns over engine copies; restore original flight identity. + + Each copy is the same trajectory run with one engine type, so the merged + columns and numeric attrs (e.g. ``total_co2`` in kg) are the mean over + engine types, not a sum. + """ import numpy as np # Use first flight as the structural template @@ -425,6 +457,47 @@ class FleetRunner(Runner): This runner keeps "continue-on-error" semantics: per-flight failures are recorded and filtered out rather than aborting the whole run. + + After ``eval()``, ``self.results`` has the shape:: + + { + "fleet_meta_data": {...}, + "flight_results": [ + { + "flight_information": { + ..., # FlightReport fields (flight_id, aobt, ...) + "co2_baseline_kg": float, # flight total CO2 [kg] + "fuel_burn_kg": float, # flight total fuel [kg] + "contrails_ef_J": float, # contrail energy forcing [J] + }, + "climate_metrics": [ + { + "species": "CO2", # CO2 first; also CH4, O3, H2O, Contrails + "value": [ + { + "horizon": 20, # years (20, 50, 100) + "EAGWP_Wm2yr": float, # [W m-2 yr] + "CO2eq_kg": float, # [kg CO2e] + }, + ..., + ], + }, + ..., + ], + }, + ..., + ], + } + + Successful flights come first, then low-altitude zero results, then + error records. Error records are ``{"flight_information": {...}, + "error": str}``. Low-altitude zero results have NaN ``co2_baseline_kg`` + and ``fuel_burn_kg`` and 0.0 metrics. When metrics come from OpenAirClim + AGWP attrs (Method D in ``gwp.py``), ``contrails_ef_J`` holds the 20-year + contrail AGWP in W m-2 yr, not an energy forcing in J (see + docs/internal/pyneats_units_audit.md). Flights with a list ``engine_uid`` + are run once per engine and merged back into one result (see + ``_merge_engine_group``). """ parser: Callable[[pd.DataFrame], Flight4D] @@ -558,6 +631,8 @@ def _load_data(self) -> Self: ) def _load_trajectories(self) -> Self: + # No unit conversion here: true_airspeed must already be m/s and q_fuel + # J/kg (the NEATS application converts the AO knots and MJ/kg values). if self.cfg.trajectory_json_filepath is not None: path = Path(self.cfg.trajectory_json_filepath) logger.info("Loading trajectories from %s", path) @@ -806,8 +881,8 @@ def results_as_dataframe(self) -> pd.DataFrame: """ Flatten self.results into one row per flight. Climate metrics are expanded into wide columns: - __EAGWP_Wm2yr - __CO2eq_kg + __EAGWP_Wm2yr (W m-2 yr; horizon in years) + __CO2eq_kg (kg CO2e) """ rows: list[dict[str, Any]] = [] diff --git a/src/pyneats/runners/flight.py b/src/pyneats/runners/flight.py index 25efcb9..eaf1dab 100644 --- a/src/pyneats/runners/flight.py +++ b/src/pyneats/runners/flight.py @@ -95,6 +95,26 @@ def __init__( bada_path: str | None = None, airport_fuel_path: str | None = None, ) -> None: + """Build the per-flight pipeline steps from the registry. + + Args: + weather: Weather provider used for the intersection step. + source: One flight as a DataFrame. Trajectory columns use the + Flight4D names, in the units the default ``neats`` parser + expects: ``latitude`` and ``longitude`` (degrees), + ``altitude`` (flight level, hundreds of ft), ``time`` + (datetime64, naive = UTC, or a string in the parser + ``date_format``), ``fuel_flow`` (kg/s, whole aircraft), + ``aircraft_mass`` (kg), ``engine_efficiency`` (dimensionless + fraction [0, 1]) and ``true_airspeed`` (m/s). Flight-level + fields (``flight_id``, ``aobt``, ``aircraft_type``, ..., + ``q_fuel`` in J/kg, ``hydrogen_content`` in wt.% (percent by + mass)) go in ``source.attrs``. A list ``engine_uid`` is + rejected by the NEATS parser; use FleetRunner for those. + cfg: Step names and per-step params. + bada_path: Root of the BADA 3 / BADA 4 coefficient files. + airport_fuel_path: Optional airport -> fuel properties file. + """ self.cfg = cfg or RunnerConfig() self._source = None # initialize backing field before using the property diff --git a/src/pyneats/runners/performance_runner.py b/src/pyneats/runners/performance_runner.py index e6a3938..a684982 100644 --- a/src/pyneats/runners/performance_runner.py +++ b/src/pyneats/runners/performance_runner.py @@ -7,8 +7,9 @@ Classes ------- FleetRunnerPerformanceOnly - Stops after _performance() and serialises fuel_flow / aircraft_mass / - true_airspeed trajectory vectors as JSON. + Stops after _performance() and serialises fuel_flow (kg/s, whole + aircraft) / aircraft_mass (kg) / true_airspeed (m/s) / engine_efficiency + (dimensionless fraction) trajectory vectors as JSON. """ from __future__ import annotations @@ -28,7 +29,9 @@ __all__ = ["FleetRunnerPerformanceOnly"] -# Columns to extract from each FlightWithPerformance trajectory +# Columns to extract from each FlightWithPerformance trajectory: +# fuel_flow kg/s (whole aircraft), aircraft_mass kg, true_airspeed m/s, +# engine_efficiency dimensionless fraction. Emitted without conversion. _PERF_VECTOR_COLS = ("fuel_flow", "aircraft_mass", "true_airspeed", "engine_efficiency") # Flight-level attrs included in the output @@ -43,7 +46,7 @@ def _to_list_or_null(series: Any) -> list[float | None]: - """Convert a pandas Series (or array-like) to a JSON-safe list.""" + """Convert a pandas Series (or array-like) to a JSON-safe list (NaN -> None).""" result: list[float | None] = [] for v in series: try: @@ -77,7 +80,7 @@ class FleetRunnerPerformanceOnly(FleetRunnerLargeEmitter): "engine_uid": "...", # only if present }, "performance": { - "timestamps": ["2025-01-01T10:00:00+00:00", ...], + "timestamps": ["2025-01-01T10:00:00", ...], "fuel_flow": [1.2, 1.3, null, ...], "aircraft_mass": [70000.0, 69950.0, ...], "true_airspeed": [240.0, 245.0, ...], @@ -88,6 +91,16 @@ class FleetRunnerPerformanceOnly(FleetRunnerLargeEmitter): ], } + Units: ``timestamps`` are ISO 8601 strings of tz-naive UTC datetimes, so + no offset is printed; ``fuel_flow`` kg/s (whole aircraft); + ``aircraft_mass`` kg; ``true_airspeed`` m/s; ``engine_efficiency`` + dimensionless fraction. NaN values are written as ``null``. + + Flights with a list ``engine_uid`` are not merged here: the engine-merge + step runs inside ``_emissions``, which this runner skips, so each engine + copy gives its own record (``flight_id`` ``"__eng"``, scalar + ``engine_uid``). + Error records (flights that failed during parsing / weather / performance) are included in ``flight_results`` with an ``"error"`` key, exactly as in the full pipeline. diff --git a/src/pyneats/steps/climate_functions/climaccf.py b/src/pyneats/steps/climate_functions/climaccf.py index c693bbd..8f29502 100644 --- a/src/pyneats/steps/climate_functions/climaccf.py +++ b/src/pyneats/steps/climate_functions/climaccf.py @@ -66,6 +66,20 @@ def make_accf_surface_view(surface: MetDataset) -> MetDataset: - variable aliasing / renaming for ClimAccf expectations - minimal attrs/unit tweaks + Notes + ----- + Only names and ``units`` attributes change; no values are converted + (see docs/internal/pyneats_units_audit.md): + + - ``toa_outgoing_longwave_flux`` (pycontrails OLR, positive upward, + W m-2) is renamed to ECMWF ``top_net_thermal_radiation`` without a + sign change, although ECMWF defines ``ttr`` as minus OLR. + - ``surface_solar_downward_radiation`` is relabelled ``W m**-2``. + ``DWDFactory.build_cache`` (default ``sdr_accumulate_dt_s=3600``) + stores this field as W m-2 x 3600 labelled ``J m**-2``; after the + relabel, pycontrails ``ACCF`` multiplies by 3600 again. The DWD + source field (``sdr``) is TOA incoming solar according to the DWD + store history, not surface radiation. """ # Work on an xr.Dataset view @@ -88,6 +102,8 @@ def make_accf_surface_view(surface: MetDataset) -> MetDataset: # 2) Align units on the surface shortwave flux if needed if SurfaceSolarDownwardRadiation.standard_name in ds: # Ensure units are consistent with the radiation flux variables expected by ClimAccf (e.g., W m-2) + # Relabels the units attribute only; values are not converted + # (see docs/internal/pyneats_units_audit.md). ds[SurfaceSolarDownwardRadiation.standard_name].attrs.update( {"units": TOAOutgoingLongwaveFlux.units} ) @@ -152,14 +168,25 @@ def run(self, flight: FlightWithEmissions) -> FlightWithNonCO2Impact: Delegates to the ``climaccf`` package to obtain per-segment aCCF outputs (K per kg NOx or K per kg fuel depending on ``accf_kwargs``), then converts them to ``ATR_20_{CH4,O3,H2O}`` columns on the flight. + ``climaccf`` is not a declared PyNEATS dependency; this backend is not + the default. + + With ``accf_kwargs["unit_K_per_kg_fuel"] = False`` (default), the + CH4/O3 aCCFs are multiplied here by the flight's ``nox_ei`` + [kg NOx (as NO2) / kg fuel]. With ``True``, the flight's ``nox_ei`` is + not used: with ``merged`` True (pycontrails default) CLIMaCCF applies + its own NOx emission index, selected by the pycontrails ``ACCF`` + ``nox_ei`` parameter (default ``"TTV"``, a fixed 0.013 kg/kg) + (see docs/internal/pyneats_units_audit.md). NaN ``fuel_burn`` and + ``nox_ei`` are treated as 0. Args: - flight: Flight with emissions data. Must expose ``fuel_burn`` and - ``nox_ei`` columns. + flight: Flight with emissions data. Must expose ``fuel_burn`` + [kg fuel per segment] and ``nox_ei`` [kg/kg fuel] columns. Returns: ``FlightWithNonCO2Impact`` — flight with ``ATR_20_CH4``, - ``ATR_20_O3``, ``ATR_20_H2O`` columns (Kelvin) consumed by + ``ATR_20_O3``, ``ATR_20_H2O`` columns (K per segment) consumed by ``GWPMetrics`` downstream. Raises: @@ -172,15 +199,15 @@ def run(self, flight: FlightWithEmissions) -> FlightWithNonCO2Impact: nox_ei = pd.to_numeric(df["nox_ei"], errors="coerce").fillna(0.0) if not self.params.accf_kwargs["unit_K_per_kg_fuel"]: - # Convert aCCF outputs from K per kg NOx to K per kg H2O and K per Fuel + # Convert O3/CH4 aCCFs from K per kg NOx (as NO2) to K per kg fuel df["aCCF_CH4"] = df["aCCF_CH4"] * nox_ei df["aCCF_O3"] = df["aCCF_O3"] * nox_ei - # Compute final ART_20 (K) from aCCF outputs (K per kg fuel) + # ATR_20 [K per segment] = fuel_burn [kg] x aCCF [K per kg fuel] - flight["ATR_20_CH4"] = fuel_burn * df["aCCF_CH4"] # K - flight["ATR_20_O3"] = fuel_burn * df["aCCF_O3"] # K - flight["ATR_20_H2O"] = fuel_burn * df["aCCF_H2O"] # K + flight["ATR_20_CH4"] = fuel_burn * df["aCCF_CH4"] # K per segment + flight["ATR_20_O3"] = fuel_burn * df["aCCF_O3"] # K per segment + flight["ATR_20_H2O"] = fuel_burn * df["aCCF_H2O"] # K per segment except KeyError as e: self.logger.error("ACCFbackend: %s", e) diff --git a/src/pyneats/steps/climate_functions/local_accf.py b/src/pyneats/steps/climate_functions/local_accf.py index 9ce50f3..c4a159e 100644 --- a/src/pyneats/steps/climate_functions/local_accf.py +++ b/src/pyneats/steps/climate_functions/local_accf.py @@ -49,12 +49,17 @@ class LocalACCFParams(ACCFParams): col_geopotential: str = "geopotential" # [m^2 s^-2] col_time: str = "time" # datetime64 col_latitude: str = "latitude" # [deg] + # [K m^2 kg^-1 s^-1] (SI, pycontrails PotentialVorticity), not PVU + # (1 PVU = 1e-6 SI). Used as supplied; see h2o_raw_formula and + # docs/internal/pyneats_units_audit.md. col_potential_vorticity: str = "potential_vorticity" - col_nox_ei: str = "nox_ei" # [kg(NOx)/kg(fuel)] - col_fuel_burn: str = "fuel_burn" # [kg(fuel)] - col_air_pressure: str = "air_pressure" + # [kg(NOx)/kg(fuel)]; pycontrails NOx EI, NOx expressed as NO2 (ICAO EDB) + col_nox_ei: str = "nox_ei" + col_fuel_burn: str = "fuel_burn" # [kg(fuel)] per segment + col_air_pressure: str = "air_pressure" # [Pa] col_phase: str = "phase" + # Version-specific divisors applied to the raw regressions (dimensionless) scale_o3: Mapping[str, float] = field(default_factory=lambda: ACCF_SCALE_03) scale_ch4: Mapping[str, float] = field( default_factory=lambda: ACCF_SCALE_CH4, @@ -62,7 +67,7 @@ class LocalACCFParams(ACCFParams): scale_h2o: Mapping[str, float] = field( default_factory=lambda: ACCF_SCALE_H2O, ) - solar_constant: float = SOLAR_CONSTANT + solar_constant: float = SOLAR_CONSTANT # [W m^-2] def __post_init__(self) -> None: if self.accf_kwargs["accf_v"] not in self.scale_o3: @@ -94,7 +99,11 @@ def _fin_rowwise( lat_series: pd.Series | FloatArray, solar_constant: float, ) -> FloatArray: - """Top-of-atmosphere incoming shortwave parametrization.""" + """Top-of-atmosphere incoming shortwave parametrization [W m^-2]. + + ``lat_series`` is in degrees; the result has the units of + ``solar_constant`` (W m^-2). + """ times = np.asarray(pd.to_datetime(time_series), dtype="datetime64[ns]") doy = (times.astype("datetime64[D]") - times.astype("datetime64[Y]")).astype( "timedelta64[D]" @@ -145,15 +154,19 @@ def _require_cols(self, flight: Flight, cols: list[str]) -> None: # ---- O3: raw formula & compute ---- @staticmethod def o3_raw_formula(t: FloatArray, gp: FloatArray) -> FloatArray: - """Returns raw (unscaled, unclamped) P-ATR20-O3 [K/kg(NOx)].""" + """Returns raw (unscaled, unclamped) P-ATR20-O3 [K/kg(NOx)]. + + ``t`` in K, ``gp`` (geopotential) in m^2 s^-2. NOx mass is expressed + as NO2 (CLIMaCCF states K/kg(NO2)). + """ return -5.20e-11 + 2.30e-13 * t + 4.85e-16 * gp - 2.04e-18 * t * gp def compute_o3(self, flight: Flight) -> FloatArray: """Compute ozone (O3) Pulse-based Average Temperature Response over 20 years (P-ATR20) for ozone formation from NOx emissions. - Uses air temperature and geopotential as input parameters - + Uses air temperature and geopotential as input parameters. + Returns K per segment. """ t = _as_np(flight[self.params.col_air_temperature]) @@ -167,10 +180,10 @@ def compute_o3(self, flight: Flight) -> FloatArray: # Scaling with version-specific factor accf /= self.params.scale_o3[self.params.accf_kwargs["accf_v"]] - # convert from K-per-kg(NOx) to K-per-kg(fuel) with NOx EI + # convert from K-per-kg(NOx as NO2) to K-per-kg(fuel) with NOx EI accf *= _as_np(flight[self.params.col_nox_ei]) - # Final conversion from K-per-kg(fuel) to K + # Final conversion from K-per-kg(fuel) to K per segment accf *= _as_np(flight[self.params.col_fuel_burn]) except KeyError as e: @@ -183,7 +196,11 @@ def compute_o3(self, flight: Flight) -> FloatArray: # ---- CH4: raw formula & compute ---- @staticmethod def ch4_raw_formula(gp: FloatArray, fin: FloatArray) -> FloatArray: - """Returns raw (unscaled, unclamped) P-ATR20-CH4 [K/kg(NOx)].""" + """Returns raw (unscaled, unclamped) P-ATR20-CH4 [K/kg(NOx)]. + + ``gp`` (geopotential) in m^2 s^-2, ``fin`` (TOA incoming solar) in + W m^-2. NOx mass is expressed as NO2 (CLIMaCCF states K/kg(NO2)). + """ return -9.83e-13 + 1.99e-18 * gp - 6.32e-16 * fin + 6.12e-21 * gp * fin @@ -191,6 +208,7 @@ def compute_ch4(self, flight: Flight) -> FloatArray: """Compute methane (CH4) Pulse-based Average Temperature Response over 20 years (P-ATR20) for methane depletion from NOx emissions. Uses geopotential and top-of-atmosphere incoming shortwave radiation as input parameters. + Returns K per segment (PMO included when ``accf_kwargs["PMO"]`` is True). """ gp = _as_np(flight[self.params.col_geopotential]) @@ -208,10 +226,10 @@ def compute_ch4(self, flight: Flight) -> FloatArray: # Scaling with version-specific factor accf /= self.params.scale_ch4[self.params.accf_kwargs["accf_v"]] - # convert from K-per-kg(NOx) to K-per-kg(fuel) with NOx EI + # convert from K-per-kg(NOx as NO2) to K-per-kg(fuel) with NOx EI accf *= _as_np(flight[self.params.col_nox_ei]) - # Final conversion from K-per-kg(fuel) to K + # Final conversion from K-per-kg(fuel) to K per segment accf *= _as_np(flight[self.params.col_fuel_burn]) except KeyError as e: @@ -219,6 +237,7 @@ def compute_ch4(self, flight: Flight) -> FloatArray: f"Unknown ACCF version for CH4: {self.params.accf_kwargs['accf_v']}" ) from e if self.params.accf_kwargs["PMO"]: + # PMO effect merged into CH4 (dimensionless factor 1.29) accf = accf * 1.29 return accf @@ -226,14 +245,23 @@ def compute_ch4(self, flight: Flight) -> FloatArray: # ---- H2O: raw formula & compute ---- @staticmethod def h2o_raw_formula(pv: FloatArray) -> FloatArray: - """Returns raw (unscaled) P-ATR20-H2O [K/kg(fuel)].""" + """Returns raw (unscaled) P-ATR20-H2O [K/kg(fuel)]. + + ``pv`` is used as supplied, i.e. in SI units (K m^2 kg^-1 s^-1). The + CLIMaCCF reference implementation applies this regression to PV in + PVU (``pv * 1e6``). Without that factor the PV term here is ~1e-6 of + the CLIMaCCF value, so the result is effectively the constant term + (see docs/internal/pyneats_units_audit.md). + """ return 4.05e-16 + 1.48e-16 * np.abs(pv) def compute_h2o(self, flight: Flight) -> FloatArray: """Compute water vapor (H2O) Pulse-based Average Temperature Response over 20 years (P-ATR20) for water vapor emissions. - Uses potential vorticity as input parameter. + Uses potential vorticity as input parameter (SI units, not PVU; see + ``h2o_raw_formula`` and docs/internal/pyneats_units_audit.md). + Returns K per segment. """ pv = _as_np(flight[self.params.col_potential_vorticity]) @@ -242,7 +270,7 @@ def compute_h2o(self, flight: Flight) -> FloatArray: # Scaling with version-specific factor accf /= self.params.scale_h2o[self.params.accf_kwargs["accf_v"]] - # Final conversion from K-per-kg(fuel) to K + # Final conversion from K-per-kg(fuel) to K per segment accf *= _as_np(flight[self.params.col_fuel_burn]) except KeyError as e: @@ -259,8 +287,8 @@ def run(self, flight: FlightWithEmissions) -> FlightWithNonCO2Impact: In-tree reimplementation of the aCCF v1.0A regressions (O3, CH4, H2O), useful when CLIMaCCF is unavailable or pinned. Applies a validity mask based on air pressure: samples where pressure exceeds - ``DEFAULT_ACCF_VALIDITY_PRESSURE`` are set to NaN, since aCCF formulas - are only valid at cruise altitudes. + ``DEFAULT_ACCF_VALIDITY_PRESSURE`` (40000 Pa) are set to NaN, since aCCF + formulas are only valid at cruise altitudes. Args: flight: Flight with emissions; required columns are @@ -269,7 +297,8 @@ def run(self, flight: FlightWithEmissions) -> FlightWithNonCO2Impact: Returns: ``FlightWithNonCO2Impact`` — flight with ``ATR_20_O3``, - ``ATR_20_CH4`` and ``ATR_20_H2O`` columns (NaN outside validity). + ``ATR_20_CH4`` and ``ATR_20_H2O`` columns in K per segment + (P-ATR20 of that segment's emissions; NaN outside validity). Raises: ACCFStepError: If the local regression evaluation fails. diff --git a/src/pyneats/steps/climate_functions/open_airclim.py b/src/pyneats/steps/climate_functions/open_airclim.py index de30d82..4b95b81 100644 --- a/src/pyneats/steps/climate_functions/open_airclim.py +++ b/src/pyneats/steps/climate_functions/open_airclim.py @@ -48,10 +48,10 @@ class OpenAirClimParams(ClimateParams): base_inventory_path: str | None = str( files("pyneats.resources").joinpath("ELK_aviation_2025_res5deg_flat.nc") ) - grid_res: float = OPEN_AIRCLIM_GRID_RES - oac_levels: tuple[int, ...] = OPEN_AIRCLIM_LEVELS - computation_horizon: int = OPEN_AIRCLIM_COMPUTATION_HORIZON - horizons: tuple[int, ...] = METRICS_HORIZONS + grid_res: float = OPEN_AIRCLIM_GRID_RES # [degrees] lat/lon bin size + oac_levels: tuple[int, ...] = OPEN_AIRCLIM_LEVELS # [hPa] + computation_horizon: int = OPEN_AIRCLIM_COMPUTATION_HORIZON # [years] + horizons: tuple[int, ...] = METRICS_HORIZONS # [years] def _create_worker_toml( @@ -65,6 +65,15 @@ def _create_worker_toml( ): """ Generate OpenAirClim configuration file for a computation worker. + + ``time.range`` covers every year from ``year`` (the inventory year, i.e. + the year of the first waypoint) to ``year + horizon`` [years]. + ``metrics.t_0`` is kept from the base TOML (2025 in the bundled + ``openairclim.toml``), and ``OpenAirClimModel.run`` reads + ``AGWP_{H}_2025``. openairclim raises when a metric window + [t_0, t_0 + H) falls outside ``time.range``, so the shipped + configuration matches the emission year only for flights dated 2025 + (see docs/internal/pyneats_units_audit.md). """ with open(base_toml_path) as f: config = toml.load(f) @@ -125,6 +134,23 @@ class OpenAirClimModel( default_params = OpenAirClimParams def _flight_to_inventory(self, flight: Flight) -> xr.Dataset: + """Bin the flight into an openairclim emission inventory. + + Variables are sums over waypoints in each (lat, lon, plev) cell: + ``fuel`` (``fuel_burn``), ``CO2``, ``H2O`` and ``NOx`` (pycontrails + ``co2``/``h2o``/``nox``) in kg, from per-segment masses; ``distance`` + in km (``segment_length`` m / 1000). ``plev`` is the nearest + ``oac_levels`` entry to ``air_pressure`` (Pa / 100, i.e. hPa); + ``lat``/``lon`` are bin centres in degrees. ``Inventory_Year`` is the + year of the first waypoint; openairclim treats the inventory as that + year's emissions (no time-evolution file). + + NaN NOx, CO2 and H2O (e.g. NOx outside the met domain) are filled + with 0. NOx is the pycontrails NOx mass (expressed as NO2 by ICAO + EDB convention), while the bundled ``openairclim.toml`` sets + ``[species] nox = "NO"``, for which openairclim applies no 30/46 + correction (see docs/internal/pyneats_units_audit.md). + """ grid_res = self.params.grid_res df = flight.to_dataframe() @@ -186,6 +212,14 @@ def get_nearest_level(p): return ds def run(self, flight: FlightWithEmissions) -> FlightWithGlobalAGWP: + """Run openairclim on the flight inventory (Method D). + + Sets ``flight.attrs["AGWP_{H}_{O3,CH4,H2O,CONT}"]`` for H = 20, 50 and + 100 years, in W m-2 yr (openairclim "W m-2 year"), without efficacy, + read from openairclim's ``AGWP_{H}_2025`` (window starting at + ``metrics.t_0``; see ``_create_worker_toml``). ``AGWP_{H}_CH4`` + includes PMO. + """ with tempfile.TemporaryDirectory(dir=self.params.tmp_base_dir_path) as temp_dir: # 1. SETUP SYMLINK FOR REPOSITORY (Must be done before chdir if using absolute source paths) local_repo_link = os.path.join(temp_dir, "repository") diff --git a/src/pyneats/steps/climate_functions/views.py b/src/pyneats/steps/climate_functions/views.py index 7b54d69..73896b8 100644 --- a/src/pyneats/steps/climate_functions/views.py +++ b/src/pyneats/steps/climate_functions/views.py @@ -21,12 +21,15 @@ "FlightWithGlobalAGWP", ] +# P-ATR20 of each segment's emissions, per species [K per segment] REQUIRED_ATR_COLS: Final[tuple[str, ...]] = ( "ATR_20_CH4", "ATR_20_O3", "ATR_20_H2O", ) +# Flight-level AGWP at H = 20/50/100 years [W m-2 yr], without efficacy +# (openairclim "W m-2 year"); CH4 includes PMO. REQUIRED_AGWP_ATTRS: Final[tuple[str, ...]] = ( "AGWP_20_CH4", "AGWP_20_O3", @@ -43,8 +46,8 @@ ) REQUIRED_CONTRAIL_COLS: Final[tuple[str, ...]] = ( - "ef", - "sac", + "ef", # CoCiP energy forcing of the contrail formed at each waypoint [J] + "sac", # Schmidt-Appleman criterion flag: 1.0 met, 0.0 not met (NaN propagates) ) diff --git a/src/pyneats/steps/climate_metrics/gwp.py b/src/pyneats/steps/climate_metrics/gwp.py index a4836e0..311a295 100644 --- a/src/pyneats/steps/climate_metrics/gwp.py +++ b/src/pyneats/steps/climate_metrics/gwp.py @@ -45,7 +45,7 @@ # Supported species for aCCFs pathway SPECIES: Final[tuple[str, ...]] = ("CH4", "O3", "H2O") -# Column name pattern expected for ATR at H0 = 20 years +# Column name pattern expected for ATR at H0 = 20 years [K per segment] ATR_COL_TEMPLATE: Final[str] = "ATR_20_{spec}" @@ -60,9 +60,9 @@ class GWPParams(BaseParams): Parameters for GWP climate metrics computation. """ - horizons: tuple[int, ...] = METRICS_HORIZONS - surface_earth: float = SURFACE_EARTH - seconds_per_year: int = SECONDS_PER_YEAR + horizons: tuple[int, ...] = METRICS_HORIZONS # [years] + surface_earth: float = SURFACE_EARTH # [m^2] + seconds_per_year: int = SECONDS_PER_YEAR # [s/yr] agwp_coeff_wm2yr_per_kg: Mapping[int, float] = field( default_factory=lambda: CO2_AGWP_COEFF_WM2YR_PER_KG ) @@ -73,6 +73,7 @@ class GWPParams(BaseParams): default_factory=lambda: CONVERSION_FACTORS_ATR_TO_RF ) efficacy: Mapping[str, float] = field(default_factory=lambda: EFFICACY) + # H0 of the input ATR columns [years] atr_ref_horizon: int = field( default_factory=lambda: DEFAULT_CLIMACCF_KWARGS.get("time_horizon", 20) ) @@ -136,7 +137,8 @@ def _eagwp_contrails_Wm2yr(self, total_ef_J: float) -> dict[int, float]: Returns ------- dict[int, float] - EAGWP for each time horizon in W·m⁻²·yr. + EAGWP for each time horizon in W·m⁻²·yr (the same value for every + horizon). """ eps = float(self.params.efficacy.get("Contrails", 1.0)) s = self.params.surface_earth @@ -144,14 +146,17 @@ def _eagwp_contrails_Wm2yr(self, total_ef_J: float) -> dict[int, float]: return dict.fromkeys(self.params.horizons, total_ef_J * eps / (s * s_yr)) def _co2eq_from_eagwp_Wm2yr(self, agwp_Wm2yr: dict[int, float]) -> dict[int, float]: - """CO2eq(H) = EAGWP(H) / C(H)""" + """CO2eq(H) [kg CO2-eq] = EAGWP(H) [W m-2 yr] / C(H) [W m-2 yr kg-1]""" return { h: agwp_Wm2yr[h] / self.params.agwp_coeff_wm2yr_per_kg[h] for h in self.params.horizons } def _eagwp_spec_Wm2yr(self, species: str, atr_H0_K: float) -> dict[int, float]: - """Scales ATR(H0) [K] to AGWP(H) [W·m⁻²·yr]""" + """Scales ATR(H0) [K] to EAGWP(H) [W·m⁻²·yr]. + + EAGWP(H) = ATR(H0) * k_AGWP(H) / k_ATR(H0) / rf_backward * efficacy**2. + """ eps_spec = float(self.params.efficacy.get(species, 1.0)) h_0 = self.params.atr_ref_horizon rf_backward = float(self.params.rf_backward_factor.get(species, 1.0)) @@ -178,9 +183,12 @@ def _compute_from_segment_atr( ) -> list[dict[str, Any]]: """ Computations for Method C - 1. Sums 'ef' column for Contrails. - 2. Sums 'ATR_20_X' columns for Non-CO2. + 1. Sums 'ef' column for Contrails [J]. + 2. Sums 'ATR_20_X' columns for Non-CO2 [K]. 3. Applies conversion factors. + + Each entry's values hold ``EAGWP_Wm2yr`` [W m-2 yr] and ``CO2eq_kg`` + [kg CO2-equivalent]. """ df = flight.to_dataframe() results = [] @@ -246,9 +254,12 @@ def _compute_from_global_agwp( ) -> list[dict[str, Any]]: """ Computations for Method D - 1. Reads 'AGWP_{h}_{spec}' directly from attributes. + 1. Reads 'AGWP_{h}_{spec}' directly from attributes [W m-2 yr]. 2. Applies Efficacy. 3. Divides by CO2 contribution (CO2eq). + + Each entry's values hold ``EAGWP_Wm2yr`` [W m-2 yr] and ``CO2eq_kg`` + [kg CO2-equivalent]. """ results = [] horizons = self.params.horizons @@ -325,7 +336,7 @@ def run(self, flight: FlightWithNonCO2Impact) -> FlightWithClimateImpact: Returns: ``FlightWithClimateImpact`` — stores the final payload at - ``flight.attrs["climate_impact"]`` with shape:: + ``flight.attrs["climate_impact"]`` with shape (units below):: { "flight_information": { @@ -341,6 +352,24 @@ def run(self, flight: FlightWithNonCO2Impact) -> FlightWithClimateImpact: ], } + Units: + + - ``co2_baseline_kg``: flight total CO2 [kg] (``attrs["total_co2"]``, + the pycontrails nansum of per-segment ``co2``). + - ``fuel_burn_kg``: flight total fuel burned [kg] (nansum of the + per-segment ``fuel_burn`` column, else the ``fuel_burn`` or + ``fuel_burn_kg`` attr). + - ``contrails_ef_J``: Method C, sum of the CoCiP per-segment energy + forcing ``ef`` [J]. Method D stores ``attrs["AGWP_20_CONT"]`` + under this key, an AGWP in W m-2 yr without efficacy, not J + (see docs/internal/pyneats_units_audit.md). + - Each ``climate_metrics`` value: ``horizon`` [years], + ``EAGWP_Wm2yr`` (efficacy-weighted time-integrated forcing, + [W m-2 yr]) and ``CO2eq_kg`` (EAGWP / C(H), [kg CO2-equivalent]). + For CO2, ``CO2eq_kg`` is the total CO2 [kg]. CH4 includes PMO in + Method D and, with ``accf_kwargs["PMO"]`` True (default), in + Method C. + Raises: ClimateImpactStepError: If ``total_co2`` is missing, if the input view shape is inconsistent, or if the non-CO2/CO2 guardrail @@ -349,6 +378,7 @@ def run(self, flight: FlightWithNonCO2Impact) -> FlightWithClimateImpact: global-AGWP. """ # 1. Common Metadata & CO2 Baseline + # total_co2: flight total CO2 mass [kg] try: total_co2_kg = float(flight.attrs["total_co2"]) except Exception: @@ -423,6 +453,9 @@ def run(self, flight: FlightWithNonCO2Impact) -> FlightWithClimateImpact: self.logger.info( "Computed metrics using Flight-Level AGWP attributes (Integrated)" ) + # Stores AGWP_20_CONT, an AGWP in W m-2 yr (without efficacy), under + # the J-named contrails_ef_J key below + # (see docs/internal/pyneats_units_audit.md). total_ef_J = float(flight.attrs.get("AGWP_20_CONT", 0.0)) else: raise TypeError(f"Unknown input view type: {type(flight)}") diff --git a/src/pyneats/steps/climate_metrics/report.py b/src/pyneats/steps/climate_metrics/report.py index 7f84586..20ecfca 100644 --- a/src/pyneats/steps/climate_metrics/report.py +++ b/src/pyneats/steps/climate_metrics/report.py @@ -98,6 +98,8 @@ class FlightReport: Subset of PyContrails Flight.attrs relevant for *per-flight* reporting. This class defines the keys you care about; extraction returns a dict. + All fields are identifiers, labels or (``aobt``) a timestamp copied from + ``Flight.attrs``; none is a numeric physical quantity, so no units apply. Attributes ---------- @@ -112,7 +114,12 @@ class FlightReport: aircraft_id : str or None Aircraft identifier. aobt : Any - Actual off-block time (datetime or str). + Departure time labelled "actual off-block time", passed through as + found in ``Flight.attrs``. The AO field it comes from is the + flight-plan Item 13 time, i.e. the EOBT (ICAO Doc 4444). + The NEATS JSON parser sets it to the input ``departure_date_time`` + string, copied verbatim and never parsed + (e.g. ``"2025-07-09T18:07:00+0000"``). model_type : str or None Model type. aircraft_type : str or None diff --git a/src/pyneats/steps/climate_metrics/views.py b/src/pyneats/steps/climate_metrics/views.py index a46abe9..e55f8df 100644 --- a/src/pyneats/steps/climate_metrics/views.py +++ b/src/pyneats/steps/climate_metrics/views.py @@ -17,6 +17,10 @@ class FlightWithClimateImpact(FlightWithEmissions): Zero-copy view guaranteeing that attrs['climate_impact'] exists. Provides convenient accessors for the climate impact metadata and results. + + The payload is written by ``GWPMetrics.run``, which documents its shape + and units (``EAGWP_Wm2yr`` in W m-2 yr, ``CO2eq_kg`` in kg + CO2-equivalent, ``co2_baseline_kg`` and ``fuel_burn_kg`` in kg). """ ''' diff --git a/src/pyneats/steps/emissions/views.py b/src/pyneats/steps/emissions/views.py index e27e2de..4a78efa 100644 --- a/src/pyneats/steps/emissions/views.py +++ b/src/pyneats/steps/emissions/views.py @@ -17,6 +17,8 @@ "FlightWithEmissions", ] +# nvpm_ei_m [kg nvPM / kg fuel]; nox_ei [kg NOx / kg fuel], NOx expressed as +# NO2 (pycontrails, ICAO EDB convention). REQUIRED_EMISSION_COLS: Final[tuple[str, ...]] = ("nvpm_ei_m", "nox_ei") @@ -24,6 +26,21 @@ class FlightWithEmissions(FlightWithPerformance): """ Zero-copy typed view for emissions-enriched flights. + Only ``nvpm_ei_m`` and ``nox_ei`` are enforced. With the ``pycontrails`` + emissions backend the flight also carries, per waypoint: + + - ``co2``, ``h2o``, ``so2``, ``sulphates``, ``oc``, ``nox``, ``co``, + ``hc``, ``nvpm_mass``: mass emitted over the segment [kg], i.e. fuel + burned (whole-aircraft ``fuel_flow`` [kg/s] x segment duration [s]) + times the emission index. + - ``nvpm_number``: number of nvPM particles emitted over the segment. + - ``nox_ei``, ``co_ei``, ``hc_ei``, ``nvpm_ei_m``: emission indices + [kg/kg fuel]; ``nvpm_ei_n``: particles per kg fuel. + - ``attrs["total_"]``: flight totals (``np.nansum`` of the + per-segment values; kg, or a count for ``total_nvpm_number``). + Segments with NaN emission indices (e.g. NOx/CO/HC outside the met + domain) contribute nothing to the totals. + Attributes ---------- REQUIRED : ClassVar[tuple[str, ...]] diff --git a/src/pyneats/steps/interpolation/pycontrails_interpolation.py b/src/pyneats/steps/interpolation/pycontrails_interpolation.py index 01b2635..00c58c4 100644 --- a/src/pyneats/steps/interpolation/pycontrails_interpolation.py +++ b/src/pyneats/steps/interpolation/pycontrails_interpolation.py @@ -72,6 +72,13 @@ def run(self, flight: Flight4D) -> Flight4D: """ Interpolate a flight trajectory using pycontrails and linear interpolation for optional columns. + The latitude, longitude and altitude columns are resampled by + ``Flight.resample_and_fill(freq=interpolation_time)`` (default + ``"1min"`` = 60 s). Optional columns (fuel_flow, aircraft_mass, + engine_efficiency, true_airspeed) are linearly interpolated in time + with ``np.interp``; their units are unchanged. A NaN input value + makes every resampled sample strictly between its two neighbours NaN. + Parameters ---------- flight : Flight4D @@ -102,6 +109,7 @@ def run(self, flight: Flight4D) -> Flight4D: ) # Interpolate optional columns + # Time as int64 ns since epoch (pycontrails stores datetime64[ns]). t_orig_num = original["time"].values.astype(np.int64) t_new_num = interpolated_flight["time"].astype(np.int64) diff --git a/src/pyneats/steps/parsing/neats_io.py b/src/pyneats/steps/parsing/neats_io.py index b56ee73..6ddbdad 100644 --- a/src/pyneats/steps/parsing/neats_io.py +++ b/src/pyneats/steps/parsing/neats_io.py @@ -26,12 +26,18 @@ def neats_json_to_flights( per-flight DataFrames in the canonical Flight4D schema. This function assumes the JSON is already decoded (list[dict]). - It does *not* perform any file I/O. + It does *not* perform any file I/O. No unit conversion is applied here; + values are copied as received, except that ``ts`` is reformatted (see + Notes). Parameters ---------- flights : Sequence[Mapping[str, Any]] - Sequence of decoded NEATS/NM JSON flight objects. + Sequence of decoded NEATS/NM JSON flight objects. Must be a list (or + other sequence) of flight objects. The AO schema also allows a bare + single object at top level; it is not accepted here, because + iterating a dict yields its keys (see + docs/internal/pyneats_units_audit.md). model_type : str, optional Model type label to assign to each flight (default is "NM"). @@ -39,6 +45,77 @@ def neats_json_to_flights( ------- list of pd.DataFrame List of per-flight DataFrames in Flight4D schema. + + Notes + ----- + Field mapping (JSON field -> Flight4D name: unit as passed on). "MRR" + is Implementing Regulation (EU) 2018/2066 as amended by (EU) 2024/2493. + + Per-point columns, from ``flight_information.trajectory.trajectory_data[*]``: + + - ``ts`` -> ``time``: string ``"%Y-%m-%d %H:%M:%S"``. Parsed with + ``"%Y-%m-%dT%H:%M:%S%z"``; the UTC offset is parsed but discarded, not + applied, so a ``+02:00`` value keeps its local clock time and is later + read as UTC by ``NeatsTrajectoryParser`` (``utc=True``). Values must + already be UTC (``Z``, ``+0000`` or ``+00:00``), as the AO schema + requires. Fractional seconds, a lowercase ``z`` or a missing offset + raise ``ValueError`` (see docs/internal/pyneats_units_audit.md). Point + spacing is not assumed (the interpolation step resamples). + - ``lat`` / ``lon`` -> ``latitude`` / ``longitude``: decimal degrees + (WGS 84), unchanged. + - ``fl`` -> ``altitude``: still flight level (hundreds of ft). + ``NeatsTrajectoryParser`` converts it to metres of ISA pressure + altitude (FL x 100 ft x 0.3048 m/ft). + - ``tas`` -> ``true_airspeed``: m/s. The AO JSON field ``tas`` is in + knots; the NEATS application converts knots to m/s (x 1852 / 3600) + before calling PyNEATS, so values received here must already be in m/s. + - ``ff`` -> ``fuel_flow``: kg/s, whole aircraft. The AO schema does not + say so, but MRR Annex IIIa §1(33) defines fuel flow as the fuel mass + passing "into the aircraft's engines per second"; pycontrails + ``Emissions`` divides by ``n_engine`` to get per-engine flow. + - ``ee`` -> ``engine_efficiency``: dimensionless fraction [0, 1]. MRR + Annex IIIa §1(34) and the Commission NEATS guidance define it as the + (overall propulsion) efficiency, useful thrust relative to fuel energy input; + the regulation words it as a "percentage". pycontrails consumes it in + the SAC mixing-line slope as F * TAS / (fuel_flow * q_fuel). + - ``am`` -> ``aircraft_mass``: kg. + + Per-flight attributes (``df.attrs``), from ``flight_information``: + + - ``departure_date_time`` -> ``aobt``: copied verbatim as a string, + never parsed (used only in flight keys and outputs). The AO schema + describes it as the flight-plan Item 13 time in UTC, which ICAO + Doc 4444 defines as the estimated off-block time (EOBT); PyNEATS + labels the attr AOBT. + - ``arrival_date_time`` is mapped but then dropped (not a Flight4D attr). + ``nwp_date_time``, ``confidential``, ``ets_exemption`` and + ``trajectory.trj_data_source`` are not read; ``model_type`` comes from + the function argument. + - ``aircraft_properties.takeoff_mass`` -> ``takeoff_weight``: kg, the + mass "at beginning of the take-off run" (MRR Annex IIIa §1(29)). + - ``aircraft_properties.load_factor`` -> ``payload_factor``: + dimensionless fraction, used as payload mass / maximum payload (BADA + MPL, or MTOW - OEW when MPL is unavailable), not a seat load factor, + as defined in MRR Annex IIIa §1(32). + - ``fuel_properties.hydrogen_content`` -> ``hydrogen_content``: wt.% + (percent by mass, e.g. 13.79), not a 0-1 fraction. + - ``fuel_properties.hydrogen_per_carbon_ratio`` -> ``h_c_ratio``: used + by ``NEATSFuel`` as an atomic (molar) H/C ratio, only when + ``hydrogen_content`` is absent; MRR Annex IIIa §1(36) defines it as + hydrogen atoms per carbon atom. + - ``fuel_properties.calorific_value`` -> ``q_fuel``: J/kg of fuel (lower + calorific value). The AO JSON field ``calorific_value`` is in MJ/kg of + fuel; the NEATS application multiplies by 1e6 before calling PyNEATS, + so values received here must already be in J/kg. + - ``fuel_properties.aromatic_content`` -> ``aromatic_content`` and + ``fuel_properties.sulphur`` -> ``sulphur_content``: stored, not used in + calculations. Basis per the defaults in MRR Annex IIIa §5(4): + aromatics % by volume, sulphur % by mass. Whether the JSON value is a + percentage or a fraction is not defined. + - ``fuel_properties.naphthalene`` is emitted under the key + ``"naphthalene"``, but ``Flight4D.ATTRS_OPTIONAL`` spells it + ``"naphtalene"``, so the ``allowed_attrs`` filter drops the value here + (see docs/internal/pyneats_units_audit.md). """ result: list[pd.DataFrame] = [] @@ -64,6 +141,8 @@ def neats_json_to_flights( time_over = None else: # "2025-11-02T10:02:00+0000" → "2025-11-02 10:02:00" + # The UTC offset is parsed but discarded (not applied); the + # parser reads the result as UTC, so ts must already be UTC. dt = datetime.strptime(ts_raw, "%Y-%m-%dT%H:%M:%S%z") time_over = dt.strftime("%Y-%m-%d %H:%M:%S") @@ -79,6 +158,10 @@ def neats_json_to_flights( "fuel_flow": p.get("ff"), "engine_efficiency": p.get("ee"), "aircraft_mass": p.get("am"), + # PyNEATS input unit: m/s. The AO JSON field `tas` is in + # knots. The NEATS application converts knots to m/s + # (x 1852 / 3600) before calling PyNEATS. Values received + # here must already be in m/s. "true_airspeed": p.get("tas"), } ) @@ -115,6 +198,11 @@ def neats_json_to_flights( "hydrogen_content": fp.get("hydrogen_content"), "h_c_ratio": fp.get("hydrogen_per_carbon_ratio"), "aromatic_content": fp.get("aromatic_content"), + # PyNEATS input unit: J/kg of fuel (lower calorific value). + # The AO JSON field `calorific_value` is in MJ/kg of fuel. + # The NEATS application multiplies by 1e6 before calling + # PyNEATS. Values received here must already be in J/kg; do + # not convert again. "q_fuel": fp.get("calorific_value"), "sulphur_content": fp.get("sulphur"), "naphthalene": fp.get("naphthalene"), @@ -248,6 +336,16 @@ def split_df_into_flights( attr_columns : list[str] or None Columns that should be stored in df.attrs instead of standard columns. If None, defaults to all non-trajectory columns except `flight_key`. + + Notes + ----- + No unit conversion is applied. Trajectory columns use the same names and + units as the ``NeatsTrajectoryParser`` input (see + ``neats_json_to_flights``): ``altitude`` in flight level, ``time`` as + datetime64 (naive values are read as UTC, tz-aware values are converted + to UTC) or as a string in the parser ``date_format``. Per-flight columns + (e.g. ``q_fuel`` in J/kg, ``hydrogen_content`` in wt.%) are moved to + ``df.attrs``. """ flight_key_cols = ["flight_id", "arrival_airport", "departure_airport", "aobt"] diff --git a/src/pyneats/steps/parsing/neats_parser.py b/src/pyneats/steps/parsing/neats_parser.py index feafb71..a53ff46 100644 --- a/src/pyneats/steps/parsing/neats_parser.py +++ b/src/pyneats/steps/parsing/neats_parser.py @@ -32,6 +32,12 @@ # Maps CSV column names → Flight4D attr names. # "naphtalene" matches the existing (typo'd) key in ATTRS_OPTIONAL. +# NEATSFuel.from_attrs reads "naphthalene", so this value is stored but not +# read by NEATSFuel (see docs/internal/pyneats_units_audit.md). +# Values are converted with float() only (no unit conversion): q_fuel is +# expected in J/kg; h_c_ratio is used by NEATSFuel as an atomic H/C ratio; +# the fraction vs percentage scale of aromatics/sulphur/naphthalene is not +# established. _FUEL_COL_MAP: dict[str, str] = { "aromatics_content": "aromatic_content", "h_c_ratio": "h_c_ratio", @@ -53,6 +59,12 @@ def _load_airport_fuel_props( - JSON: ``{airport: float}`` (treated as q_fuel-only) or ``{airport: {attr: value, ...}}`` (multi-property) + Values are converted with ``float()`` only; no unit conversion is applied. + This file is read by PyNEATS directly, so the NEATS application's x 1e6 + conversion of the AO ``calorific_value`` does not apply: ``q_fuel`` must be + given in J/kg (bare JSON numbers are taken as ``q_fuel``). There is no + ``hydrogen_content`` column mapping. + Returns ``{}`` on missing/invalid file (with a warning logged). """ if path is None: @@ -128,18 +140,24 @@ class NeatsTrajectoryParserParams(TrajectoryParserParams): Attributes ---------- date_format : str - Format string for parsing dates. + Format string for parsing string timestamps. The default has no UTC + offset; values are parsed with ``utc=True``, i.e. read as UTC. + Strings that do not match become NaT and their rows are dropped. timezone : str - Output timezone; parsing is done as UTC then converted. + Timezone the parsed timestamps are converted to before the ``Flight`` + is built. It has no effect on the result: pycontrails ``Flight`` + converts any tz-aware time column to tz-naive UTC, so ``Flight.time`` + is always naive UTC (see docs/internal/pyneats_units_audit.md). airport_fuel_path : str | None Optional path to a CSV or JSON file mapping airport codes to fuel properties - (q_fuel, h_c_ratio, aromatic_content, sulphur_content, naphtalene). + (q_fuel in J/kg, h_c_ratio, aromatic_content, sulphur_content, + naphtalene). Values are used as-is, without unit conversion. Operator-provided values always take precedence; airport values fill only the gaps not supplied by the operator. """ date_format: str = "%Y-%m-%d %H:%M:%S" - timezone: str = "UTC" # output tz; parsing is done as UTC then converted + timezone: str = "UTC" # no effect on output: Flight.time is naive UTC airport_fuel_path: str | None = None @@ -162,7 +180,7 @@ class NeatsTrajectoryParser( Ensures: - All required columns are present - Numeric values are valid - - Timestamps are properly formatted and timezone-aware + - Timestamps are parsed as UTC (output ``Flight.time`` is tz-naive UTC) - No duplicate timestamps exist - No missing values in required columns """ @@ -182,16 +200,18 @@ def run(self, flight: pd.DataFrame) -> Flight4D: Performs the standard NEATS trajectory normalisation: 1. Verifies all required columns are present (``Flight4D.REQUIRED``). - 2. Converts altitude from flight-level (FL) to metres. - 3. Parses timestamps in ``params.date_format`` (tz-aware, converts to - ``params.timezone`` if set). + 2. Converts altitude from flight level (FL) to metres of ISA pressure + altitude (FL x 100 ft x 0.3048 m/ft). + 3. Parses timestamps in ``params.date_format`` as UTC. The output + ``Flight.time`` is tz-naive UTC whatever ``params.timezone`` is. 4. Cleans and deduplicates rows on the ``time`` column. 5. Rejects trajectories spanning more than 24 hours (likely a concatenation or timezone bug). 6. Builds flight attributes, optionally filling fuel properties from the per-airport fallback map if provided. 7. Rejects low-altitude trajectories that never cross the configured - ``max_pressure_level`` (no contribution to non-CO2 climate impact). + ``max_pressure_level`` in hPa (no contribution to non-CO2 climate + impact). Args: flight: Raw NEATS DataFrame containing the ``Flight4D.REQUIRED`` @@ -214,7 +234,9 @@ def run(self, flight: pd.DataFrame) -> Flight4D: if missing: raise TrajectoryParserStepError(f"missing required columns: {missing}") - # 2) FL → meters + # 2) FL → meters: FL x 100 ft x 0.3048 m/ft (= FL x 30.48 m). The + # result is ISA pressure altitude, not geometric altitude; + # pycontrails derives level/air_pressure from it with the ISA. try: alt_ft = ( pd.to_numeric(flight["altitude"], errors="coerce") @@ -229,7 +251,8 @@ def run(self, flight: pd.DataFrame) -> Flight4D: f"altitude conversion failed: {e}" ) from e - # 3) Parse time (tz-aware) + # 3) Parse time as UTC. Any tz conversion below is undone by + # pycontrails Flight, which stores time as tz-naive UTC. try: ts = pd.to_datetime( df["time"], @@ -282,6 +305,7 @@ def run(self, flight: pd.DataFrame) -> Flight4D: attrs[k] = attrs_input[k] # 5.5) Airport-based fuel-property fallback (operator values always win) + # Values are used as-is (no unit conversion; q_fuel in J/kg). if self._airport_fuel_props: airport = str(attrs.get("departure_airport", "")).upper() if airport in self._airport_fuel_props: @@ -312,6 +336,7 @@ def run(self, flight: pd.DataFrame) -> Flight4D: base = Flight(data=data_req, attrs=attrs, fuel=fuel_obj) # 8.5) Reject flights that never reach high enough altitude + # max_pressure_level and Flight.level are both in hPa. if self.params.max_pressure_level is not None: min_level = float(base.level.min()) if min_level > self.params.max_pressure_level: diff --git a/src/pyneats/steps/parsing/params.py b/src/pyneats/steps/parsing/params.py index 94bc655..1fcb44d 100644 --- a/src/pyneats/steps/parsing/params.py +++ b/src/pyneats/steps/parsing/params.py @@ -24,7 +24,9 @@ class TrajectoryParserParams(BaseParams): max_pressure_level : float | None Maximum pressure level (hPa) allowed for the flight's highest point. Flights whose *minimum* pressure level exceeds this threshold (i.e. never - reach high enough altitude) are rejected. ``500 hPa`` ≈ FL180. + reach high enough altitude) are rejected. The default, 550 hPa + (``DEFAULT_MAX_PRESSURE_LEVEL``), is about FL160 ISA pressure + altitude. Set to ``None`` to disable the filter. """ diff --git a/src/pyneats/steps/parsing/views.py b/src/pyneats/steps/parsing/views.py index 9179ce6..caccb52 100644 --- a/src/pyneats/steps/parsing/views.py +++ b/src/pyneats/steps/parsing/views.py @@ -21,27 +21,33 @@ ] REQUIRED_4D_COLS: Final[tuple[str, ...]] = ( - "latitude", - "longitude", + "latitude", # degrees (decimal, WGS 84) + "longitude", # degrees (decimal, WGS 84) + # Flight level (hundreds of ft) at parser input; metres of ISA pressure + # altitude after NeatsTrajectoryParser. "altitude", + # Parsed as UTC; tz-naive UTC in the parsed Flight. "time", ) OPTIONAL_COLS: Final[tuple[str, ...]] = ( - "fuel_flow", - "aircraft_mass", + "fuel_flow", # kg/s, treated as whole aircraft + "aircraft_mass", # kg + # Dimensionless fraction; pycontrails uses it as overall propulsion efficiency. "engine_efficiency", - "true_airspeed", + "true_airspeed", # m/s ) ATTRS_OPTIONAL: Final[tuple[str, ...]] = ( - "takeoff_weight", - "payload_factor", - "hydrogen_content", - "h_c_ratio", - "q_fuel", + "takeoff_weight", # kg + "payload_factor", # dimensionless fraction: payload / maximum payload + "hydrogen_content", # wt.% (percent by mass, e.g. 13.79) + "h_c_ratio", # read as atomic (molar) H/C; used only if hydrogen_content absent + "q_fuel", # J/kg of fuel (lower calorific value) "aircraft_series", "engine_uid", - "aromatic_content", - "sulphur_content", + "aromatic_content", # % by volume basis, scale undefined; unused in calculations + "sulphur_content", # % by mass basis, scale undefined; unused in calculations + # Spelled "naphtalene" here; neats_io emits and NEATSFuel.from_attrs reads + # "naphthalene" (see docs/internal/pyneats_units_audit.md). "naphtalene", ) @@ -51,7 +57,7 @@ "departure_airport", "arrival_airport", "model_type", - "aobt", + "aobt", # string copied verbatim from the input; never parsed ) diff --git a/src/pyneats/steps/performance/altitude_filter.py b/src/pyneats/steps/performance/altitude_filter.py index 49f075a..2119b46 100644 --- a/src/pyneats/steps/performance/altitude_filter.py +++ b/src/pyneats/steps/performance/altitude_filter.py @@ -27,8 +27,10 @@ def filter_low_altitude_points( ) -> FlightWithWeather: """Filter trajectory points below *min_altitude_fl* and return a new view. - At this pipeline stage altitude is already in **metres**, so the FL - threshold is converted: ``FL x 100 x 0.3048``. + At this pipeline stage the ``altitude`` column is already in **metres** + (pressure altitude; the NEATS parser converts FL x 100 ft to m), so the + FL threshold is converted: ``ft_to_m(FL x 100)`` = ``FL x 100 x 0.3048`` + m (the default FL20 is 609.6 m). Parameters ---------- @@ -37,8 +39,8 @@ def filter_low_altitude_points( min_altitude_fl : float Minimum Flight Level (inclusive). Default ``DEFAULT_MIN_ALTITUDE_FL``. max_filter_ratio : float - Maximum fraction of points that may be removed. If exceeded the - flight is considered un-recoverable. + Maximum fraction of points that may be removed (dimensionless + fraction [0, 1]). If exceeded the flight is considered un-recoverable. Raises ------ diff --git a/src/pyneats/steps/performance/bada_adapters.py b/src/pyneats/steps/performance/bada_adapters.py index 3abc4ce..7b8c3d6 100644 --- a/src/pyneats/steps/performance/bada_adapters.py +++ b/src/pyneats/steps/performance/bada_adapters.py @@ -133,6 +133,12 @@ class BaseBADAAdapter(AircraftProtocol, Protocol): - rocd_mps: meters/second (positive = climb) - accel_mps2: meters/second^2 (time derivative of TAS) - delta_tau: K (delta temperature vs ISA) + + ``thrust_fuel_segment`` returns (fuel flow [kg/s, whole aircraft], + thrust [N], phase, thrust segment label "LIDL"/"MCMB"/"TOTAL"). + ``v_stall_cas`` returns a CAS stall speed [m/s]. Metadata: MTOW, OEW and + MPL in kg. The ROCD phase threshold is given in ft/min and converted to + m/s at construction. """ def thrust_fuel_segment( @@ -173,7 +179,17 @@ def _phase_from_rocd(self, rocd_mps: float) -> FlightPhase: return "Cruise" def v_stall_cas(self, mass: float, config: str) -> float | None: - """Return CAS stall speed (m/s) or ``None`` when data is unavailable.""" + """Return CAS stall speed (m/s) or ``None`` when data is unavailable. + + ``mass`` in kg. BADA3: pyBADA returns kt2ms(Vstall) x + sqrt(mass / reference mass). BADA4: pyBADA 0.1.12 + ``FlightEnvelope.VStall`` has the signature + ``(mass, HLid, LG, nz=1.0, **kwargs)`` (BADA3: ``(mass, config)``), so + this call with ``config`` and no ``HLid``/``LG`` raises ``TypeError`` + for ``BADA4Adapter``; the caller in ``run_by_bada_version`` catches it + and continues without the speed filter + (see docs/internal/pyneats_units_audit.md). + """ result = self._obj.flightEnvelope.VStall(mass=mass, config=config) if result is None: return None @@ -514,6 +530,8 @@ def span(self) -> float | None: value = getattr(self._obj, "span", None) return cast("float | None", value) + # MTOW/OEW [kg]: from OPF files pyBADA converts tonnes to kg (x 1000) and + # sets MTOW = OPF maximum mass, OEW = OPF minimum mass. @property def MTOW(self) -> float | None: value = getattr(self._obj, "MTOW", None) @@ -521,7 +539,10 @@ def MTOW(self) -> float | None: @property def MPL(self) -> float | None: - return None # BADA3 has no MPL attribute ... + # Always None. pyBADA 0.1.12 Bada3Aircraft does set MPL (from OPF + # files: max payload x 1000, kg), but it is not read here, so mass + # estimation uses MTOW - OEW (see docs/internal/pyneats_units_audit.md). + return None @property def OEW(self) -> float | None: diff --git a/src/pyneats/steps/performance/bada_model.py b/src/pyneats/steps/performance/bada_model.py index 41a9bcb..a1a771f 100644 --- a/src/pyneats/steps/performance/bada_model.py +++ b/src/pyneats/steps/performance/bada_model.py @@ -146,9 +146,13 @@ class BADAPerformanceModelParams(PerformanceModelParams): Attributes ---------- true_air_speed_smoothing_window : int - Window size for true airspeed smoothing. + Savitzky-Golay window length (number of points) passed to pycontrails + ``segment_true_airspeed``; used only when TAS is computed from winds. reference_q_fuel : float - Reference fuel heat value. + Reference lower calorific value [J/kg of fuel]. BADA fuel flow is + divided by ``q_fuel_used / reference_q_fuel``. The default + (``REFERENCE_Q_FUEL`` = 43.13e6 J/kg) is the pycontrails JetA value, + not the per-aircraft BADA4 LHV (see docs/internal/pyneats_units_audit.md). delta_tau_compute_method : {"point", "zero"} Method for delta-tau computation. delta_tau_fill_method : {"bffill", "none", "zero"} @@ -162,15 +166,34 @@ class BADAPerformanceModelParams(PerformanceModelParams): bada3_root_path : str Root path for BADA3 data. payload_factor : float - Payload factor for mass estimation. + Payload mass as a fraction of the BADA maximum payload (MPL): + ZFW = OEW + payload_factor x MPL [kg], with MPL taken as MTOW - OEW + when the adapter gives no MPL (always the case with ``BADA3Adapter``). + The first-pass constant-mass guess is + OEW + payload_factor x (MTOW - OEW). Not a seat load factor; not + clipped. Overridden by the ``payload_factor`` flight attribute and + used only when neither an ``aircraft_mass`` column nor a + ``takeoff_weight`` attribute is given. fuel_reserve_fraction : float - Fuel reserve fraction for mass estimation. + Dimensionless fraction: fuel on board = trip fuel x + (1 + fuel_reserve_fraction) in the initial mass estimate. max_rel_mass_diff : float - Maximum relative mass difference for convergence. + Convergence threshold on the relative change of the initial mass + between iterations (dimensionless fraction). max_mass_estimation_iter : int Maximum number of mass estimation iterations. rocd_phase_threshold : float - Threshold for rate of climb/descent phase detection. + Rate of climb/descent threshold for phase detection [ft/min]; the + adapters convert it to m/s. + fuel_burn_threshold : float + Guardrail: the flight is rejected when total fuel burn [kg] exceeds + fuel_burn_threshold x (MTOW - OEW) [kg] (dimensionless multiple). + + Notes + ----- + ``ff_outlier_threshold`` [kg/s, whole aircraft] is an upper bound checked + only against the BADA-computed (q_fuel-corrected) fuel flow, never against + a provided ``fuel_flow`` column. """ true_air_speed_smoothing_window: int = DEFAULT_TRUE_AIR_SPEED_SMOOTHING_WINDOW @@ -303,7 +326,8 @@ def run(self, flight: FlightWithWeather) -> FlightWithPerformance: Fallback chain: 1. Normal BADA (4 or 3) execution (with speed filter applied) 2. Retry with BADA3 forced (if first attempt used BADA4) - 3. Altitude fallback — filter points below FL15, retry + 3. Altitude fallback — filter points below DEFAULT_MIN_ALTITUDE_FL + (FL20 = 2000 ft = 609.6 m), retry """ # --- Preprocessing + aircraft data extraction --- @@ -460,6 +484,8 @@ def run_by_bada_version( ) # 1b) speed filter — remove low-TAS points before perf computation + # segment_duration [s], rocd and acceleration are not recomputed + # after rows are removed (see docs/internal/pyneats_units_audit.md). try: df = filter_low_speed_points(df, adapter) except Exception: @@ -519,6 +545,7 @@ def _build_result( ) -> FlightWithPerformance: """Apply guardrails, compute efficiency, and build the output flight.""" # Fuel burn guardrail: reject if total fuel > (MTOW - OEW) * threshold + # Total fuel burn in kg (pandas sum skips NaN); MTOW - OEW in kg. if adapter.MTOW is not None and adapter.OEW is not None: useful_payload = adapter.MTOW - adapter.OEW total_fuel = float(df["fuel_burn"].sum()) @@ -566,6 +593,9 @@ def _extract_aircraft_attrs( engine_id: str | None = ( flight.attrs.get("engine_uid") if hasattr(flight, "attrs") else None ) + # J/kg (LHV), already converted from the AO MJ/kg upstream; see the + # NEATSFuel q_fuel docstring. Never None in practice: a pycontrails Fuel + # always carries q_fuel (NEATSFuel default DEFAULT_Q_FUEL = 42.8e6 J/kg). q_fuel: float | None = flight.fuel.q_fuel return icao, series, engine_id, q_fuel @@ -616,9 +646,15 @@ def _preprocess(self, flight: Flight) -> pd.DataFrame: def _compute_ground_and_true_airspeed( self, flight: Flight, window: int ) -> pd.DataFrame: - """Compute ground speed and true airspeed columns.""" + """Compute ground speed and true airspeed columns [m/s]. + + A provided ``true_airspeed`` column is used as is (m/s). Otherwise TAS + is computed by pycontrails ``segment_true_airspeed`` from the ground + speed (Savitzky-Golay smoothed over ``window`` points) and the + ``u_wind``/``v_wind`` columns (m/s). + """ df = flight.dataframe.copy(deep=False) - df["ground_speed"] = flight.segment_groundspeed() + df["ground_speed"] = flight.segment_groundspeed() # m/s if "true_airspeed" not in df.columns: u: NDArray[np.floating] = df["u_wind"].to_numpy(dtype=float, copy=False) @@ -631,6 +667,9 @@ def _compute_ground_and_true_airspeed( polyorder=1, ) else: + # PyNEATS input unit: m/s. The AO JSON field `tas` is in knots. + # The NEATS application converts knots to m/s (x 1852 / 3600) before + # calling PyNEATS. Values received here must already be in m/s. self.logger.debug("'true_airspeed' column already provided, keeping it") # ensure numeric dtype early @@ -638,7 +677,14 @@ def _compute_ground_and_true_airspeed( return df def _compute_kinematics(self, flight: Flight, df: pd.DataFrame) -> None: - """Compute kinematic columns: rocd, acceleration, segment_duration.""" + """Compute kinematic columns: rocd, acceleration, segment_duration. + + Forward differences: ``segment_duration`` [s] is t[i+1] - t[i] (NaN at + the last point); ``rocd`` [ft/min] comes from pycontrails + ``segment_rocd`` (converted to m/s at the adapter call); + ``acceleration`` [m/s^2] is the TAS forward difference (0 at the last + point, NaN replaced by 0). + """ df["segment_duration"] = flight.segment_duration() df["rocd"] = flight.segment_rocd() df["segment_duration"] = pd.to_numeric(df["segment_duration"], errors="coerce") @@ -657,7 +703,12 @@ def _compute_delta_tau( method: Literal["point", "zero"], fill: Literal["bffill", "zero", "none"], ) -> None: - """Compute delta-tau (air temperature deviation from ISA).""" + """Compute delta-tau (air temperature deviation from ISA) [K]. + + ``air_temperature`` [K] minus the ISA temperature at the ``altitude`` + column (pressure altitude, m). Passed to the adapters as pyBADA + ``deltaTemp``. + """ if method == "point": T_isa: NDArray[np.floating] = m_to_T_isa( df["altitude"].to_numpy(dtype=float, copy=False) @@ -680,7 +731,12 @@ def _early_exit_if_fuel_and_efficiency( engine_id: str, flight: Flight, ) -> FlightWithPerformance | None: - """If AO provided fuel_flow, engine_efficiency, and aircraft_mass, skip BADA model.""" + """If AO provided fuel_flow, engine_efficiency, and aircraft_mass, skip BADA model. + + ``fuel_burn`` = fuel_flow [kg/s, treated as whole aircraft] x + segment_duration [s], i.e. kg per forward segment. No fuel-burn + guardrail and no fuel-flow outlier check are applied on this path. + """ if ( "fuel_flow" in df.columns and "engine_efficiency" in df.columns @@ -707,10 +763,13 @@ def _derive_fuel_flow_from_mass_if_needed(self, df: pd.DataFrame) -> None: Derive fuel_flow from aircraft_mass trajectory if mass provided but fuel_flow not. When airlines provide aircraft_mass along the trajectory, we can infer fuel consumption - from the decrease in mass over time: - fuel_flow[i] = -(mass[i+1] - mass[i]) / segment_duration[i] + from the decrease in mass over time. The code computes: + fuel_flow[i] = -(mass[i] - mass[i-1]) / segment_duration[i] - where segment_duration[i] is in seconds, and fuel_flow is in kg/s. + i.e. a backward mass difference [kg] divided by the forward + segment_duration[i] = t[i+1] - t[i] [s], giving kg/s for the whole + aircraft, clipped at 0. fuel_flow[0] is 0 and fuel_flow[-1] is NaN + (segment_duration[-1] is NaN) (see docs/internal/pyneats_units_audit.md). The derived fuel_flow will be preserved in the output (like airline-provided fuel_flow), while engine efficiency will still be computed using BADA thrust and fuel_flow for @@ -730,7 +789,9 @@ def _derive_fuel_flow_from_mass_if_needed(self, df: pd.DataFrame) -> None: mass_vals = df["aircraft_mass"].values dt_vals = df["segment_duration"].values - # fuel_flow[i] = -(mass[i+1] - mass[i]) / segment_duration[i] + # mass_diff[i] = mass[i] - mass[i-1] (backward, kg), divided by the + # forward segment_duration[i] = t[i+1] - t[i] (s); last value is NaN + # (see docs/internal/pyneats_units_audit.md). # Use prepend to handle first waypoint (assumes no mass change at t=0) mass_diff = np.diff(mass_vals, prepend=mass_vals[0]) # First segment: diff=0 fuel_flow_derived = -mass_diff / dt_vals # Negative sign because mass decreases @@ -821,7 +882,15 @@ def _single_pass_performance( default_q_fuel: float, ) -> PerfOutput: """Single pass performance calculation with given or initial mass. - Now applies q_fuel correction during calculation.""" + Now applies q_fuel correction during calculation. + + With ``initial_mass`` [kg, at the first retained trajectory point] the + mass is integrated with the BADA fuel flow + (mass -= ff_corrected [kg/s] x segment_duration [s]), even when a + ``fuel_flow`` column is provided. Without it, the ``aircraft_mass`` + column [kg] is used point by point. Returned fuel flow is kg/s for the + whole aircraft and thrust is N. + """ if initial_mass is None and "aircraft_mass" not in df.columns: raise PerformanceStepError( @@ -837,7 +906,10 @@ def _single_pass_performance( max_consecutive_failures = self.params.max_consecutive_failures ff_outlier_threshold = self.params.ff_outlier_threshold - # Calculate correction ratio + # Calculate correction ratio (dimensionless: J/kg / J/kg). + # ff_corrected [kg/s] = ff_BADA / (q_fuel_used / reference_q_fuel), so + # ff_corrected x q_fuel_used = ff_BADA x reference_q_fuel (see + # ``reference_q_fuel`` in BADAPerformanceModelParams). q_fuel_ratio = q_fuel_used / default_q_fuel mass_arr: list[float] = [] @@ -867,8 +939,8 @@ def _single_pass_performance( _as_float(pt.true_airspeed), ft_to_m(_as_float(pt.rocd)) / 60, # Convert ft/min to m/s. The adapter requires SI - _as_float(pt.acceleration), - _as_float(pt.delta_tau), + _as_float(pt.acceleration), # m/s^2 + _as_float(pt.delta_tau), # K ) # Apply q_fuel correction immediately @@ -943,7 +1015,13 @@ def _estimate_initial_mass_iterative( q_fuel_used: float, default_q_fuel: float, ) -> PerfOutput: - """Iterative initial mass estimation using BADA performance model.""" + """Iterative initial mass estimation using BADA performance model. + + All masses in kg. Initial mass = min(MTOW, ZFW + trip fuel x + (1 + fuel_reserve_fraction)), with ZFW = OEW + payload_factor x MPL + (or x (MTOW - OEW) when MPL is None) and trip fuel = sum of BADA fuel + flow [kg/s] x segment_duration [s] (NaN skipped). + """ payload_factor: float | None = flight.attrs.get("payload_factor") @@ -1041,7 +1119,15 @@ def _compute_engine_efficiency_if_missing( q_fuel_used: float, ) -> None: """Compute engine_efficiency column if missing in flight. - Now uses the correct q_fuel_used value.""" + Now uses the correct q_fuel_used value. + + Overall propulsion efficiency (dimensionless) from pycontrails + ``overall_propulsion_efficiency``: thrust [N] x TAS [m/s] / + (BADA fuel flow [kg/s] x q_fuel_used [J/kg]), clipped to [0, 0.5]. + By design it uses the BADA fuel flow, not a provided ``fuel_flow``. + Runs only when the column is absent: NaN values in a provided column + are not filled, and provided values are not clipped. + """ if "engine_efficiency" in df.columns: self.logger.debug("'engine_efficiency' column already provided, keeping it") @@ -1063,6 +1149,14 @@ def _compute_engine_efficiency_if_missing( ) def _finalize_columns(self, df: pd.DataFrame, perf: PerfOutput) -> None: + """Write BADA outputs to ``df``. + + ``thrust`` [N]; ``aircraft_mass`` [kg] is the input column when given, + otherwise the mass integrated with the BADA fuel flow (the constant + first-pass guess if ``max_mass_estimation_iter`` <= 1); a provided or + mass-derived ``fuel_flow`` [kg/s] is kept as is; ``fuel_burn`` = + fuel_flow x segment_duration [kg per segment]. + """ df["thrust"] = perf["thrust"] df["aircraft_mass"] = perf["mass"] df["phase"] = perf["phase"] @@ -1086,6 +1180,9 @@ def _attach_output_attrs( nb_eng: int, engine_id: str, ) -> None: + # n_engine: engine count. pycontrails Emissions divides fuel_flow by + # n_engine to get per-engine flow, so fuel_flow is treated as + # whole-aircraft flow. out.attrs["n_engine"] = nb_eng out.attrs["wingspan"] = adapter.span out.attrs["bada_version"] = bada_version diff --git a/src/pyneats/steps/performance/speed_filter.py b/src/pyneats/steps/performance/speed_filter.py index 906d68d..70b2330 100644 --- a/src/pyneats/steps/performance/speed_filter.py +++ b/src/pyneats/steps/performance/speed_filter.py @@ -25,6 +25,10 @@ def filter_low_speed_points( ) -> pd.DataFrame: """Filter trajectory points with TAS below VStall. + Drops points whose ``true_airspeed`` [m/s] is below the BADA stall speed + [m/s, CAS] at MTOW in the "TO" configuration. The comparison is TAS + against CAS; the two coincide only near sea level. + Parameters ---------- df : pd.DataFrame @@ -32,7 +36,8 @@ def filter_low_speed_points( adapter : BaseBADAAdapter BADA adapter providing ``v_stall_cas()`` and ``MTOW``. max_filter_ratio : float - Maximum fraction of points that may be removed. + Maximum fraction of points that may be removed (dimensionless + fraction [0, 1]). Returns ------- diff --git a/src/pyneats/steps/weather/weather_factory.py b/src/pyneats/steps/weather/weather_factory.py index a35aa87..c139d91 100644 --- a/src/pyneats/steps/weather/weather_factory.py +++ b/src/pyneats/steps/weather/weather_factory.py @@ -7,8 +7,11 @@ Features -------- - Flexible caching strategies (disk/Zarr) -- Automatic unit conversions -- Variable standardization +- Variable standardization (DWD names mapped to pycontrails MetVariables) +- Limited DWD unit handling: level Pa -> hPa (heuristic on the first level), + clc percent -> fraction, and SDR x ``sdr_accumulate_dt_s`` when building the + Zarr cache. Every other variable is relabelled with the pycontrails unit + without converting its values (see docs/internal/pyneats_units_audit.md) - Chunked data access - Thread-safe operations - Comprehensive error handling @@ -113,7 +116,10 @@ class DWDZarrCacheSpec: met_chunks: Mapping[str, int] | None = None rad_chunks: Mapping[str, int] | None = None wind_chunks: Mapping[str, int] | None = None - # Unit fix: convert SDR [W m-2] to [J m-2] by multiplying by dt (seconds) + # Accumulation step [s]. build_cache multiplies the DWD 'sdr' values + # (assumed W m-2) by this and labels them 'J m**-2'. The live NetCDF + # path relabels them 'J m**-2' without multiplying, so the two paths + # differ by this factor (see docs/internal/pyneats_units_audit.md). sdr_accumulate_dt_s: int | None = ( DEFAULT_SDR_ACCUMULATE_DT_S # 1 hour by default; None to skip ) @@ -136,9 +142,9 @@ class WeatherFactoryParams: # for ERA5 and for DWD live reading data_dir: str - horizontal_resolution: float = DEFAULT_HORIZONTAL_RES_DEG - pressure_levels: tuple[float, ...] = DEFAULT_PRESSURE_LEVELS_HPA - weather_offset_hours: int = DEFAULT_WEATHER_OFFSET_H + horizontal_resolution: float = DEFAULT_HORIZONTAL_RES_DEG # degrees + pressure_levels: tuple[float, ...] = DEFAULT_PRESSURE_LEVELS_HPA # hPa + weather_offset_hours: int = DEFAULT_WEATHER_OFFSET_H # h, window length # unified optional cache config cache: WeatherCacheConfig | None = None @@ -147,6 +153,11 @@ class WeatherFactoryParams: chunks: Mapping[str, int] | None = None def time_bounds(self, dt: datetime) -> tuple[str, str]: + """Return ``(dt, dt + weather_offset_hours)`` as "%Y-%m-%d %H:%M:%S". + + ``strftime`` drops any time zone without converting, so ``dt`` must be + naive UTC or UTC-aware. + """ later = dt + timedelta(hours=self.weather_offset_hours) fmt = "%Y-%m-%d %H:%M:%S" return dt.strftime(fmt), later.strftime(fmt) @@ -239,11 +250,16 @@ class DWDFactory(WeatherFactoryProtocol): "geopot": Geopotential, "clc": CloudAreaFractionInAtmosphereLayer, "rhi": RelativeHumidity, + # pv: labelled K m**2 kg**-1 s**-1 (SI); values are not rescaled. "pv": PotentialVorticity, } _rad_map: ClassVar[Mapping[str, MetVariable]] = { "tsr": TOANetDownwardShortwaveFlux, "olr": TOAOutgoingLongwaveFlux, + # The DWD 'sdr' field is TOA incoming solar radiation (store history: + # swflx_dn, "solar incoming radiation at TOA") but is mapped to the + # pycontrails surface variable (J m**-2); see + # docs/internal/pyneats_units_audit.md. "sdr": SurfaceSolarDownwardRadiation, } @@ -259,6 +275,14 @@ def __init__(self, params: WeatherFactoryParams) -> None: def __call__( self, asofdate: datetime, hour: int | None = None ) -> WeatherProviderProtocol: + """Open (or build) the DWD weather for one forecast run. + + The run is selected from ``asofdate`` (date) and ``hour`` (run hour, + defaults to ``asofdate.hour``). The file names and, on the Zarr path, + the ``[asofdate, asofdate + weather_offset_hours]`` window are built + with ``strftime`` and no time-zone conversion, so ``asofdate`` must be + naive UTC or UTC-aware. + """ zc = self.params.cache.zarr if self.params.cache else None run_hour = asofdate.hour if hour is None else hour date_str = asofdate.strftime("%Y%m%d") @@ -402,7 +426,8 @@ def build_cache( os.makedirs(met_dir, exist_ok=True) os.makedirs(rad_dir, exist_ok=True) - # optional SDR accumulation (W m-2 -> J m-2) + # optional SDR accumulation: values (assumed W m-2) x + # sdr_accumulate_dt_s [s], labelled J m**-2 if (zc.sdr_accumulate_dt_s is not None) and ( "surface_solar_downward_radiation" in rad_ds_xr ): @@ -511,14 +536,20 @@ def _load_and_standardize_live( try: ds = xr.open_mfdataset(files, combine="by_coords") ds = self._standardize_dims(ds) - # Convert Pa → hPa on level coord if needed + # Convert Pa → hPa on level coord if needed (tests the first + # level only; the level units attribute is dropped, not set to hPa) if "level" in ds and float(ds["level"].values[0]) > 2000.0: ds = ds.assign_coords(level=ds["level"].astype("float64") / 100.0) # Unit fixes BEFORE mapping + # clc: percent -> dimensionless fraction [0, 1], unconditional + # (assumes the source is in percent; not checked). if "clc" in ds: ds["clc"] = ds["clc"] / 100.0 ds["clc"].attrs["units"] = "1" + # rhi: relabel only; '%' (or '-') becomes '1' without dividing by + # 100, so no mismatch warning follows (see + # docs/internal/pyneats_units_audit.md). if "rhi" in ds and ds["rhi"].attrs.get("units", "-") in ("-", "%"): ds["rhi"].attrs["units"] = "1" @@ -572,6 +603,12 @@ def _norm_units(u: str | None) -> str: def _standardize_vars( self, ds: xr.Dataset, var_map: Mapping[str, MetVariable] ) -> xr.Dataset: + """Rename DWD variables to pycontrails names and set their attrs. + + ``attrs['units']`` is overwritten with the MetVariable unit. On a unit + mismatch only a warning is logged; values are not converted (see + docs/internal/pyneats_units_audit.md). + """ for raw, mv in var_map.items(): if raw not in ds: logger.debug("%s not found in dataset", raw)