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