From 0ea148097c827215a12ae15ba9689415b2db32f6 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Thu, 10 Sep 2026 13:56:09 -0500 Subject: [PATCH] Remove the PhiDP system phase offset before Bringi KDP processing The BNF C-SAPR2 raw differential phase sits on a system phase offset of roughly 220 degrees. The bringi branch passed it to csu_kdp unchanged, so corrected_differential_phase kept the offset, and Py-ART's _prepare_phidp smeared it along each ray with np.maximum.accumulate. Z-PHI then read that plateau as propagation phase: the self consistency number came out at 10.34 instead of order 0.05 to 0.2, and the reflectivity fed to the integral was inflated by about 16 dB. The result was specific differential attenuation up to 2.2 dB/km, roughly 20 times the physical C-band value, with radial streaking where the plateau differed between adjacent rays. Estimate the offset and subtract it, undoing the 0-360 fold in the same step, before calling calc_kdp_bringi. get_sys_phase takes the median over the first meteorological gates of each ray then the median across rays; Py-ART's estimators take a per-ray minimum, which sits well below the true offset when PhiDP noise is order 30 degrees, and both raise IndexError on the BNF RHI files whose sweep_end_ray_index disagrees with the ray count. They are kept as fallbacks. The value used is recorded as a system_phase_offset attribute and can be pinned with the new phidp_sys_phase configuration key. On bnfcsapr2cfrS3.a1.20250520.000009 at 3.5 degrees this takes corrected_differential_phase from 196..257 to -26..36 degrees, ADP max from 2.206 to 0.100 dB/km, and PIDA max from 4.639 to 0.543 dB. Two related corrections: Py-ART's field metadata carries valid_min/valid_max on specific_attenuation but none on specific_differential_attenuation, so AH above the cap was silently masked on read while the ADP derived from that same AH survived, leaving the pair disagreeing by over a thousand gates. Set both, deriving the ADP bound from the AH bound through the same ADP = c * AH ** d relation. calculate_attenuation_zphi defaults to temp_ref='temperature', so the iso0_field argument was silently ignored and the melting layer came from a sounding_temperature copy whose masked gates had been filled with +9999 degC, marking the entire ice region as liquid. Pass temp_ref explicitly and drop the fill. This changes the melting layer mask on 26 percent of gates volume wide but no output gates, since those gates are already excluded by phase_proc_gates. Co-Authored-By: Claude Opus 5 (1M context) --- cmac/cmac_processing.py | 68 +++++++++++++++++++++++++++++++++++++++++ cmac/cmac_radar.py | 64 +++++++++++++++++++++++++++++++++----- cmac/default_config.py | 3 ++ configs/bnfcsapr2.yaml | 3 ++ 4 files changed, 130 insertions(+), 8 deletions(-) diff --git a/cmac/cmac_processing.py b/cmac/cmac_processing.py index d994272..5456961 100644 --- a/cmac/cmac_processing.py +++ b/cmac/cmac_processing.py @@ -11,6 +11,7 @@ import netCDF4 import numpy as np import pyart +from pyart.correct.phase_proc import det_sys_phase, det_sys_phase_gf from scipy import integrate from scipy import ndimage, interpolate import skfuzzy as fuzz @@ -335,6 +336,73 @@ def get_melt(radar, melt_cat=None, fzl_ceiling=5000.0, fzl = radar.gate_altitude['data'].min() return fzl +def get_sys_phase(radar, gatefilter=None, phidp_field='differential_phase', + ncp_field='normalized_coherent_power', + rhv_field='copol_correlation_coeff', first_gate=30, + n_gates=25, min_gates=10, min_rays=10): + """ + Estimate the PhiDP system phase offset in degrees. + + Takes the median PhiDP over the first n_gates meteorological gates of + each ray, then the median of those across rays. The median is used in + preference to the per-ray minimum that Py-ART's estimators take because + PhiDP here is noisy enough (gate to gate standard deviation of order 30 + degrees) that a minimum sits well below the true offset. + + The Py-ART estimators are kept as fallbacks. Both walk the rays of the + first sweep, so they raise IndexError on volumes whose + sweep_end_ray_index disagrees with the actual ray count, which is the + case for some of the BNF RHI files. + + Returns None if no estimate can be made, in which case the caller should + leave PhiDP alone. + + """ + phidp = radar.fields[phidp_field]['data'][:, first_gate:] + if gatefilter is not None: + included = gatefilter.gate_included[:, first_gate:] + else: + included = ~np.ma.getmaskarray(phidp) + + values = np.ma.filled(phidp, np.nan) + phases = [] + for ray in range(values.shape[0]): + good = values[ray][included[ray] & np.isfinite(values[ray])] + if good.size >= min_gates: + phases.append(np.median(good[:n_gates])) + if len(phases) >= min_rays: + return float(np.median(phases)) + + for estimator in ( + lambda: det_sys_phase_gf(radar, gatefilter, + phidp_field=phidp_field, + first_gate=int(first_gate)), + lambda: det_sys_phase(radar, ncp_field=ncp_field, + rhv_field=rhv_field, + phidp_field=phidp_field)): + try: + sys_phase = estimator() + except (KeyError, IndexError, ValueError, TypeError): + continue + if sys_phase is not None and np.isfinite(sys_phase): + return float(sys_phase) + return None + + +def remove_sys_phase(phidp, sys_phase, wrap_min=-90.0): + """ + Subtract the system phase offset from PhiDP and undo the 0-360 fold. + + PhiDP is measured modulo 360, so subtracting the offset and taking the + remainder recovers the propagation phase in a single step. The result is + wrapped into [wrap_min, wrap_min + 360) rather than [0, 360) so that + noise about zero at close range stays slightly negative instead of + folding up to nearly 360 degrees. + + """ + return np.mod(phidp - sys_phase - wrap_min, 360.0) + wrap_min + + def fix_phase_fields(orig_kdp, orig_phidp, rrange, happy_kdp, max_kdp=15.0): diff --git a/cmac/cmac_radar.py b/cmac/cmac_radar.py index 7886a43..127aa14 100644 --- a/cmac/cmac_radar.py +++ b/cmac/cmac_radar.py @@ -5,6 +5,7 @@ import copy import json import sys +import warnings import xarray as xr import numpy as np @@ -13,7 +14,7 @@ from .cmac_processing import ( do_my_fuzz, get_melt, get_texture, fix_phase_fields, gen_clutter_field_from_refl, beam_block, - snow_rate, rain_rate) + snow_rate, rain_rate, get_sys_phase, remove_sys_phase) from .config import (get_cmac_values, get_field_names, get_metadata, get_zs_relationships, get_default_metadata) from . import csu_kdp @@ -336,6 +337,28 @@ def cmac(radar, sonde, config, geotiff=None, flip_velocity=False, fill_value = -9999. dp = radar.fields[field_config['input_phidp_field']]['data'] dz = radar.fields[field_config['reflectivity']]['data'] + # calc_kdp_bringi assumes PhiDP has the system phase offset removed + # and is already unfolded. Feeding it the raw phase carries the + # offset (~220 degrees at BNF) through to the Z-PHI attenuation, + # where it is read as propagation phase and inflates the specific + # attenuation by one to two orders of magnitude. + sys_phase = cmac_config.get('phidp_sys_phase', None) + if sys_phase is None: + sys_phase = get_sys_phase( + radar, gatefilter=kdp_gates, + phidp_field=field_config['input_phidp_field'], + ncp_field=field_config['normalized_coherent_power'], + rhv_field=field_config['cross_correlation_ratio']) + if sys_phase is None: + warnings.warn( + 'Unable to estimate the PhiDP system phase offset; the raw ' + 'phase will be processed unchanged. Set phidp_sys_phase in ' + 'the configuration to supply it manually.', UserWarning) + else: + dp = remove_sys_phase(dp, sys_phase) + if verbose: + print('## PhiDP system phase offset removed: ' + '%.2f degrees' % sys_phase) dp = dp.filled(fill_value) dz = dz.filled(fill_value) rng = np.tile(radar.range['data'], (radar.nrays, 1)) / 1e3 @@ -349,6 +372,8 @@ def cmac(radar, sonde, config, geotiff=None, flip_velocity=False, kdp = pyart.config.get_metadata("corrected_specific_differential_phase") phidp["data"] = phidp_data kdp["data"] = kdp_data + if sys_phase is not None: + phidp["system_phase_offset"] = np.round(sys_phase, 4) print("Processed phase") @@ -417,22 +442,32 @@ def cmac(radar, sonde, config, geotiff=None, flip_velocity=False, radar, 'differential_phase', gatefilter=phase_proc_gates, size=cmac_config.get('phidp_despeckle_size', 49)) - radar.fields['sounding_temperature_filled'] = copy.deepcopy(radar.fields['sounding_temperature']) - radar.fields['sounding_temperature_filled']['data'] = np.where( - radar.fields['sounding_temperature_filled']['data'] > -100., - radar.fields['sounding_temperature_filled']['data'], 9999.) + # calculate_attenuation_zphi defaults to temp_ref='temperature', which + # ignores iso0_field altogether. Ask for the iso0 reference explicitly so + # that the height_over_iso0 field built above is what sets the melting + # layer, and fall back to the gate sounding temperature if the 0 degC + # level could not be located. + if np.isfinite(np.ma.filled(iso0, np.nan)): + temp_ref = 'height_over_iso0' + temp_field = None + else: + warnings.warn( + 'Could not locate the 0 degC level in the sounding; falling back ' + 'to the gate sounding temperature to set the melting layer.', + UserWarning) + temp_ref = 'temperature' + temp_field = 'sounding_temperature' (spec_at, pia_dict, cor_z, spec_diff_at, pida_dict, cor_zdr) = pyart.correct.calculate_attenuation_zphi( - radar, temp_field='sounding_temperature_filled', + radar, temp_field=temp_field, zdr_field=field_config['zdr_field'], pia_field=field_config['pia_field'], iso0_field='height_over_iso0', phidp_field=field_config['phidp_field'], refl_field=field_config['refl_field'], c=c_coef, d=d_coef, a_coef=attenuation_a_coef, beta=beta_coef, - gatefilter=phase_proc_gates) + gatefilter=phase_proc_gates, temp_ref=temp_ref) # cor_zdr['data'] += cmac_config['zdr_offset'] Now taken care of at start - del radar.fields['sounding_temperature_filled'] radar.add_field('specific_attenuation', spec_at, replace_existing=True) radar.add_field('path_integrated_attenuation', pia_dict, replace_existing=True) @@ -443,6 +478,19 @@ def cmac(radar, sonde, config, geotiff=None, flip_velocity=False, replace_existing=True) radar.add_field('corrected_differential_reflectivity', cor_zdr, replace_existing=True) + + # Py-ART's field metadata puts valid_min/valid_max on specific_attenuation + # but none on specific_differential_attenuation, so AH above the cap is + # silently masked on read while the ADP derived from that same AH is not. + # Set both, deriving the ADP bound from the AH bound through the + # ADP = c * AH ** d relation used to compute it, so the pair always masks + # the same gates. + ah_valid_max = float(cmac_config.get('specific_attenuation_valid_max', 1.0)) + radar.fields['specific_attenuation']['valid_min'] = 0.0 + radar.fields['specific_attenuation']['valid_max'] = ah_valid_max + radar.fields['specific_differential_attenuation']['valid_min'] = 0.0 + radar.fields['specific_differential_attenuation']['valid_max'] = float( + np.round(c_coef * ah_valid_max ** d_coef, 4)) radar.fields['corrected_velocity']['units'] = 'm/s' if 'valid_min' not in radar.fields['corrected_velocity'].keys(): diff --git a/cmac/default_config.py b/cmac/default_config.py index 59585ca..7969f64 100644 --- a/cmac/default_config.py +++ b/cmac/default_config.py @@ -1019,6 +1019,9 @@ 'rain_rate_a_coef_Kdp': 25.1, 'rain_rate_b_coef_Kdp': 0.777, 'kdp_method': "bringi", + # PhiDP system phase offset in degrees. None estimates it from + # each volume; set a number to pin it. + 'phidp_sys_phase': None, 'beam_width': 1.0, 'radar_height_offset': 10.0,}, # We expect clutter corrected fields now diff --git a/configs/bnfcsapr2.yaml b/configs/bnfcsapr2.yaml index 90cc023..2325863 100644 --- a/configs/bnfcsapr2.yaml +++ b/configs/bnfcsapr2.yaml @@ -99,6 +99,9 @@ cmac_values: rain_rate_a_coef_Kdp: 25.1 rain_rate_b_coef_Kdp: 0.777 kdp_method: bringi + # PhiDP system phase offset in degrees. Leave as null to estimate it + # from each volume; set a number to pin it. + phidp_sys_phase: null beam_width: 1.0 radar_height_offset: 10.0 # Fuzzy-logic membership functions. Each class lists every field used