From d4603c54e224c2a5db953a4c5800cc092be2ee8a Mon Sep 17 00:00:00 2001 From: Steve Morley Date: Thu, 2 Apr 2020 18:05:37 -0600 Subject: [PATCH] ensure noisy gap filling is correctly applied - propagated fillval from the downloaded CDF to fill_gaps - made sure fill_gaps' output data wasn't overwriting the input --- advect_imf.py | 10 ++++++++-- missing.py | 5 +++-- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/advect_imf.py b/advect_imf.py index d4fdae0..e4b850d 100644 --- a/advect_imf.py +++ b/advect_imf.py @@ -48,7 +48,10 @@ def load_acedata(tstart, tend, noise=True, proxy=None): for dataset, local_name, cdaweb_name in [(mag_data, 'b', 'BGSM'), (swepam_data, 'u', 'V_GSM'), (swepam_data, '', 'SC_pos_GSM')]: - t, values = dataset['Epoch'], fill_gaps(dataset[cdaweb_name][:, i], noise=noise) + t = dataset['Epoch'] + values = fill_gaps(dataset[cdaweb_name][:, i], + fillval=dataset[cdaweb_name].attrs['FILLVAL'], + noise=noise) # Grab the appropriate component from # VALIDMIN and VALIDMAX attributes @@ -99,7 +102,10 @@ def load_dscovr(tstart, tend, noise=True, proxy=None): (mag_data, 'b', 'B1GSE', 'Epoch1'), (plasma_data, 'u', 'V_GSE', 'Epoch'), (orbit_data, '', 'GSE_POS', 'Epoch')]: - t, values = dataset[cdaweb_time_var], fill_gaps(dataset[cdaweb_name][:, i], noise=noise) + t = dataset[cdaweb_time_var] + values = fill_gaps(dataset[cdaweb_name][:, i], + fillval=dataset[cdaweb_name].attrs['FILLVAL'], + noise=noise) for attr in ('VALIDMIN', 'VALIDMAX'): diff --git a/missing.py b/missing.py index 84e5f5a..f61c678 100644 --- a/missing.py +++ b/missing.py @@ -60,11 +60,12 @@ def fill_gaps(data, fillval=9999999, sigma=5, winsor=0.05, noise=False, constrai # draw fluctuations from CDF and apply to linearly filled gaps for gap in gaps: for i in range(gap[1]-gap[0]+1): - data[gap[0]+i] += dx[p.searchsorted(random.random())] + series[gap[0]+i] += dx[p.searchsorted(random.random())] # cap variable if it should be strictly positive (e.g. number density) # use lowest measured value as floor if constrain and series.min() > 0.0: - data[data < series.min()] = series.min() + series[series < series.min()] = series.min() + return series return data