#!/usr/bin/env python3 """ atmprf2bufr.py -- UCAR/CDAAC atmPrf (netCDF) -> WMO BUFR Ed.4, template 3 10 026 Converts a CDAAC atmPrf radio-occultation file into a GNSS-RO BUFR message laid out like UCAR's operational bfrPrf product: * Step 1b : 247-level impact-height grid, 3 "frequencies" per level (L1 = 1.5 GHz, L2 = 1.2 GHz, ionosphere-corrected = 0), bending-angle error from Bend_ang_stdv, percent confidence from Bend_ang_conf * Step 2a : refractivity on the same grid (MSL altitude of each level) * Step 2b : optional, from the matching wetPf/wetPf2 file (--wet) * Step 2c : surface block, missing (percent confidence 0) like UCAR * Header : time = start_time (GPS->UTC), time increment = occpt_offset, centre of curvature rotated J2000 -> Earth-fixed at the occultation time, quality flags from irs / freq2 / bad * Master table version selectable (default 43); optional satellite sub-identifier 0 01 016 placed before 3 10 026 (needs v39+) Two encoders: --engine builtin pure-Python BUFR writer (default, no ecCodes needed) --engine eccodes ecCodes (needs tables for the chosen master version) Commands encode atmPrf.nc -o out.bufr [--wet wetPf2.nc] [--arh] [...] decode file.bufr (built-in decoder, prints a summary) compare new.bufr reference.bufr (field-by-field differences) Requires numpy + scipy (netCDF3 classic) or netCDF4; ecCodes optional. """ import argparse import datetime as dt import math import os import re import struct import sys import numpy as np # =========================================================================== # BUFR Table B parameters used by 3 10 026: code -> (scale, reference, width) # (WMO Table B; operator-modified occurrences are passed explicitly) # =========================================================================== P31 = -1073741824 TB = { 1016: (0, 0, 16), 1007: (0, 0, 10), 2019: (0, 0, 11), 1033: (0, 0, 8), 2172: (0, 0, 8), 25060: (0, 0, 14), 8021: (0, 0, 5), 4001: (0, 0, 12), 4002: (0, 0, 4), 4003: (0, 0, 6), 4004: (0, 0, 5), 4005: (0, 0, 6), 4006: (0, 0, 6), 33039: (0, 0, 16), 33007: (0, 0, 7), 27031: (2, P31, 31), 28031: (2, P31, 31), 10031: (2, P31, 31), 1041: (5, P31, 31), 1042: (5, P31, 31), 1043: (5, P31, 31), 2020: (0, 0, 9), 1050: (0, 0, 17), 4016: (0, -4096, 13), 5001: (5, -9000000, 25), 6001: (5, -18000000, 26), 10035: (1, 62000000, 22), 5021: (2, 0, 16), 10036: (2, -15000, 15), 31001: (0, 0, 8), 31002: (0, 0, 16), 2121: (-8, 0, 7), 7040: (1, 62000000, 22), 15037: (8, -100000, 23), 8023: (0, 0, 6), 7007: (0, -1000, 17), 15036: (3, 0, 19), 7009: (0, -1000, 17), 10004: (-1, 0, 14), 12001: (1, 0, 12), 13001: (5, 0, 14), 8003: (0, 0, 6), } REPLICATORS = (31001, 31002) def item(code, value, dscale=0, dwidth=0): """One data element with operator 2 01 / 2 02 changes applied.""" s, r, w = TB[code] return (code, value, s + dscale, r, w + dwidth) # =========================================================================== # Grid: the ROPP 247-level impact-height grid (ropp_thin_asglog-247.dat, # 10 cm precision), nudged by <=2.5 cm where needed so that float32(RoC + # undulation + h) reproduces UCAR's impact parameters exactly in all six # reference files (COSMIC-2 R03/G23, YAM-8 E06, GNOMES GN04/GN05, Spire S263). # --grid-file FILE reads any ROPP thinner .dat file instead. UCAR247 = np.array([ 2094.200, 2211.900, 2330.093, 2449.700, 2569.300, 2690.900, 2812.500, 2935.527, 3058.110, 3183.110, 3308.110, 3433.593, 3560.527, 3687.500, 3815.900, 3944.300, 4073.200, 4203.110, 4334.000, 4465.300, 4597.200, 4730.000, 4862.800, 4997.093, 5131.300, 5266.110, 5401.900, 5538.093, 5675.300, 5813.000, 5951.200, 6090.300, 6230.000, 6370.110, 6510.700, 6652.800, 6794.900, 6938.000, 7082.500, 7226.093, 7371.593, 7517.093, 7664.093, 7811.027, 7959.000, 8107.900, 8257.300, 8407.200, 8559.093, 8710.900, 8863.300, 9017.093, 9171.400, 9326.700, 9483.400, 9639.610, 9797.400, 9956.027, 10116.700, 10276.900, 10438.500, 10600.600, 10764.200, 10927.700, 11093.300, 11260.300, 11426.800, 11595.200, 11765.100, 11935.500, 12106.900, 12279.800, 12453.600, 12627.900, 12804.200, 12981.000, 13159.200, 13338.900, 13519.500, 13701.700, 13884.300, 14068.800, 14254.900, 14441.900, 14629.400, 14818.400, 15010.300, 15201.200, 15395.500, 15590.800, 15786.600, 15984.400, 16182.600, 16383.800, 16585.900, 16789.100, 16994.600, 17200.700, 17408.200, 17617.700, 17828.600, 18041.000, 18255.400, 18470.700, 18689.000, 18907.700, 19128.900, 19350.100, 19574.200, 19800.300, 20027.300, 20256.800, 20486.800, 20718.800, 20953.100, 21189.500, 21426.300, 21665.000, 21905.800, 22149.400, 22393.600, 22639.200, 22887.700, 23136.700, 23388.700, 23642.100, 23897.000, 24153.300, 24412.100, 24671.900, 24933.100, 25196.800, 25462.400, 25728.500, 25997.600, 26267.100, 26539.100, 26813.000, 27088.400, 27364.700, 27643.100, 27922.400, 28203.600, 28486.300, 28770.000, 29056.600, 29342.800, 29631.800, 29921.400, 30212.400, 30503.900, 30797.400, 31092.800, 31387.700, 31684.600, 31982.400, 32281.200, 32581.100, 32881.300, 33182.600, 33484.900, 33787.600, 34091.300, 34396.000, 34700.700, 35005.900, 35311.000, 35617.700, 35923.300, 36230.000, 36537.600, 36843.800, 37150.900, 37458.500, 37765.100, 38072.300, 38379.400, 38687.000, 38994.600, 39301.300, 39608.400, 39915.500, 40223.100, 40530.800, 40836.900, 41144.500, 41452.100, 41758.800, 42065.900, 42373.500, 42681.200, 42988.300, 43295.400, 43603.000, 43910.600, 44218.300, 44524.400, 44832.000, 45139.600, 45447.300, 45754.400, 46061.500, 46369.100, 46676.800, 46983.900, 47291.500, 47599.100, 47906.700, 48214.400, 48520.500, 48828.100, 49135.700, 49443.400, 49750.500, 50058.100, 50365.700, 50673.300, 50980.500, 51288.100, 51595.700, 51903.300, 52210.000, 52517.600, 52825.200, 53132.800, 53439.900, 53747.600, 54055.200, 54362.800, 54669.900, 54977.500, 55285.200, 55592.800, 55899.900, 56207.500, 56515.100, 56822.800, 57129.900, 57437.500, 57745.100, 58052.700, 58360.400, 58667.500, 58975.100, 59282.700, 59590.300, 59896.500, ]) # Low-pass FIR filter (20 m taps, -500..+500 m in impact height) that UCAR # applies to bending angles and refractivity on the native atmPrf grid before # sampling them at the 247 levels. Estimated by regularised least squares from # six atmPrf/bfrPrf pairs (2-12 km): symmetric, main lobe +-~210 m, small # negative side lobes, 50 % amplitude response at ~420 m vertical wavelength. # Error, confidence, lat/lon, azimuth and MSL altitude are NOT filtered. UCAR_FILTER = np.array([ +0.00062, +0.00085, +0.00104, +0.00101, +0.00058, -0.00033, -0.00170, -0.00343, -0.00537, -0.00719, -0.00848, -0.00890, -0.00819, -0.00615, -0.00258, +0.00284, +0.01042, +0.02029, +0.03213, +0.04516, +0.05832, +0.07048, +0.08065, +0.08816, +0.09270, +0.09417, +0.09257, +0.08796, +0.08048, +0.07045, +0.05843, +0.04535, +0.03236, +0.02052, +0.01055, +0.00279, -0.00276, -0.00632, -0.00823, -0.00886, -0.00849, -0.00728, -0.00547, -0.00343, -0.00158, -0.00021, +0.00061, +0.00097, +0.00101, +0.00086, +0.00061, ]) # PlanetiQ satellites (all fly Pyxis-A, instrument 534): LEO id -> # (satellite identifier 0 01 007, satellite sub-identifier 0 01 016) # GNOMES-1..-5 -> 267, sub-id 1..5 # GNOMES-10..-13 -> 268, sub-id 10..13 # YAM-8 -> 768, sub-id 8 PLANETIQ = {f"GN{n:02d}": (267, n) for n in range(1, 6)} PLANETIQ.update({f"GN{n:02d}": (268, n) for n in range(10, 14)}) PLANETIQ["YM08"] = (768, 8) PYXIS_A = 534 def planetiq(leo): """(satid, instrument, subid) for a PlanetiQ LEO id; exit if unknown.""" if leo not in PLANETIQ: sys.exit(f"PlanetiQ LEO id {leo!r} has no satellite id/sub-id assigned: " f"add it to PLANETIQ or pass --satid/--instrument/--subid") sat, sub = PLANETIQ[leo] return sat, PYXIS_A, sub # Mission prefix in CDAAC fileStamp -> # (satellite id C-5, instrument C-8[, sub-id 0 01 016, CCCC, Section 1 centre]) MISSIONS = { # PlanetiQ GNOMES (GNnn) and YAM (YMnn): ids/sub-ids from PLANETIQ, KWSP r"(?:GN|YM)\d\d": lambda m: planetiq(m.group(0)) + ("KWSP",), # Spire Lemur (CDAAC prefix Snnn): C-5 269, SGNOS-A 530, sub-id = Spire # FM number nnn, heading KWSS, Section 1 originating centre 178 (Spire) r"S(\d{3})": lambda m: (269, 530, int(m.group(1)), "KWSS", 178), r"C2E(\d)": lambda m: (749 + int(m.group(1)), 104), # COSMIC-2 E1..E6 r"PAZ": lambda m: (44, 104), r"TSX|TSRX": lambda m: (42, 103), r"TDX|TDMX": lambda m: (43, 103), r"KOM5": lambda m: (825, 103), r"GRCC|GRC1": lambda m: (803, 104), r"GRCD|GRC2": lambda m: (804, 104), r"MTPA|MTPB": lambda m: (3, 202), r"MTPC": lambda m: (5, 202), # Commercial: PlanetiQ GNOMES-A/Pyxis-A 267/534 (B: 268/535), Spire 269/530. # Their CDAAC prefixes and sub-ids vary; pass --satid/--instrument/--subid. } L2_MIN_IMPACT_HEIGHT = 20600.0 # m, between UCAR247 levels 112/113 REFR_MAX_IMPACT_HEIGHT = 59435.0 # m, between UCAR247 levels 244/245 GNSS_CLASS = {"G": 401, "R": 402, "E": 403, "C": 404, "J": 405} LEAP_SECONDS = 18 # GPS-UTC, valid since 2017-01-01 # =========================================================================== # netCDF reading # =========================================================================== def read_nc(path): """Return (attrs dict, vars dict of float arrays with NaN for fill).""" try: from scipy.io import netcdf_file f = netcdf_file(path, "r", mmap=False) attrs = {k: _attr(v) for k, v in f._attributes.items()} var = {k: np.array(v.data, dtype=float) for k, v in f.variables.items()} f.close() except Exception: import netCDF4 with netCDF4.Dataset(path) as f: attrs = {k: _attr(f.getncattr(k)) for k in f.ncattrs()} var = {k: np.ma.filled(np.ma.asarray(v[:], dtype=float), np.nan) for k, v in f.variables.items()} for k, a in var.items(): a[a <= -999.0] = np.nan return attrs, var def _attr(v): if isinstance(v, bytes): return v.decode(errors="replace").strip("\x00 ") a = np.asarray(v) return a.item() if a.size == 1 else a.astype(float) # =========================================================================== # Time and reference-frame helpers # =========================================================================== GPS0 = dt.datetime(1980, 1, 6) def gps_to_utc(gps_seconds, leap): return GPS0 + dt.timedelta(seconds=float(gps_seconds) - leap) _AS = math.pi / 180 / 3600 def _R1(a): c, s = math.cos(a), math.sin(a) return np.array([[1, 0, 0], [0, c, s], [0, -s, c]]) def _R2(a): c, s = math.cos(a), math.sin(a) return np.array([[c, 0, -s], [0, 1, 0], [s, 0, c]]) def _R3(a): c, s = math.cos(a), math.sin(a) return np.array([[c, s, 0], [-s, c, 0], [0, 0, 1]]) def j2000_to_ecef(t_utc): """IAU-76 precession + IAU-80 nutation (main terms) + GAST. ~cm level for vectors of a few km (enough for the centre of curvature).""" days = (t_utc - dt.datetime(2000, 1, 1, 12)).total_seconds() / 86400.0 T = (days + 69.184 / 86400.0) / 36525.0 # TT z = (2306.2181 * T + 1.09468 * T**2 + 0.018203 * T**3) * _AS zeta = (2306.2181 * T + 0.30188 * T**2 + 0.017998 * T**3) * _AS th = (2004.3109 * T - 0.42665 * T**2 - 0.041833 * T**3) * _AS P = _R3(-z) @ _R2(th) @ _R3(-zeta) om = math.radians(125.04452 - 1934.136261 * T) L = math.radians(280.4665 + 36000.7698 * T) Lp = math.radians(218.3165 + 481267.8813 * T) dpsi = (-17.20 * math.sin(om) - 1.32 * math.sin(2 * L) - 0.23 * math.sin(2 * Lp) + 0.21 * math.sin(2 * om)) * _AS deps = (9.20 * math.cos(om) + 0.57 * math.cos(2 * L) + 0.10 * math.cos(2 * Lp) - 0.09 * math.cos(2 * om)) * _AS eps = (84381.448 - 46.8150 * T) * _AS N = _R1(-(eps + deps)) @ _R3(-dpsi) @ _R1(eps) Tu = days / 36525.0 gmst = math.radians((280.46061837 + 360.98564736629 * days + 0.000387933 * Tu**2 - Tu**3 / 38710000.0) % 360) return _R3(gmst + dpsi * math.cos(eps)) @ N @ P # =========================================================================== # Profile preparation # =========================================================================== def interp_valid(xg, x, y): """Linear interpolation of y(x) at xg using only valid samples, no extrapolation (NaN outside the valid range).""" ok = ~(np.isnan(x) | np.isnan(y)) if ok.sum() < 2: return np.full(len(xg), np.nan) xs, ys = x[ok], y[ok] o = np.argsort(xs) xs, ys = xs[o], ys[o] out = np.interp(xg, xs, ys) out[(xg < xs[0]) | (xg > xs[-1])] = np.nan return out def ropp_sg_coeff(npoints, order): """Savitzky-Golay coefficients exactly as ropp_io_thin_sg (NR savgol, symmetric window of npoints left/right, polynomial order <= npoints).""" m = min(order, npoints) k = np.arange(-npoints, npoints + 1, dtype=float) A = np.vander(k, m + 1, increasing=True) b = np.zeros(m + 1) b[0] = 1.0 return A @ np.linalg.solve(A.T @ A, b) def ropp_io_thin_sg(val, npoints=1, order=1): """Python port of ROPP ropp_io_thin_sg.f90 (ROPP rev 3551). val: values sorted by level, NaN = missing. Returns the smoothed array with ROPP's masking: the first and last npoints+1 values and every value within npoints of a missing input are set missing (NaN).""" val = np.asarray(val, dtype=float) nn = len(val) c = ropp_sg_coeff(npoints, order) good = ~np.isnan(val) out = np.convolve(np.where(good, val, 0.0), c[::-1], "same") out[:1 + npoints] = np.nan out[nn - npoints - 1:] = np.nan for i in np.where(~good)[0]: out[max(0, i - npoints):min(nn, i + npoints + 1)] = np.nan return out def ropp_io_thin_fixed(full_lev, full_val, thin_lev, method="ASGLOG", npoints=None, sigma=False): """Python port of ROPP ropp_io_thin_fixed.f90 (ROPP rev 5714), including the caller's removal of missing values (ropp_io_thin_select). method: ASGLOG, ASGLIN, SGLOG, SGLIN, LOG, LIN. For ASG*, the S-G half-width is ROPP's adaptive choice unless npoints is given: fullres = |max(diff(full_lev))| (input order kept, as in ROPP) thinres = |min(diff(thin_lev))| nsmooth = 2*(max(2, floor(thinres/max(fullres,0.1)))//2); npoints = nsmooth/2 Returns (thin_val with NaN = missing, npoints used, sigmafac).""" ok = ~(np.isnan(full_lev) | np.isnan(full_val)) lev, val = np.asarray(full_lev, float)[ok], np.asarray(full_val, float)[ok] thin_lev = np.asarray(thin_lev, float) out = np.full(len(thin_lev), np.nan) if len(lev) < 2: return out, 0, 1.0 sigmafac, used = 1.0, 0 if method.startswith(("SGLOG", "SGLIN")): val, used = ropp_io_thin_sg(val, 1, 1), 1 elif method.startswith(("ASGLOG", "ASGLIN")): if npoints is None: fullres = abs(np.max(lev[1:] - lev[:-1])) thinres = abs(np.min(thin_lev[1:] - thin_lev[:-1])) nsmooth = max(2, int(np.floor(thinres / max(fullres, 0.1)))) nsmooth = 2 * (nsmooth // 2) npoints = nsmooth // 2 if sigma: sigmafac = math.sqrt(2 * npoints + 1.0) val, used = ropp_io_thin_sg(val, npoints, 2), npoints o = np.argsort(lev, kind="stable") lev, val = lev[o], val[o] log = method.startswith(("ASGLOG", "SGLOG", "LOG")) n = len(lev) for k, tl in enumerate(thin_lev): if tl < lev[0] or tl > lev[-1]: continue j = 1 while j < n - 1 and lev[j] < tl: j += 1 i = j - 1 thinint = abs(thin_lev[k] - thin_lev[k - 1]) if k == len(thin_lev) - 1 \ else abs(thin_lev[k + 1] - thin_lev[k]) if np.isnan(val[i]) or np.isnan(val[j]) or abs(lev[i] - lev[j]) < 0.001 \ or abs(lev[i] - lev[j]) > thinint: continue v1, v2 = val[i], val[j] lg = log and v1 > 0 and v2 > 0 if lg: v1, v2 = math.log(v1), math.log(v2) r = v2 + (tl - lev[j]) / (lev[i] - lev[j]) * (v1 - v2) out[k] = (math.exp(r) if lg else r) * sigmafac return out, used, sigmafac def sg_smooth(x, y, half=13, deg=2): """ROPP S-G smoothing of y(x) on the native samples (sorted by x).""" o = np.argsort(x) s = ropp_io_thin_sg(y[o], half, deg) r = np.empty_like(s) r[o] = s return r def interp_log(xg, x, y): """ROPP LOG interpolation: linear in log(y); falls back to linear where either bracketing value is <= 0. No extrapolation.""" ok = ~(np.isnan(x) | np.isnan(y)) if ok.sum() < 2: return np.full(len(xg), np.nan) xs, ys = x[ok], y[ok] o = np.argsort(xs) xs, ys = xs[o], ys[o] j = np.clip(np.searchsorted(xs, xg), 1, len(xs) - 1) w = (xg - xs[j - 1]) / (xs[j] - xs[j - 1]) a, b = ys[j - 1], ys[j] with np.errstate(all="ignore"): lg = np.exp(np.log(a) + w * (np.log(b) - np.log(a))) out = np.where((a > 0) & (b > 0), lg, a + w * (b - a)) out[(xg < xs[0]) | (xg > xs[-1])] = np.nan return out def smooth(x, y, width): """Smooth y(x) before interpolation. width = 'ucar' (UCAR_FILTER on the native 20 m grid, weights renormalised at data edges), 'none'/0, or a number = triangular window width in metres.""" if width in (None, 0, "0", "none"): return y if width == "ucar": o = np.argsort(x) v = y[o] ok = ~np.isnan(v) num = np.convolve(np.where(ok, v, 0.0), UCAR_FILTER[::-1], "same") den = np.convolve(ok * 1.0, UCAR_FILTER[::-1], "same") with np.errstate(invalid="ignore", divide="ignore"): f = np.where(ok, num / den, np.nan) out = np.empty_like(f) out[o] = f return out width = float(width) ok = ~(np.isnan(x) | np.isnan(y)) out = np.full_like(y, np.nan) xs, ys = x[ok], y[ok] for i in np.where(ok)[0]: w = np.clip(1 - np.abs(xs - x[i]) / width, 0, None) out[i] = (w * ys).sum() / w.sum() return out def nanv(v): return None if v is None or (isinstance(v, float) and math.isnan(v)) else v def prepare(atm_path, args): A, V = read_nc(atm_path) p = {"A": A} leap = args.leap # --- times -------------------------------------------------------------- t0 = gps_to_utc(A["start_time"], leap) p["t0"] = t0 p["toffset"] = A.get("occpt_offset") p["tocc"] = t0 + dt.timedelta(seconds=float(p["toffset"] or 0.0)) p["tnom"] = dt.datetime(int(A["year"]), int(A["month"]), int(A["day"]), int(A["hour"]), int(A["minute"])) # GTS heading # --- identification ------------------------------------------------------ stamp = str(A.get("fileStamp", os.path.basename(atm_path))) p["leo"] = re.split(r"[._]", stamp)[0] satid = instr = subid = cccc = hdr_centre = None for pat, fn in MISSIONS.items(): m = re.match(pat, stamp) if m: ids = fn(m) + (None, None, None) satid, instr, subid, cccc, hdr_centre = ids[:5] break p["centre"] = args.centre if args.centre is not None else (hdr_centre or 60) p["subid"] = args.subid if args.subid is not None else subid p["cccc"] = args.cccc or cccc or "KWBC" satid = args.satid if args.satid is not None else satid instr = args.instrument if args.instrument is not None else instr if satid is None or instr is None: sys.exit(f"unknown mission '{stamp}': give --satid and --instrument") p["satid"], p["instr"] = satid, instr m = re.search(r"\.([GREJC])(\d{2})(?:\.|$)", stamp) p["gclass"] = GNSS_CLASS[m.group(1)] if m else 401 p["prn"] = int(m.group(2)) if m else int(A.get("occulting_sat_id")) # --- quality flags (0 33 039, BUFR bit 1 = MSB) -------------------------- bits = set() p["bad"] = str(A.get("bad", "0")) not in ("0", "", "None") if p["bad"]: # UCAR: summary + bending/refr/met bits bits.update((1, 5, 6, 7)) if args.offline: bits.add(2) if int(A.get("irs", 1)) == -1: # CDAAC: irs = -1 -> rising bits.add(3) f2 = str(A.get("freq2", "")).upper() # set unless semi-codeless L2 P(Y)/W if f2 and f2 not in ("L2P", "L2W", "L2Y", "L2D"): bits.add(10) p["qcflags"] = sum(1 << (16 - b) for b in bits) # --- local Earth parameters ---------------------------------------------- p["roc"] = A["rfict"] * 1000.0 p["und"] = A["rgeoid"] * 1000.0 p["lat"], p["lon"] = A["lat"], A["lon"] p["azim"] = (A["azim"] + 180.0) % 360.0 # GNSS -> LEO line of sight curv = np.asarray(A.get("curv"), dtype=float) * 1000.0 if curv.size == 3 and args.curv_frame == "ecef": curv = j2000_to_ecef(p["tocc"]) @ curv p["curv"] = curv if curv.size == 3 else None if args.float32: # mimic UCAR precision for k in ("roc", "lat", "lon"): p[k] = float(np.float32(p[k])) # --- vertical grid --------------------------------------------------------- ih = V["Impact_height"] * 1000.0 # m above RoC+geoid base = A["rfict"] * 1000.0 + A["rgeoid"] * 1000.0 if args.grid_from: hg = reference_grid(args.grid_from) - base elif args.grid_file: txt = open(args.grid_file).read().split("Hlevels=")[1].split() hg = np.array([float(v) for v in txt]) elif args.grid == "ucar247": hg = UCAR247.copy() elif args.grid == "native": hg = np.sort(ih[~np.isnan(ih)]) else: dz = float(args.grid.split(":")[1]) hg = np.arange(math.ceil(np.nanmin(ih) / dz) * dz, args.top + 1, dz) a = base + hg if args.float32: a = a.astype(np.float32).astype(float) p["impact"] = a S = args.smooth if S in ("ropp-asglog", "ropp-asglin"): # Exact ROPP thin_fixed: adaptive npoints from the input/output level # spacings (or --sg-half if given explicitly), ROPP edge masking. npts = args.sg_half if args.sg_half_set else None meth = S.split("-")[1].upper() def thin(v): out, used, _ = ropp_io_thin_fixed(ih, v, hg, meth, npts) p.setdefault("ropp_npoints", used) if args.sg_edge == "ucar": fill = interp_valid(hg, ih, smooth(ih, v, "ucar")) out = np.where(np.isnan(out), fill, out) return out elif S in ("asglog", "sglog", "sglin"): # ROPP ropp_io_thin_sg (npoints=--sg-half, order=--sg-order), then LOG # (asglog/sglog) or linear (sglin) interpolation to the fixed levels. # Levels ROPP masks (within npoints+1 samples of the profile ends) are # filled from the empirical UCAR filter unless --sg-edge mask. def thin(v): s = sg_smooth(ih, v, args.sg_half, args.sg_order) out = (interp_valid if S == "sglin" else interp_log)(hg, ih, s) if args.sg_edge == "ucar": fill = interp_valid(hg, ih, smooth(ih, v, "ucar")) out = np.where(np.isnan(out), fill, out) return out else: thin = lambda v: interp_valid(hg, ih, smooth(ih, v, S)) p["bend"] = [thin(V[k]) for k in ("Bend_ang1", "Bend_ang2", "Bend_ang")] p["bend_err"] = interp_valid(hg, ih, V.get("Bend_ang_stdv", np.full_like(ih, np.nan))) p["pc"] = interp_valid(hg, ih, V.get("Bend_ang_conf", np.full_like(ih, np.nan))) p["tp_lat"] = interp_valid(hg, ih, V["Lat"]) p["tp_lon"] = interp_valid(hg, ih, V["Lon"]) p["tp_azim"] = (interp_valid(hg, ih, V["Azim"]) + 180.0) % 360.0 p["z"] = interp_valid(hg, ih, V["MSL_alt"] * 1000.0) p["refr"] = thin(V["Ref"]) # Empirical UCAR cut-offs (identical in both reference files checked): # L2 bending reported only above ~20.6 km impact height, refractivity # only up to ~59.4 km impact height. Disable with --no-ucar-cutoffs. if args.ucar_cutoffs: for b in p["bend"]: # UCAR drops bending below -1e-4 rad b[b < -1e-4] = np.nan p["bend"][1][hg < L2_MIN_IMPACT_HEIGHT] = np.nan p["refr"][hg > REFR_MAX_IMPACT_HEIGHT] = np.nan # --- optional Step 2b from wetPf/wetPf2 --------------------------------- p["l2b"] = read_wet(args.wet) if args.wet else None return p def read_wet(path): _, W = read_nc(path) def pick(*names): for n in names: if n in W: return W[n] return None gph = pick("gph", "Geop", "GPH", "geop") if gph is None: gph = pick("MSL_alt") gph = gph * 1000.0 if np.nanmax(gph) < 1000 else gph # km -> m pres, temp = pick("Pres"), pick("Temp") q = pick("sph", "shum", "Sph", "q") if temp is not None and np.nanmax(temp) < 150: # C -> K temp = temp + 273.15 if q is not None and np.nanmax(q) > 0.2: # g/kg q = q / 1000.0 o = np.argsort(gph) return dict(gph=gph[o], p=None if pres is None else pres[o] * 100.0, t=None if temp is None else temp[o], q=None if q is None else q[o]) # =========================================================================== # Sequence 3 10 026 (ROM SAF BUFR spec v2.8, Table 6) # =========================================================================== FREQS = (1.5e9, 1.2e9, 0.0) # as encoded by UCAR (7-bit, 1e8 Hz) def build(p, args): S = [] add = lambda *x: S.append(item(*x)) f = lambda v: nanv(float(v)) if v is not None else None t0 = p["t0"] if p["subid"] is not None: add(1016, p["subid"]) add(1007, p["satid"]); add(2019, p["instr"]); add(1033, args.proc_centre) add(2172, 2); add(25060, args.software_id); add(8021, 17) add(4001, t0.year); add(4002, t0.month); add(4003, t0.day) add(4004, t0.hour); add(4005, t0.minute) add(4006, t0.second + t0.microsecond * 1e-6, 3, 10) # 202131 201138 hdr_pc = p.get("overall_pc", 0 if p["bad"] else 100) add(33039, p["qcflags"]); add(33007, None if hdr_pc is None else int(round(hdr_pc))) v3 = lambda k: [f(x) for x in p[k]] if p.get(k) is not None else [None] * 3 for c, v in zip((27031, 28031, 10031, 1041, 1042, 1043), # LEO POD v3("leo_pos") + v3("leo_vel")): add(c, v) add(2020, p["gclass"]); add(1050, p["prn"]) for c, v in zip((27031, 28031, 10031), v3("gns_pos")): # GNSS, 202127 add(c, v, -1) for c, v in zip((1041, 1042, 1043), v3("gns_vel")): add(c, v) add(4016, f(p["toffset"]), 3, 5) # 202131 201133 add(5001, f(p["lat"])); add(6001, f(p["lon"])) cv = p["curv"] if p["curv"] is not None else [None] * 3 for c, v in zip((27031, 28031, 10031), cv): add(c, f(v)) add(10035, f(p["roc"])); add(5021, f(p["azim"])); add(10036, f(p["und"])) n1 = len(p["impact"]) add(31002, n1) for i in range(n1): add(5001, f(p["tp_lat"][i])); add(6001, f(p["tp_lon"][i])) add(5021, f(p["tp_azim"][i])) add(31001, 3) errs = p.get("bend_errs") or [p["bend_err"]] * 3 for k, hz in enumerate(FREQS): add(2121, hz); add(7040, f(p["impact"][i])) add(15037, f(p["bend"][k][i])); add(8023, 13) add(15037, f(errs[k][i]), 0, -3); add(8023, None) # 201125 pc = f(p["pc"][i]) add(33007, None if pc is None else int(math.floor(pc + 0.5))) add(31002, n1) rpc = p.get("refr_pc", p["pc"]) for i in range(n1): pc = f(rpc[i]) rerr = p.get("refr_err") add(7007, f(p["z"][i])); add(15036, f(p["refr"][i])); add(8023, 13) add(15036, None if rerr is None else f(rerr[i]), 0, -5) # 201123 add(8023, None) add(33007, None if pc is None else int(math.floor(pc + 0.5))) w = p["l2b"] n3 = 0 if w is None else len(w["gph"]) add(31002, n3) for i in range(n3): g = lambda k: None if w[k] is None else f(w[k][i]) add(7009, g("gph")); add(10004, g("p")); add(12001, g("t")) add(13001, g("q")); add(8023, 13) add(10004, None, 0, -8); add(12001, None, 0, -6) add(13001, None, 0, -5); add(8023, None) add(33007, 100) add(8003, 0); add(7009, None); add(10004, None); add(8023, 13) add(10004, None, 0, -8); add(8023, None); add(33007, p.get("l2c_pc", 0)) return S # =========================================================================== # Built-in BUFR Edition 4 writer # =========================================================================== class BitWriter: def __init__(self): self.v, self.n = 0, 0 def put(self, value, width): self.v = (self.v << width) | (value & ((1 << width) - 1)) self.n += width def bytes(self, even=True): pad = (-self.n) % 8 v, n = self.v << pad, self.n + pad b = v.to_bytes(n // 8, "big") if even and len(b) % 2: b += b"\0" return b _clip_log = {} def encode_value(code, value, scale, ref, width): missing = (1 << width) - 1 if value is None: return missing raw = int(math.floor(value * 10.0**scale + 0.5)) - ref if raw < 0 or raw >= missing: _clip_log[code] = _clip_log.get(code, 0) + 1 return missing return raw def write_builtin(seq, sec1, descriptors): bw = BitWriter() for code, value, scale, ref, width in seq: if code in REPLICATORS: bw.put(int(value), width) else: bw.put(encode_value(code, value, scale, ref, width), width) data = bw.bytes() s4 = (len(data) + 4).to_bytes(3, "big") + b"\0" + data def desc(d): f, x, y = d // 100000, (d // 1000) % 100, d % 1000 return ((f << 14) | (x << 8) | y).to_bytes(2, "big") body3 = b"\0" + (1).to_bytes(2, "big") + bytes([0x80]) + b"".join(map(desc, descriptors)) if (len(body3) + 3) % 2: body3 += b"\0" s3 = (len(body3) + 3).to_bytes(3, "big") + body3 s = sec1 body1 = (bytes([0]) + s["centre"].to_bytes(2, "big") + s["subcentre"].to_bytes(2, "big") + bytes([0, 0, 3, 50, 14, s["master"], 0]) + s["t"].year.to_bytes(2, "big") + bytes([s["t"].month, s["t"].day, s["t"].hour, s["t"].minute, s["t"].second])) s1 = (len(body1) + 3).to_bytes(3, "big") + body1 total = 8 + len(s1) + len(s3) + len(s4) + 4 return b"BUFR" + total.to_bytes(3, "big") + bytes([4]) + s1 + s3 + s4 + b"7777" # =========================================================================== # ecCodes writer (optional; key names discovered from the template itself) # =========================================================================== def write_eccodes(seq, sec1, descriptors): import eccodes as ec elems = [(c, v) for c, v, *_ in seq if c not in REPLICATORS] reps = [int(v) for c, v, *_ in seq if c in REPLICATORS] ext = [v for (c, v, *_) in seq if c == 31002] dly = [v for (c, v, *_) in seq if c == 31001] b = ec.codes_bufr_new_from_samples("BUFR4") try: t = sec1["t"] for k, v in dict(edition=4, masterTableNumber=0, bufrHeaderCentre=sec1["centre"], bufrHeaderSubCentre=sec1["subcentre"], updateSequenceNumber=0, dataCategory=3, internationalDataSubCategory=50, dataSubCategory=14, masterTablesVersionNumber=sec1["master"], localTablesVersionNumber=0, typicalYear=t.year, typicalMonth=t.month, typicalDay=t.day, typicalHour=t.hour, typicalMinute=t.minute, typicalSecond=t.second, numberOfSubsets=1, observedData=1, compressedData=0).items(): ec.codes_set(b, k, v) ec.codes_set_array(b, "inputExtendedDelayedDescriptorReplicationFactor", ext) if dly: ec.codes_set_array(b, "inputDelayedDescriptorReplicationFactor", dly) ec.codes_set_array(b, "unexpandedDescriptors", descriptors) exp = [int(d) for d in ec.codes_get_array(b, "expandedDescriptors")] exp = [d for d in exp if d < 100000 and d not in (31000,) + REPLICATORS] if exp != [c for c, _ in elems]: raise RuntimeError("ecCodes expansion differs from the built sequence") names = {} it = ec.codes_bufr_keys_iterator_new(b) while ec.codes_bufr_keys_iterator_next(it): k = ec.codes_bufr_keys_iterator_get_name(it) if "->" in k: continue try: names.setdefault(int(str(ec.codes_get(b, k + "->code"))), re.sub(r"^#\d+#", "", k)) except Exception: pass ec.codes_bufr_keys_iterator_delete(it) rank = {} for c, v in elems: n = names[c] rank[n] = rank.get(n, 0) + 1 key = f"#{rank[n]}#{n}" s, r, w = TB[c] if s == 0 and r == 0 and c not in (4006,): ec.codes_set_long(b, key, ec.CODES_MISSING_LONG if v is None else int(v)) else: ec.codes_set_double(b, key, ec.CODES_MISSING_DOUBLE if v is None else float(v)) ec.codes_set(b, "pack", 1) return ec.codes_get_message(b) finally: ec.codes_release(b) # =========================================================================== # GTS abbreviated heading (NCEP-style wrapper as in UCAR bfrPrf files) # =========================================================================== def area_letter(lat, lon): lon = (lon + 180.0) % 360.0 - 180.0 col = 0 if -90 <= lon < 0 else 1 if lon < -90 else 2 if lon >= 90 else 3 row = 0 if lat > 30 else 2 if lat < -30 else 1 return "ABCDEFGHIJKL"[row * 4 + col] # NCO-assigned GTS headings (TTAAii, CCCC) per LEO, matched on the LEO id # (CDAAC fileStamp prefix or ROPP leo_id). First match wins; edit as NCO # assigns new headings. Missions not listed fall back to the WMO-No. 386 # Table C3 area designator (IUT14) and the centre's CCCC. GTS_HEADINGS = [ (r"^GN05", "IUTM14", "KWSP"), # PlanetiQ GNOMES-5 (r"^GN04", "IUTN14", "KWSP"), # PlanetiQ GNOMES-4 (r"^YM08", "IUTO14", "KWSP"), # PlanetiQ YAM-8 (r"^S\d", "IUTM14", "KWSS"), # Spire (all satellites) (r"^C2E\d", "IUTM14", "KWBC"), # COSMIC-2 E1..E6 ] def gts_heading(leo, lat, lon, cccc_default, mode="nco"): """Return (TTAAii, CCCC) for a LEO id. mode='nco' uses GTS_HEADINGS, mode='wmo' always uses the Table C3 geographical designator.""" if mode == "nco" and leo: for pat, ttaaii, cccc in GTS_HEADINGS: if re.match(pat, leo.strip()): return ttaaii, cccc return f"IUT{area_letter(lat, lon)}14", cccc_default def wrap_arh(msg, lat, lon, t, cccc, ttaaii=None): """Wrap one BUFR message as a GTS bulletin with the 19-byte WMO/GTS flag field separator (FFS): '****nnnnnnnnnn****\n' 19 bytes, n = bytes that follow it 'TTAAii CCCC YYGGgg\r\r\n' 21 bytes abbreviated heading BUFR ... 7777 """ ttaaii = ttaaii or f"IUT{area_letter(lat, lon)}14" if not re.fullmatch(r"[A-Z]{4}\d{2}", ttaaii) or not re.fullmatch(r"[A-Z]{4}", cccc): raise ValueError(f"bad GTS heading {ttaaii!r} {cccc!r}") head = f"{ttaaii} {cccc} {t:%d%H%M}\r\r\n".encode("ascii") assert len(head) == 21 ffs = f"****{len(head) + len(msg):010d}****\n".encode("ascii") assert len(ffs) == 19 return ffs + head + msg def check_gts(data): """Parse and validate the FFS + abbreviated heading in front of a BUFR message. Returns a dict; 'ok' is False with a reason if anything is off.""" r = {"ok": True, "problems": []} i = data.find(b"BUFR") if i < 0: return {"ok": False, "problems": ["no BUFR message found"]} pre = data[:i] r["prefix_len"] = len(pre) if not pre: r["problems"].append("no GTS prefix (bare BUFR)") r["ok"] = False return r m = re.match(rb"\*\*\*\*(\d{10})\*\*\*\*\n", pre) if not m: r["problems"].append("first 19 bytes are not '****nnnnnnnnnn****\\n'") r["ok"] = False else: r["ffs_length"] = int(m.group(1)) r["actual_length"] = len(data) - 19 if r["ffs_length"] != r["actual_length"]: r["problems"].append(f"FFS length {r['ffs_length']} != bytes after FFS {r['actual_length']}") r["ok"] = False h = pre[19:] if m else pre mh = re.fullmatch(rb"([A-Z]{4}\d{2}) ([A-Z]{4}) (\d{6})( [A-Z]{3})?\r\r\n", h) if not mh: r["problems"].append(f"abbreviated heading not 'TTAAii CCCC YYGGgg': {h!r}") r["ok"] = False else: r["ttaaii"], r["cccc"], r["yygggg"] = (x.decode() for x in mh.groups()[:3]) r["heading_len"] = len(h) blen = int.from_bytes(data[i + 4:i + 7], "big") if data[i + blen - 4:i + blen] != b"7777": r["problems"].append("BUFR length / 7777 end marker inconsistent") r["ok"] = False if len(data) != i + blen: r["problems"].append(f"{len(data) - i - blen} extra bytes after 7777") return r # =========================================================================== # Built-in decoder (for verify / compare) # =========================================================================== # WMO Table B names/units of the elements used by 3 10 026 (+ 0 01 016) ELEMENT_INFO = { 1007: ("satelliteIdentifier", "Satellite identifier", "Code table"), 1016: ("satelliteSubIdentifier", "Satellite sub-identifier", "Numeric"), 1033: ("centre", "Identification of originating/generating centre", "Code table"), 1041: ("absolutePlatformVelocityFirstComponent", "Absolute platform velocity - first component", "m/s"), 1042: ("absolutePlatformVelocitySecondComponent", "Absolute platform velocity - second component", "m/s"), 1043: ("absolutePlatformVelocityThirdComponent", "Absolute platform velocity - third component", "m/s"), 1050: ("platformTransmitterIdNumber", "Platform transmitter ID number", "Numeric"), 2019: ("satelliteInstruments", "Satellite instruments", "Code table"), 2020: ("satelliteClassification", "Satellite classification", "Code table"), 2121: ("meanFrequency", "Mean frequency", "Hz"), 2172: ("productTypeForRetrievedAtmosphericGases", "Product type for retrieved atmospheric gases", "Code table"), 4001: ("year", "Year", "a"), 4002: ("month", "Month", "mon"), 4003: ("day", "Day", "d"), 4004: ("hour", "Hour", "h"), 4005: ("minute", "Minute", "min"), 4006: ("second", "Second", "s"), 4016: ("timeIncrement", "Time increment", "s"), 5001: ("latitude", "Latitude (high accuracy)", "deg"), 5021: ("bearingOrAzimuth", "Bearing or azimuth", "deg true"), 6001: ("longitude", "Longitude (high accuracy)", "deg"), 7007: ("height", "Height", "m"), 7009: ("geopotentialHeight", "Geopotential height", "gpm"), 7040: ("impactParameter", "Impact parameter", "m"), 8003: ("verticalSignificanceSatelliteObservations", "Vertical significance (satellite observations)", "Code table"), 8021: ("timeSignificance", "Time significance", "Code table"), 8023: ("firstOrderStatistics", "First-order statistics", "Code table"), 10004: ("pressure", "Pressure", "Pa"), 10031: ("distanceFromEarthCentreInDirectionOfNorthPole", "In direction of the North Pole, distance from the Earth's centre", "m"), 10035: ("earthLocalRadiusOfCurvature", "Earth's local radius of curvature", "m"), 10036: ("geoidUndulation", "Geoid undulation", "m"), 12001: ("airTemperature", "Temperature/air temperature", "K"), 13001: ("specificHumidity", "Specific humidity", "kg/kg"), 15036: ("atmosphericRefractivity", "Atmospheric refractivity", "N-units"), 15037: ("bendingAngle", "Bending angle", "rad"), 25060: ("softwareIdentification", "Software identification", "Numeric"), 27031: ("distanceFromEarthCentreInDirectionOf0DegreesLongitude", "In direction of 0 degrees longitude, distance from the Earth's centre", "m"), 28031: ("distanceFromEarthCentreInDirection90DegreesEast", "In direction 90 degrees East, distance from the Earth's centre", "m"), 31001: ("delayedDescriptorReplicationFactor", "Delayed descriptor replication factor", "Numeric"), 31002: ("extendedDelayedDescriptorReplicationFactor", "Extended delayed descriptor replication factor", "Numeric"), 33007: ("percentConfidence", "Per cent confidence", "%"), 33039: ("radioOccultationDataQualityFlags", "Quality flags for radio occultation data", "Flag table"), } # Meanings for code/flag table values that occur in RO messages CODE_MEANING = { 2020: {401: "GPS", 402: "GLONASS", 403: "Galileo", 404: "BeiDou", 405: "QZSS"}, 2172: {2: "limb sounding"}, 8021: {17: "start of phenomenon", 25: "nominal reporting time"}, 8023: {13: "root-mean-square"}, 8003: {0: "surface"}, 1033: {60: "NCAR/UCAR", 160: "NOAA/NESDIS", 94: "Copenhagen (DMI/ROM SAF)", 74: "Exeter (Met Office)", 78: "Offenbach (DWD/GFZ)", 254: "EUMETSAT", 38: "Beijing (CMA)", 28: "New Delhi", 178: "Spire"}, 1007: {750: "COSMIC-2 E1", 751: "COSMIC-2 E2", 752: "COSMIC-2 E3", 753: "COSMIC-2 E4", 754: "COSMIC-2 E5", 755: "COSMIC-2 E6", 3: "Metop-B", 4: "Metop-A", 5: "Metop-C", 44: "PAZ", 267: "PlanetiQ GNOMES-1..5", 268: "PlanetiQ GNOMES-10..13", 269: "Spire Lemur", 768: "PlanetiQ YAM-8", 825: "KOMPSAT-5"}, 2019: {104: "Tri-G", 202: "GRAS", 530: "SGNOS-A", 534: "Pyxis-A", 535: "Pyxis-B", 103: "IGOR"}, } RO_QC_BITS = {1: "non-nominal quality", 2: "offline product", 3: "rising occultation", 4: "excess phase non-nominal", 5: "bending angle non-nominal", 6: "refractivity non-nominal", 7: "meteorological product non-nominal", 8: "open-loop data", 9: "surface reflection", 10: "civil L2 signal (L2C)", 16: "missing"} def meaning(code, value): """Human-readable meaning of a code/flag table value (or '').""" if value is None: return "missing" if code == 33039: bits = [b for b in range(1, 17) if int(value) >> (16 - b) & 1] return ", ".join(f"{b}:{RO_QC_BITS.get(b, 'reserved')}" for b in bits) or "nominal (no bits set)" if code == 2121 and value is not None: return {1.5e9: "L1", 1.2e9: "L2", 0.0: "corrected (dummy)"}.get(value, "") return CODE_MEANING.get(code, {}).get(int(value), "") if code in CODE_MEANING else "" def bufr_messages(path): """Yield (offset, prefix_bytes, message_bytes) for every BUFR message in a file (bare, or each preceded by a GTS FFS/heading).""" d = open(path, "rb").read() pos, last_end = 0, 0 while True: i = d.find(b"BUFR", pos) if i < 0 or i + 8 > len(d): return n = int.from_bytes(d[i + 4:i + 7], "big") if n < 8 or d[i + n - 4:i + n] != b"7777": pos = i + 4 # 'BUFR' inside data, keep looking continue yield i, d[last_end:i], d[i:i + n] pos = last_end = i + n def decode(path, msg=0): """Decode message number msg (0-based) of a file; see decode_message.""" for k, (off, prefix, b) in enumerate(bufr_messages(path)): if k == msg: r = decode_message(b) r.update(prefix=prefix, offset=off, index=k) return r raise ValueError(f"{path}: message {msg} not found") def decode_message(b): """Decode one 3 10 026 message. Returns elements [(code, value)] and details [dict(code, value, raw, scale, ref, width, role, level, freq)].""" l1 = int.from_bytes(b[8:11], "big") o = 8 + l1 l3 = int.from_bytes(b[o:o + 3], "big") nsub = int.from_bytes(b[o + 4:o + 6], "big") flags = b[o + 6] descs = [int.from_bytes(b[o + j:o + j + 2], "big") for j in range(7, l3 - 1, 2)] descs = [(x >> 14) * 100000 + ((x >> 8) & 63) * 1000 + (x & 255) for x in descs] descs = [x for x in descs if x] if 310026 not in descs or nsub != 1 or flags & 0x40: raise ValueError(f"only uncompressed single-subset 3 10 026 supported " f"(descriptors {descs}, subsets {nsub}, flags {flags:#x})") o += l3 l4 = int.from_bytes(b[o:o + 3], "big") raw = int.from_bytes(b[o + 4:o + l4], "big") nbit, pos = (l4 - 4) * 8, [0] ctx = {"role": "", "level": None, "freq": None} def get(w): v = (raw >> (nbit - pos[0] - w)) & ((1 << w) - 1) pos[0] += w return v out, det = [], [] def rd(code, ds=0, dw=0, role=""): s, r, w = TB[code] s, w = s + ds, w + dw v = get(w) val = None if v == (1 << w) - 1 else ((v + r) / 10.0**s if s else v + r) out.append((code, val)) det.append(dict(code=code, value=val, raw=v, scale=s, ref=r, width=w, role=role or ctx["role"], level=ctx["level"], freq=ctx["freq"])) def rep(code, role): v = get(TB[code][2]) out.append((code, v)) det.append(dict(code=code, value=v, raw=v, scale=0, ref=0, width=TB[code][2], role=role, level=ctx["level"], freq=ctx["freq"])) return v ctx["role"] = "header" if 1016 in descs: rd(1016) for c in (1007, 2019, 1033, 2172, 25060, 8021, 4001, 4002, 4003, 4004, 4005): rd(c, role="start time" if 4001 <= c <= 4005 else "header") rd(4006, 3, 10, role="start time") rd(33039, role="summary quality"); rd(33007, role="summary quality") for c in (27031, 28031, 10031): rd(c, role="LEO position") for c in (1041, 1042, 1043): rd(c, role="LEO velocity") rd(2020, role="GNSS"); rd(1050, role="GNSS") for c in (27031, 28031, 10031): rd(c, -1, role="GNSS position") for c in (1041, 1042, 1043): rd(c, role="GNSS velocity") rd(4016, 3, 5, role="time of reference point") rd(5001, role="reference point"); rd(6001, role="reference point") for c in (27031, 28031, 10031): rd(c, role="centre of curvature") rd(10035, role="reference point"); rd(5021, role="line of sight azimuth") rd(10036, role="reference point") fname = {0: "L1", 1: "L2", 2: "corrected"} for lv in range(rep(31002, "Step 1b levels")): ctx.update(level=lv + 1, freq=None, role="1b tangent point") rd(5001); rd(6001); rd(5021) n0 = rep(31001, "1b frequencies") for k in range(n0): ctx.update(freq=fname.get(k, str(k + 1)) if n0 == 3 else "corrected") ctx["role"] = "1b" rd(2121); rd(7040); rd(15037, role="1b value"); rd(8023) rd(15037, 0, -3, role="1b error"); rd(8023, role="1b (stats off)") ctx.update(freq=None) rd(33007, role="1b confidence") ctx.update(level=None) for lv in range(rep(31002, "Step 2a levels")): ctx.update(level=lv + 1, role="2a") rd(7007); rd(15036, role="2a value"); rd(8023) rd(15036, 0, -5, role="2a error"); rd(8023, role="2a (stats off)") rd(33007, role="2a confidence") ctx.update(level=None) for lv in range(rep(31002, "Step 2b levels")): ctx.update(level=lv + 1, role="2b value") rd(7009); rd(10004); rd(12001); rd(13001); rd(8023) ctx["role"] = "2b error" rd(10004, 0, -8); rd(12001, 0, -6); rd(13001, 0, -5) rd(8023, role="2b (stats off)"); rd(33007, role="2b confidence") ctx.update(level=None, role="2c surface") for c in (8003, 7009, 10004, 8023): rd(c) rd(10004, 0, -8, role="2c error"); rd(8023, role="2c (stats off)") rd(33007, role="2c confidence") return dict(length=len(b), sec1=list(b[8:8 + l1]), descriptors=descs, elements=out, details=det, bits=pos[0], available=nbit) def reference_grid(path): """Impact parameters (m) of the Step 1b levels of an existing message.""" el = decode(path)["elements"] a, k = [], 0 while el[k][0] != 31002: k += 1 n1 = el[k][1] k += 1 for _ in range(n1): n0 = el[k + 3][1] a.append(el[k + 5][1]) k += 5 + 6 * n0 return np.array(a, dtype=float) # =========================================================================== # Commands # =========================================================================== def cmd_encode(args): p = prepare(args.input, args) if p["subid"] is not None and args.master < 39: sys.exit("satellite sub-identifier 0 01 016 requires master table >= 39") seq = build(p, args) descriptors = ([1016] if p["subid"] is not None else []) + [310026] sec1 = dict(centre=p["centre"], subcentre=args.subcentre, master=args.master, t=p["t0"]) writer = write_eccodes if args.engine == "eccodes" else write_builtin msg = writer(seq, sec1, descriptors) if args.arh: tt, cc = gts_heading(p.get("leo"), p["lat"], p["lon"], p["cccc"], args.heading) msg = wrap_arh(msg, p["lat"], p["lon"], p["tnom"], args.cccc or cc, args.ttaaii or tt) out = args.output or re.sub(r"^atmPrf_", "bfrPrf_", os.path.basename(args.input)).replace("_nc", "_bufr") with open(out, "wb") as fo: fo.write(msg) if "ropp_npoints" in p: print(f"ROPP S-G npoints = {p['ropp_npoints']} ({2 * p['ropp_npoints'] + 1}-point window)") print(f"wrote {out}: {len(msg)} bytes, {len(p['impact'])} Step-1b/2a levels, " f"{0 if p['l2b'] is None else len(p['l2b']['gph'])} Step-2b levels, " f"master table v{args.master}, engine={args.engine}") if _clip_log: print("values outside the BUFR range were set to missing:", _clip_log) def print_gts(path): g = check_gts(open(path, "rb").read()) if "ttaaii" in g: print(f"GTS heading: {g['ttaaii']} {g['cccc']} {g['yygggg']} " f"FFS length {g.get('ffs_length')} (bytes after FFS {g.get('actual_length')})") print("GTS wrapper:", "OK" if g["ok"] else "; ".join(g["problems"])) return g def cmd_gts_check(args): bad = 0 for f in args.files: print(f"--- {f}") bad += not print_gts(f)["ok"] sys.exit(1 if bad else 0) def cmd_decode(args): """Decode every message in a file and list every element with its WMO name, description, units, role (header / level / frequency) and value.""" import csv print_gts(args.file) rows = [] msgs = list(bufr_messages(args.file)) for k, (off, prefix, b) in enumerate(msgs): if args.message is not None and k != args.message: continue r = decode_message(b) s1 = r["sec1"] print(f"\n=== message {k + 1}/{len(msgs)} at byte {off}: {r['length']} bytes, " f"centre {s1[4] * 256 + s1[5]}, master table v{s1[13]}, " f"descriptors {r['descriptors']}, bits {r['bits']}/{r['available']}") if prefix.strip(): print(f" prefix {prefix!r}") kw, dw = (55, 66) if args.wide else (34, 40) cut = lambda t, n: t if len(t) <= n else t[:n - 1] + "~" if not args.csv: print(f"{'#':>5} {'desc':>6} {'name (ecCodes key)':<{kw}} {'description':<{dw}} " f"{'units':<10} {'level':>5} {'freq':<9} {'role':<22} {'value':>16} meaning") for n, e in enumerate(r["details"], 1): key, desc, units = ELEMENT_INFO.get(e["code"], ("?", "unknown", "")) if args.level is not None and e["level"] not in (None, args.level): continue mean = meaning(e["code"], e["value"]) if (units in ("Code table", "Flag table") or e["code"] in (2121,)) else "" rows.append(dict(message=k + 1, index=n, descriptor=f"{e['code']:06d}", key=key, description=desc, units=units, level=e["level"] or "", freq=e["freq"] or "", role=e["role"], value="missing" if e["value"] is None else e["value"], meaning=mean, raw=e["raw"], scale=e["scale"], reference=e["ref"], width=e["width"])) if not args.csv: val = "missing" if e["value"] is None else f"{e['value']:.10g}" print(f"{n:5d} {e['code']:06d} {cut(key, kw):<{kw}} {cut(desc, dw):<{dw}} " f"{units[:10]:<10} {e['level'] or '':>5} {e['freq'] or '':<9} " f"{e['role'][:22]:<22} {val:>16} {mean}") if args.csv: with open(args.csv, "w", newline="") as fo: w = csv.DictWriter(fo, fieldnames=list(rows[0].keys())) w.writeheader() w.writerows(rows) print(f"wrote {len(rows)} elements to {args.csv}") def split_blocks(el): """Split decoded elements into header / 1b / 2a / 2b / 2c blocks.""" idx = [i for i, (c, _) in enumerate(el) if c == 31002] names = ["header", "step1b", "step2a", "step2b"] cuts = [0] + idx blocks = {n: el[a:b] for n, a, b in zip(names, cuts, cuts[1:] + [None])} n3 = el[idx[2]][1] end = idx[2] + 1 + 10 * n3 blocks["step2b"], blocks["step2c"] = el[idx[2]:end], el[end:] return blocks def cmd_compare(args): A, B = decode(args.new), decode(args.ref) print("Section 1 new:", A["sec1"], "\n ref:", B["sec1"]) print("Prefix new:", A["prefix"], "\n ref:", B["prefix"]) ba, bb = split_blocks(A["elements"]), split_blocks(B["elements"]) for blk in ba: ea, eb = ba[blk], bb[blk] same = [c for c, _ in ea] == [c for c, _ in eb] print(f"\n[{blk}] {len(ea)} vs {len(eb)} elements" + ("" if same else " -- STRUCTURE DIFFERS, compared up to the shorter")) stats = {} for (c, va), (c2, vb) in zip(ea, eb): if c != c2: break s = stats.setdefault(c, [0, 0, 0, 0.0, 0.0]) s[0] += 1 if (va is None) != (vb is None): s[1] += 1 elif va is not None and va != vb: d = abs(va - vb) s[2] += 1 s[3] = max(s[3], d) if vb: s[4] = max(s[4], d / abs(vb)) print(f" {'desc':>6} {'count':>6} {'miss-mismatch':>13} {'val-diff':>8} {'max|diff|':>11} {'max rel':>9}") for c, (k, m, dv, mx, rl) in sorted(stats.items()): if m or dv: print(f" {c:06d} {k:6d} {m:13d} {dv:8d} {mx:11.5g} {rl:9.2e}") if all(not (m or dv) for _, m, dv, _, _ in stats.values()): print(" all compared values identical") def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) sp = ap.add_subparsers(dest="cmd", required=True) e = sp.add_parser("encode", help="atmPrf netCDF -> BUFR") e.add_argument("input") e.add_argument("-o", "--output") e.add_argument("--wet", help="matching wetPf/wetPf2 file for Step 2b") e.add_argument("--engine", choices=("builtin", "eccodes"), default="builtin") e.add_argument("--master", type=int, default=43, help="master table version") e.add_argument("--centre", type=int, help="Section 1 originating centre (default per mission, else 60)") e.add_argument("--subcentre", type=int, default=0) e.add_argument("--proc-centre", type=int, default=60, help="0 01 033") e.add_argument("--software-id", type=int, default=1, help="0 25 060") e.add_argument("--satid", type=int, help="0 01 007 override (C-5)") e.add_argument("--instrument", type=int, help="0 02 019 override (C-8)") e.add_argument("--subid", type=int, help="0 01 016 satellite sub-identifier") e.add_argument("--grid", default="ucar247", help="ucar247 | native | uniform:DZ (metres)") e.add_argument("--grid-from", help="copy the Step 1b grid from a reference BUFR") e.add_argument("--grid-file", help="ROPP thinner .dat file with Hlevels") e.add_argument("--sg-half", type=int, help="ROPP S-G npoints (half-width in native samples). Default 13 " "(27 pts) for sglin/sglog/asglog; ROPP's adaptive rule for " "ropp-asglog/ropp-asglin") e.add_argument("--sg-order", type=int, default=2, help="ROPP S-G polynomial order") e.add_argument("--sg-edge", choices=("ucar", "mask"), default="ucar", help="levels masked by ROPP S-G: fill from UCAR filter, or leave missing") e.add_argument("--top", type=float, default=60000.0, help="top for uniform grid (m)") e.add_argument("--smooth", default="ucar", help="filter for bending/refractivity before interpolation: " "ucar (default, empirical UCAR low-pass), sglin / sglog / asglog " "(ROPP ropp_io_thin_sg, 27 pts + linear / log interpolation), " "ropp-asglog / ropp-asglin (exact ROPP ropp_io_thin_fixed with its " "adaptive window), none, or a number = " "triangular width in m") e.add_argument("--curv-frame", choices=("ecef", "j2000"), default="ecef", help="frame for centre of curvature (UCAR uses Earth-fixed)") e.add_argument("--no-float32", dest="float32", action="store_false", help="don't round impact parameter, RoC, lat/lon to float32 like UCAR") e.add_argument("--no-ucar-cutoffs", dest="ucar_cutoffs", action="store_false", help="keep L2 below ~20.6 km and refractivity above ~59.4 km") e.add_argument("--leap", type=int, default=LEAP_SECONDS, help="GPS-UTC seconds") e.add_argument("--offline", action="store_true", help="set QC bit 2 (offline)") e.add_argument("--arh", action="store_true", help="add GTS heading like UCAR files") e.add_argument("--cccc", help="override CCCC of the GTS heading") e.add_argument("--ttaaii", help="override TTAAii of the GTS heading (e.g. IUTM14)") e.add_argument("--heading", choices=("nco", "wmo"), default="nco", help="nco: NCO-assigned headings (GTS_HEADINGS); wmo: Table C3 designator") e.set_defaults(func=cmd_encode) d = sp.add_parser("decode", help="list all elements with name/units/meaning") d.add_argument("file") d.add_argument("--message", type=int, help="only message N (0-based)") d.add_argument("--level", type=int, help="only header + profile level N (1-based)") d.add_argument("--wide", action="store_true", help="do not shorten long element names/descriptions") d.add_argument("--csv", help="write all elements to a CSV file instead") d.set_defaults(func=cmd_decode) g = sp.add_parser("gts-check", help="validate FFS + abbreviated heading") g.add_argument("files", nargs="+") g.set_defaults(func=cmd_gts_check) c = sp.add_parser("compare") c.add_argument("new") c.add_argument("ref") c.set_defaults(func=cmd_compare) args = ap.parse_args() import signal signal.signal(signal.SIGPIPE, signal.SIG_DFL) # quiet when piped to head/less if getattr(args, "cmd", "") == "encode": args.sg_half_set = args.sg_half is not None if args.sg_half is None: args.sg_half = 13 args.func(args) if __name__ == "__main__": main()