#!/usr/bin/env python3 """ ropp_nc2bufr.py -- ROPP-format netCDF (e.g. NESDIS RONA_L1C_*.nc, "atmPrf" in ROPP format) -> WMO BUFR 3 10 026, laid out like UCAR's bfrPrf product and with the fixes for the open issues B1-B11 / N2-N3 in STAR_RO_Open_Issues_3_Tables. It re-implements the ROPP chain (ropp_io_thin -> ropp_io_thin_select -> ropp_io_thin_fixed -> ropp_io_thin_sg, ConvertROPPtoBUFR, eum2bufr_eccodes). With --ropp-compat it reproduces today's NESDIS BUFR (verified bit-for-bit on the bending angles of RONA_L1C_C2E6_..._E26_A04); without it, the fixes below are applied: B1 satellite sub-identifier 0 01 016 before 3 10 026 (Spire Snnn -> nnn, PlanetiQ GNnn/YMnn -> nn), configurable B2 master table version 43 (default) B3 originating centre from processing_centre (160 NESDIS), --centre to override B4 GTS heading: NCO-assigned TTAAii/CCCC per mission (GTS_HEADINGS in atmprf2bufr.py), or --heading wmo for the Table C3 designator; 19-byte FFS '****nnnnnnnnnn****\n' in front B5 LEO/GNSS positions/velocities taken from r_leo_ref/v_leo_ref/r_gns_ref/ v_gns_ref when the netCDF provides them (N3 must be fixed upstream) B6 tangent-point azimuths interpolated on the circle (sin/cos) B7 tangent-point lat/lon interpolated as 3-D unit vectors (date line, poles) B8 levels masked by the S-G end effects are filled from the unsmoothed data B9 'offline' flag cleared for NRT (--nrt, default) ; non-nominal bits and percent confidence set from finalQA (N2); rising bit from irs; civil-L2 bit from the second frequency B10 software identification from software_version or --software-id B11 time increment from time_offset, else (time - start_time), allowed over the full BUFR range (-4.096 .. 258.046 s); outside it the reported start time is shifted so that start + increment is still the occultation time E1 missing bending-angle errors stay missing (ROPP writes 0.0) E2 247 levels are impact heights above RoC + undulation (ROPP grid / UCAR); ROPP's eum2bufr calls ropp_io_thin without impactalt=.TRUE. Usage python ropp_nc2bufr.py RONA_L1C_...nc -o out.bufr [--arh] [--thin ropp_thin_asglog-247.dat] python ropp_nc2bufr.py RONA_L1C_...nc --ropp-compat --master 39 -o ropp_like.bufr Yong Chen on 11/03/2026 This converter is working as expect and address EUMETSAT feedback (with fix) Final version: adding NCO header, using corrected bangle with generic instead of opt, and nrt as default python ropp_nc2bufr.py RONA_L1C_...nc -o out.bufr --heading nco --arh --corrected generic decode the bufr to check: python atmprf2bufr.py decode out.bufr > decode_out.txt Needs numpy + scipy (or netCDF4) and atmprf2bufr.py in the same directory. """ import argparse import datetime as dt import math import os import re import sys import numpy as np _argv, sys.argv = sys.argv, sys.argv[:1] import atmprf2bufr as ab # encoder, decoder, ROPP ports sys.argv = _argv EPOCH2000 = dt.datetime(2000, 1, 1) DEFAULT_GRID = os.path.join(os.path.dirname(os.path.abspath(__file__)), "ropp_thin_asglog-247.dat") # --- IDs (roppbufrcodes.nl defaults; --codes FILE reads the namelist) ------- LEO = {"C2E1": (750, 104), "C2E2": (751, 104), "C2E3": (752, 104), "C2E4": (753, 104), "C2E5": (754, 104), "C2E6": (755, 104), "META": (4, 202), "METB": (3, 202), "METC": (5, 202), "PAZ1": (44, 103), "TSRX": (42, 103), "TDMX": (43, 103), "KOM5": (825, 103), "GRAC": (803, 104), "GRAD": (804, 104), "FY3C": (522, 958), "FY3D": (523, 958), "FY3E": (524, 958), "SE6A": (66, 104), "S0": (269, 530), "S1": (269, 530), "S2": (269, 530), "GO01": (265, 526), "GO02": (265, 526), "GO03": (265, 526)} GNS = {"G": 401, "R": 402, "E": 403, "C": 404, "J": 405} ORG = {"NESDIS": (160, 0, "KNES"), "UCAR": (60, 0, "KWBC"), "DMI": (94, 0, "EKMI"), "METO": (74, 0, "EGRR"), "EUMETSAT": (254, 0, "EUMS"), "GFZ": (78, 173, "EDZW"), "CMA": (38, 0, "BAWX"), "ISRO": (28, 0, "DEMS")} # B1: PlanetiQ ids/sub-ids come from ab.PLANETIQ (GN01-05 -> 267/1-5, # GN10-13 -> 268/10-13, YM08 -> 768/8, all Pyxis-A 534); Spire Snnn -> 269, # SGNOS-A 530, sub-id = FM number nnn LEO.update({k: (v[0], ab.PYXIS_A) for k, v in ab.PLANETIQ.items()}) SPIRE_SUBID = r"^S(\d{1,3})$" def leo_subid(leo): """Satellite sub-identifier (0 01 016) for a LEO id, or None.""" if leo in ab.PLANETIQ: return ab.PLANETIQ[leo][1] mm = re.match(SPIRE_SUBID, leo) return int(mm.group(1)) if mm else None def read_namelist(path): """Minimal reader for roppbufrcodes.nl (LEOlist/LEOcode/InsCode, ORG...).""" txt = re.sub(r"!.*", "", open(path).read()) def arr(name, kind): mm = re.search(name + r"\s*=\s*(.*?)(?=\n\s*[A-Za-z]+\s*=|/)", txt, re.S) if not mm: return [] body = mm.group(1) return re.findall(r'"([^"]*)"', body) if kind == "s" else \ [int(x) for x in re.findall(r"-?\d+", body)] leo, lc, ic = arr("LEOlist", "s"), arr("LEOcode", "i"), arr("InsCode", "i") for n, a, b in zip(leo, lc, ic): if a >= 0: LEO[n.strip()] = (a, b) ol, oc, sc, ic2 = arr("ORGlist", "s"), arr("ORGcode", "i"), arr("Subcode", "i"), arr("ICAOcode", "s") for n, a, b, c in zip(ol, oc, sc, ic2): if a >= 0: ORG[n.strip()] = (a, max(b, 0), c) # --- netCDF ------------------------------------------------------------------- def read_ropp_nc(path): try: from scipy.io import netcdf_file f = netcdf_file(path, "r", mmap=False) attrs = {k: ab._attr(v) for k, v in f._attributes.items()} raw = {k: np.array(v.data) for k, v in f.variables.items()} f.close() except Exception: import netCDF4 with netCDF4.Dataset(path) as f: attrs = {k: ab._attr(f.getncattr(k)) for k in f.ncattrs()} raw = {k: np.ma.filled(v[:]) for k, v in f.variables.items()} V = {} for k, a in raw.items(): if a.dtype.kind in "SU": V[k] = b"".join(np.atleast_2d(a)[0]).decode(errors="ignore").strip("\x00 ") else: a = a.astype(float)[0] if a.ndim else a.astype(float) #V[k] = np.where(a < -9.0e4, np.nan, a) if np.ndim(a) else (np.nan if a < -9e4 else float(a)) V[k] = np.where(a <= -99900000.0, np.nan, a) if np.ndim(a) else ( np.nan if a <= -99900000.0 else float(a)) return attrs, V # --- thinning (ROPP port + fixes) -------------------------------------------- def thin(lev, val, grid, method, fix): """ropp_io_thin_select/fixed; with fix.edge_fill the S-G end masking (B8) is replaced by the unsmoothed values.""" out, npts, _ = ab.ropp_io_thin_fixed(lev, val, grid, method) if fix and method.startswith(("ASG", "SG")): lin, _, _ = ab.ropp_io_thin_fixed(lev, val, grid, "LOG" if "LOG" in method else "LIN") out = np.where(np.isnan(out), lin, out) return out, npts def thin_geoloc(lev, lat, lon, azi, grid, fix): if not fix: # ROPP non-polar branch (lon unwrap test, azimuth linear) tl, _ = thin(lev, lat, grid, "ASGLIN", False) lo = lon.copy() ok = ~np.isnan(lo) if np.any(np.abs(np.diff(lo[ok])) > 300): lo = np.where(lo < 0, lo + 360, lo) tn, _ = thin(lev, lo, grid, "ASGLIN", False) tn = np.where(tn > 180, tn - 360, tn) ta, _ = thin(lev, azi, grid, "ASGLIN", False) return tl, tn, ta la, lo, az = np.radians(lat), np.radians(lon), np.radians(azi) comp = [np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la), np.cos(az), np.sin(az)] x, y, z, c, s = (thin(lev, v, grid, "ASGLIN", True)[0] for v in comp) tl = np.degrees(np.arctan2(z, np.hypot(x, y))) tn = np.degrees(np.arctan2(y, x)) ta = np.degrees(np.arctan2(s, c)) % 360.0 return tl, tn, ta # --- main conversion ---------------------------------------------------------- def prepare(path, args): A, V = read_ropp_nc(path) fix = not args.ropp_compat p = {"bad": False} # identification leo, gns = V["leo_id"].strip(), V["gns_id"].strip() if re.match(r"^(GN|YM)\d\d$", leo) and leo not in ab.PLANETIQ: ab.planetiq(leo) # exits: no id/sub-id assigned key = leo if leo in LEO else leo[:2] satid, instr = LEO.get(key, (None, None)) p["satid"] = args.satid if args.satid is not None else satid p["instr"] = args.instrument if args.instrument is not None else instr if p["satid"] is None: sys.exit(f"unknown LEO id {leo!r}: add it to roppbufrcodes.nl or use --satid") p["gclass"] = GNS.get(gns[:1]) p["prn"] = int(gns[1:]) if gns[1:].isdigit() else None p["subid"] = None if args.subid is not None: p["subid"] = args.subid if args.subid >= 0 else None elif fix: p["subid"] = leo_subid(leo) pc_name = str(A.get("processing_centre", "UNKNOWN")).upper() org = next((v for k, v in ORG.items() if pc_name.startswith(k[:3])), (None, 0, "ZZZZ")) p["centre"] = args.centre if args.centre is not None else org[0] p["subcentre"], p["cccc"] = org[1], org[2] # times t0 = dt.datetime(int(V["year"]), int(V["month"]), int(V["day"]), int(V["hour"]), int(V["minute"]), int(V["second"])) + dt.timedelta(milliseconds=V["msec"]) toff = V.get("time_offset", np.nan) if fix and np.isnan(toff) and not np.isnan(V.get("time", np.nan)) \ and not np.isnan(V.get("start_time", np.nan)): toff = V["time"] - V["start_time"] # B11 lo_lim, hi_lim = (-4.096, 258.046) if fix else (0.0, 240.0) if not np.isnan(toff) and not lo_lim <= toff <= hi_lim: if fix: # B11 re-reference # smallest whole-second shift that brings the increment into range shift = math.ceil(toff - hi_lim) if toff > hi_lim else math.floor(toff - lo_lim) t0 += dt.timedelta(seconds=shift) toff -= shift else: toff = np.nan p["t0"], p["toffset"] = t0, (None if np.isnan(toff) else toff) p["tocc"] = t0 + dt.timedelta(seconds=float(p["toffset"] or 0)) # header geolocation p["lat"], p["lon"] = V["lat"], ((V["lon"] + 180) % 360) - 180 p["roc"], p["und"], p["azim"] = V["roc"], V["undulation"], V["azimuth"] p["curv"] = V["r_coc"] if np.all(np.isfinite(V["r_coc"])) else None for k, nc in (("leo_pos", "r_leo_ref"), ("leo_vel", "v_leo_ref"), ("gns_pos", "r_gns_ref"), ("gns_vel", "v_gns_ref")): v = np.asarray(V.get(nc, [np.nan] * 3), float) p[k] = v if np.all(np.isfinite(v)) and np.any(np.abs(v) > 1) else None # B5 # quality (ROPP PCD bit i -> BUFR 0 33 039 bit i+1, MSB first) ival = lambda k: int(V[k]) if np.isfinite(V.get(k, np.nan)) else 0 pcd = ival("pcd") bad = str(A.get("bad", "0")) not in ("0", "") or ival("finalQA") == 1 if fix: if args.nrt: pcd &= ~(1 << 1) # B9 offline if bad: pcd |= (1 << 0) | (1 << 4) | (1 << 5) | (1 << 6) # N2 like UCAR if ival("irs") == -1: pcd |= 1 << 2 # rising f2 = V.get("occfreq2", np.nan) civil = args.civil_l2 == "yes" or (args.civil_l2 == "auto" and not (abs(f2 - 1227.60e6) < 1e3 and gns[:1] == "G" and not re.match(r"^(C2|S\d|GN|YM)", leo))) if civil and np.isfinite(f2): pcd |= 1 << 9 # civil L2 signal p["qcflags"] = sum(1 << (15 - i) for i in range(16) if pcd >> i & 1) p["bad"] = bad if fix else False oq = V.get("overall_qual", np.nan) p["overall_pc"] = (0 if (fix and bad) else (None if np.isnan(oq) else oq)) # software id (B10) sw = args.software_id if sw is None: mm = re.search(r"(\d+(?:\.\d+)?)", str(A.get("software_version", ""))) sw = int(round(float(mm.group(1)) * 1000)) if mm else None args.software_id = sw # grid (E2) G = np.array([float(x) for x in open(args.thin).read().split("Hlevels=")[1].split()]) base = V["roc"] + (V["undulation"] if (fix and not args.grid_roc_only) else 0.0) grid = base + G p["impact"] = grid # Level 1b meth = args.method p["bend"], errs = [], [] for imp, ba, sg in (("impact_L1", "bangle_L1", "bangle_L1_sigma"), ("impact_L2", "bangle_L2", "bangle_L2_sigma"), ("impact_opt" if args.corrected == "opt" else "impact", "bangle_opt" if args.corrected == "opt" else "bangle", "bangle_opt_sigma" if args.corrected == "opt" else "bangle_sigma")): b, npts = thin(V[imp], V[ba], grid, meth, fix) p["bend"].append(b) s = V.get(sg) if s is None or np.all(np.isnan(s)): s = V.get("bangle_L1_sigma") # ROPP falls back to the L1 error s = np.where(s <= 0, np.nan, s) if not fix and sg == "bangle_L2_sigma": # ROPP: ASGLIN + sigma inflation (E3) e = ab.ropp_io_thin_fixed(V[imp], s, grid, "ASGLIN", sigma=True)[0] else: e, _ = thin(V[imp], s, grid, "LIN", False) errs.append(e if fix else np.where(np.isnan(e), 0.0, e)) # E1 p["bend_errs"], p["bend_err"], p["npoints"] = errs, errs[2], npts q = V.get("bangle_opt_qual" if args.corrected == "opt" else "bangle_qual") p["pc"] = thin(V["impact"], q, grid, "ASGLIN", fix)[0] if q is not None \ else np.full(len(grid), np.nan) p["tp_lat"], p["tp_lon"], p["tp_azim"] = thin_geoloc( V["impact"], V["lat_tp"], V["lon_tp"], V["azimuth_tp"], grid, fix) # Level 2a: N thinned on impact parameter a = (1 + 1e-6 N)(h + RoC + U) N, h = V["refrac"], V["alt_refrac"] a_n = (1 + 1e-6 * N) * (h + V["roc"] + V["undulation"]) p["refr"], _ = thin(a_n, N, grid, meth, fix) p["refr_err"] = thin(a_n, V.get("refrac_sigma"), grid, "LIN", False)[0] \ if V.get("refrac_sigma") is not None else None p["refr_pc"] = thin(a_n, V.get("refrac_qual"), grid, "ASGLIN", fix)[0] \ if V.get("refrac_qual") is not None else p["pc"] nz, _ = thin(a_n, N, grid, meth, fix) p["z"] = grid / (1 + 1e-6 * nz) - V["roc"] - V["undulation"] p["l2b"] = None p["l2c"] = None p["l2c_pc"] = 0 if fix else None # UCAR writes 0; ROPP writes missing p["leo"], p["gns"] = leo, gns return p def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("input") ap.add_argument("-o", "--output") ap.add_argument("--thin", default=DEFAULT_GRID, help="ROPP thinner .dat (247 levels)") ap.add_argument("--method", default="ASGLOG", choices=("ASGLOG", "ASGLIN", "SGLOG", "SGLIN", "LOG", "LIN")) ap.add_argument("--corrected", choices=("opt", "generic"), default="opt", help="corrected bending angle source (NESDIS uses bangle_opt)") ap.add_argument("--master", type=int, default=43) ap.add_argument("--centre", type=int) ap.add_argument("--cccc", help="override CCCC of the GTS heading") ap.add_argument("--ttaaii", help="override TTAAii of the GTS heading (e.g. IUTM14)") ap.add_argument("--heading", choices=("nco", "wmo"), default="nco", help="nco: NCO-assigned headings; wmo: Table C3 designator") ap.add_argument("--arh-time", choices=("start", "occ"), default="start", help="YYGGgg from occultation start (ROPP/NESDIS) or occultation time") ap.add_argument("--satid", type=int) ap.add_argument("--instrument", type=int) ap.add_argument("--subid", type=int, help="force sub-id (-1 = none)") ap.add_argument("--software-id", type=int, help="0 25 060 value") ap.add_argument("--nrt", action=argparse.BooleanOptionalAction, default=True, help="NRT product: clear the offline bit (default)") ap.add_argument("--civil-l2", choices=("auto", "yes", "no"), default="auto") ap.add_argument("--grid-roc-only", action="store_true", help="levels above RoC only (ROPP default, not UCAR)") ap.add_argument("--codes", help="roppbufrcodes.nl to load LEO/ORG tables from") ap.add_argument("--arh", action="store_true", help="add GTS heading like UCAR files") ap.add_argument("--ropp-compat", action="store_true", help="reproduce today's ROPP/NESDIS behaviour (no fixes)") args = ap.parse_args() if args.codes: read_namelist(args.codes) p = prepare(args.input, args) args.proc_centre = p["centre"] seq = ab.build(p, args) desc = ([1016] if p["subid"] is not None else []) + [310026] sec1 = dict(centre=p["centre"], subcentre=p["subcentre"], master=args.master, t=p["t0"]) msg = ab.write_builtin(seq, sec1, desc) if args.arh: tt, cc = ab.gts_heading(p["leo"], p["lat"], p["lon"], p["cccc"], args.heading) t_arh = p["t0"] if args.arh_time == "start" else p["tocc"] msg = ab.wrap_arh(msg, p["lat"], p["lon"], t_arh, args.cccc or cc, args.ttaaii or tt) out = args.output or os.path.basename(args.input).rsplit(".", 1)[0] + ".bufr" open(out, "wb").write(msg) print(f"wrote {out}: {len(msg)} bytes | LEO {p['leo']} -> {p['satid']}" f"{'' if p['subid'] is None else ' sub-id ' + str(p['subid'])} | " f"S-G npoints {p['npoints']} | QC flags {p['qcflags']} | master v{args.master}") if __name__ == "__main__": main()