From 7756f9ce78d41bbb05a4769e9f1c49ac3f728268 Mon Sep 17 00:00:00 2001 From: Russell Desiderio Date: Sat, 13 May 2017 10:43:37 -0700 Subject: [PATCH] botsflu_2017-05-12 --- ion_functions/data/prs_functions.py | 229 +++++++++++------- ion_functions/data/test/test_prs_functions.py | 150 +++++++++++- 2 files changed, 287 insertions(+), 92 deletions(-) diff --git a/ion_functions/data/prs_functions.py b/ion_functions/data/prs_functions.py index 74283bf..007cba6 100755 --- a/ion_functions/data/prs_functions.py +++ b/ion_functions/data/prs_functions.py @@ -42,9 +42,9 @@ BOTSFLU: anchor_bin - calc_daydepth_plus calc_meandepth_plus calculate_sliding_means + calculate_sliding_slopes Functions calculating event notifications; they return either True or False. @@ -290,7 +290,7 @@ def prs_botsflu_time15s(timestamp, botpres): OOI (2015). Data Product Specification for Seafloor Uplift and Subsidence (BOTSFLU) from the BOTPT instrument. Document Control Number 1341-00080. """ - ###### the second calling argument is a placeholder + ###### [in development]: the second calling argument is a placeholder time15s, _, _ = anchor_bin_rawdata_to_15s(timestamp, botpres) return time15s @@ -591,12 +591,12 @@ def prs_botsflu_time24h(time15s): (BOTSFLU) from the BOTPT instrument. Document Control Number 1341-00080. """ # the second calling argument is a placeholder - time24h, _ = anchor_bin_detided_to_24h(time15s, None) + time24h, _ = anchor_bin_detided_to_24h(time15s, None, None) return time24h -def prs_botsflu_daydepth(timestamp, botpres): +def prs_botsflu_daydepth(timestamp, botpres, dday_coverage=0.90): """ Description: @@ -607,8 +607,8 @@ def prs_botsflu_daydepth(timestamp, botpres): 2015-01-14: Russell Desiderio. Initial code. 2017-05-05: Russell Desiderio. Changed time24h time base to span the entire - dataset including data gaps. Therefore added - lines to insert Nans at the data gaps. + dataset including data gaps. + Usage @@ -629,12 +629,23 @@ def prs_botsflu_daydepth(timestamp, botpres): OOI (2015). Data Product Specification for Seafloor Uplift and Subsidence (BOTSFLU) from the BOTPT instrument. Document Control Number 1341-00080. """ - daydepth = calc_daydepth_plus(timestamp, botpres) + # calculate 15sec bin timestamps and de-tided depth. + time15s, meandepth, _ = calc_meandepth_plus(timestamp, botpres) + + # bin the 15sec data into 24 hour bins so that the timestamps are at midnight. + # to calculate daydepth, don't need the time24h timestamps. + _, daydepth = anchor_bin_detided_to_24h(time15s, meandepth, dday_coverage) + + # downstream data products no longer require the mask_nonzero variable return daydepth + #daydepth = calc_daydepth_plus(timestamp, botpres) + # + #return daydepth + -def prs_botsflu_4wkrate(timestamp, botpres): +def prs_botsflu_4wkrate(timestamp, botpres, dday_coverage=0.9, rate_coverage=0.75): """ Description: @@ -669,11 +680,11 @@ def prs_botsflu_4wkrate(timestamp, botpres): (BOTSFLU) from the BOTPT instrument. Document Control Number 1341-00080. """ # calculate daydepth - data_w_gaps = calc_daydepth_plus(timestamp, botpres) + daydepth = prs_botsflu_daydepth(timestamp, botpres, dday_coverage) # 4 weeks of data window_size = 29 - botsflu_4wkrate = calculate_sliding_slopes(data_w_gaps, window_size) + botsflu_4wkrate = calculate_sliding_slopes(daydepth, window_size, rate_coverage) # convert units: # the units of the slopes are [y]/[x] = meters/day; # to get units of cm/yr, multiply by 100cm/m * 365 days/yr @@ -682,7 +693,7 @@ def prs_botsflu_4wkrate(timestamp, botpres): return botsflu_4wkrate -def prs_botsflu_8wkrate(timestamp, botpres): +def prs_botsflu_8wkrate(timestamp, botpres, dday_coverage=0.9, rate_coverage=0.75): """ Description: @@ -716,12 +727,12 @@ def prs_botsflu_8wkrate(timestamp, botpres): OOI (2015). Data Product Specification for Seafloor Uplift and Subsidence (BOTSFLU) from the BOTPT instrument. Document Control Number 1341-00080. """ - # calculate daydepth and the mask of nonzero data bins. - data_w_gaps = calc_daydepth_plus(timestamp, botpres) + # calculate daydepth + daydepth = prs_botsflu_daydepth(timestamp, botpres, dday_coverage) # 8 weeks of data window_size = 57 - botsflu_8wkrate = calculate_sliding_slopes(data_w_gaps, window_size) + botsflu_8wkrate = calculate_sliding_slopes(daydepth, window_size, rate_coverage) # convert units: # the units of the slopes are [y]/[x] = meters/day; # to get units of cm/yr, multiply by 100cm/m * 365 days/yr @@ -843,7 +854,7 @@ def anchor_bin_rawdata_to_15s(time, data): return bin_timestamps, binned_data, mask_nonzero -def anchor_bin_detided_to_24h(time, data): +def anchor_bin_detided_to_24h(time, data, dday_coverage): """ Description: @@ -925,6 +936,7 @@ def anchor_bin_detided_to_24h(time, data): """ bin_duration = 86400.0 # number of seconds in a day half_bin = bin_duration/2.0 + max_count = 86400.0/15.0 # maximum number of values in a day's bin # anchor time-centered bins by determining the start time to be half a bin # before the first 'anchor timestamp', which will an integral number of @@ -936,8 +948,11 @@ def anchor_bin_detided_to_24h(time, data): bin_number = np.floor(time_elapsed).astype(int) # the number of elements in each bin is given by bin_count = np.bincount(bin_number).astype(float) - # create a logical mask of non-zero bin_count values - mask_nonzero = (bin_count != 0) + # bin_count is used as a divisor to calculate mean values at each bin + # bins with bincounts below the threshold value will have a nan value + bin_count[bin_count/max_count < dday_coverage] = np.nan + # replace 0 with nan + bin_count[bin_count == 0] = np.nan # directly calculate bin timestamp, units of [sec]: # the midpoint of the data interval is used. @@ -945,58 +960,12 @@ def anchor_bin_detided_to_24h(time, data): # sum the values in each time bin, and put into the variable binned_data binned_data = np.bincount(bin_number, data) - # divide the values in non-empty bins by the number of values in each bin - binned_data = binned_data[mask_nonzero]/bin_count[mask_nonzero] - - # re-constitute the original data, with data gaps (if present) represented by nans. - daydepth = np.zeros(mask_nonzero.size) + np.nan - daydepth[mask_nonzero] = binned_data + # divide the values in each bin by the number of values in each bin + daydepth = binned_data/bin_count return bin_timestamps, daydepth -def calc_daydepth_plus(timestamp, botpres): - """ - Description: - - Worker function to calculate the botsflu data product daydepth plus an - additional boolean mask required to calculate other botsflu data products - downstream from daydepth. - - Implemented by: - - 2015-01-14: Russell Desiderio. Initial code. - - Usage - - daydepth, mask_nonzero = calc_daydepth_plus(timestamp, botpres) - - where - - daydepth = BOTSFLU-DAYDEPTH_L2 [m] - mask_nonzero = boolean of positions of non-empty 24 hr bins - timestamp = OOI system timestamps [sec since 01-01-1900] - botpres = BOTPRES_L1 [psia] - - Notes: - - References: - - OOI (2015). Data Product Specification for Seafloor Uplift and Subsidence - (BOTSFLU) from the BOTPT instrument. Document Control Number 1341-00080. - """ - # calculate 15sec bin timestamps and de-tided depth. - time15s, meandepth, _ = calc_meandepth_plus(timestamp, botpres) - - # bin the 15sec data into 24 hour bins so that the timestamps are at midnight. - # to calculate daydepth, don't need the time24h timestamps. - - _, daydepth = anchor_bin_detided_to_24h(time15s, meandepth) - - # downstream data products no longer require the mask_nonzero variable - return daydepth - - def calc_meandepth_plus(timestamp, botpres): """ Description: @@ -1068,6 +1037,7 @@ def calculate_sliding_means(data, window_size): Implemented by: 2015-01-13: Russell Desiderio. Initial code. + 2017-05-06: Russell Desiderio. Updated to make work with odd window_size. Usage @@ -1077,24 +1047,29 @@ def calculate_sliding_means(data, window_size): means = 1D array of sliding means data = 1D array of data - window_size = even integer + window_size = window size, integer data type + + Notes + + The unit test values were calculated in Matlab, so that the python convolution + result for even sized windows is shifted by 1 element to match the matlab result. + """ kk = np.ones(window_size) / window_size means = np.convolve(data, kk, 'same') - # matlab and numpy behave differently for even window sizes, so - means = np.roll(means, -1) - # in this application, window_size is always even. - hfwin = window_size/2 - 1 # nan out data with boundary effects at edges. - # the values in the array means will be of type float because kk is a float, - # so that the np.nan assignment statements will work as intended. - means[0:hfwin] = np.nan - means[-hfwin-1:] = np.nan + # integer arithmetic will 'truncate' 5/2 to 2 and -5/2 to -3. + means[:window_size/2] = np.nan + means[-((window_size-1)/2):] = np.nan + + # matlab and numpy behave differently for even window sizes, so roll the python + # result to mimic matlab + means = np.roll(means, -np.mod(window_size+1, 2)) # roll only if window_size is even return means -def calculate_sliding_slopes(data, window_size): +def calculate_sliding_slopes(data, window_size, coverage_threshold): """ Description: @@ -1107,6 +1082,9 @@ def calculate_sliding_slopes(data, window_size): 2017-05-03: Russell Desiderio. Initial code. Replaces the much faster Moore-Penrose pseudo-inverse method so that nan-masking can be incorporated. + 2017-05-08: Russell Desiderio. Added the 70% coverage criterion: if the number of non-Nan + data points in a window is greater than or equal to 70% of + the maximum possible window points calculate the slope. Usage @@ -1117,35 +1095,68 @@ def calculate_sliding_slopes(data, window_size): slopes = 1D array of sliding slopes data = 1D array of data window_size = integer + coverage = fractional window fill threshold for calculation of slope values Notes - Used equations 14.2.15-14.2.17 in the Numerical Recipes reference below. + The robust regression equations are taken from equations 14.2.15-14.2.17 in the Numerical + Recipes reference below. - The slopes are backwards-looking, not centered. The first non-nan value occurs - at (python) index window_size, and is the slope of a regression of the first - window_size points. + Before the May 2017 modifications, just one Nan within a window would result in a Nan value + for the slope. The routine now will calculate slopes by ignoring Nans, and if the coverage + is 70% or better, a calculated value will result. If the coverage is less than 70%, then the + output value is Nan. + + The data vector is padded so that bins within a window of the beginning of the data record + that satisfy the coverage criterion will have non-Nan data product values. References: Press, Flannery, Teukolsky and Vetterling. Numerical Recipes, 1986; 1987 reprint. Cambridge University Press. page 507. """ - # odd window sizes are expected; if even, increment by one + # ODD WINDOW SIZES are expected; if even, increment by one window_size = window_size + 1 - np.mod(window_size, 2) half_window = window_size / 2 # this will 'floor' when window_size is odd as desired - slopes = np.zeros(data.size) + np.nan - x = np.arange(window_size).astype('float') - for ii in range(half_window, data.size - half_window): - y = data[(ii-half_window):(ii+half_window+1)] - #x = x[~np.isnan(y)] # LATER! - #y = y[~np.isnan(y)] # LATER! + # pad front end of data with Nans (because for loop is backwards-looking) + npts = 2 * half_window + data.size + padded_data = np.zeros(npts) + np.nan + padded_data[window_size-1:] = data + + # first calculate values for all sliding windows + slopes = np.zeros(npts) + np.nan + abscissa = np.arange(npts).astype('float') + for ii in range(window_size-1, npts): + x = abscissa[(ii-2*half_window):(ii+1)] # rather than np.arange-ing in each iteration + y = padded_data[(ii-2*half_window):(ii+1)] + x = x[~np.isnan(y)] + y = y[~np.isnan(y)] + if x.size < 2: + continue # trap out cases of either 0 or 1 valid datapoint in window tti = x - x.mean(axis=0) stt = np.sum(tti * tti) slopes[ii] = np.sum(tti * y) / stt - - slopes = np.roll(slopes, half_window) + # get rid of padding + slopes = slopes[2*half_window:] + + # now determine the fraction of good values per window (vectorized): + # first change non-nan values to 1, then nan values to 0, and take the average + padded_data[~np.isnan(padded_data)] = 1.0 + padded_data[np.isnan(padded_data)] = 0.0 + fraction_good = calculate_sliding_means(padded_data, window_size) + # convert to backwards-looking + fraction_good = np.roll(fraction_good, half_window) + # get rid of padding + fraction_good = fraction_good[2*half_window:] + + # nan out fractional values (means) in sliding windows less than coverage threshold + # there can be issues involving roundoff error when considering 100% coverage, so: + machine_epsilon = np.finfo(np.float).eps + fraction_good = fraction_good + 100 * machine_epsilon + # avoid a python warning message by trapping out nans in the conditional + fraction_good[np.isnan(fraction_good)] = -999.0 # any negative number will work as intended + slopes[fraction_good < coverage_threshold] = np.nan return slopes @@ -1372,6 +1383,48 @@ def anchor_bin_DEPRECATED(time, data, bin_duration, mode): return bin_timestamps, binned_data, mask_nonzero +def calc_daydepth_plus_deprecated(timestamp, botpres): + """ + Description: + + Worker function to calculate the botsflu data product daydepth plus an + additional boolean mask required to calculate other botsflu data products + downstream from daydepth. + + Implemented by: + + 2015-01-14: Russell Desiderio. Initial code. + + Usage + + daydepth, mask_nonzero = calc_daydepth_plus(timestamp, botpres) + + where + + daydepth = BOTSFLU-DAYDEPTH_L2 [m] + mask_nonzero = boolean of positions of non-empty 24 hr bins + timestamp = OOI system timestamps [sec since 01-01-1900] + botpres = BOTPRES_L1 [psia] + + Notes: + + References: + + OOI (2015). Data Product Specification for Seafloor Uplift and Subsidence + (BOTSFLU) from the BOTPT instrument. Document Control Number 1341-00080. + """ + # calculate 15sec bin timestamps and de-tided depth. + time15s, meandepth, _ = calc_meandepth_plus(timestamp, botpres) + + # bin the 15sec data into 24 hour bins so that the timestamps are at midnight. + # to calculate daydepth, don't need the time24h timestamps. + + _, daydepth = anchor_bin_detided_to_24h(time15s, meandepth) + + # downstream data products no longer require the mask_nonzero variable + return daydepth + + def calculate_sliding_slopes__MoorePenrose(data, window_size): # DEPRECATED because nan-masking cannot be implemented with this algorithm # """ diff --git a/ion_functions/data/test/test_prs_functions.py b/ion_functions/data/test/test_prs_functions.py index efbf764..7677688 100755 --- a/ion_functions/data/test/test_prs_functions.py +++ b/ion_functions/data/test/test_prs_functions.py @@ -286,7 +286,8 @@ def test_calculate_sliding_means(self): even window sizes; specifically, make sure edge effects are nan'd out appropriately. - Implemented by Russell Desiderio, January 2015. + Russell Desiderio, January 2015. Original code. + May 2017: Added case of odd window-size for use in sliding_slope with Nanmasking """ xpctd = np.array([np.nan, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, np.nan, np.nan]) data = np.arange(10.0) @@ -294,6 +295,12 @@ def test_calculate_sliding_means(self): calc = prsfunc.calculate_sliding_means(data, window_size) np.testing.assert_array_equal(calc, xpctd) + xpctd = np.array([np.nan, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, np.nan]) + data = np.arange(10.0) + window_size = 3 + calc = prsfunc.calculate_sliding_means(data, window_size) + np.testing.assert_array_equal(calc, xpctd) + def test_calculate_sliding_slopes(self): """ Test the calculation of backwards-looking rolling slopes. @@ -306,8 +313,37 @@ def test_calculate_sliding_slopes(self): xpctd = np.array([np.nan, np.nan, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) data = np.arange(10.0) window_size = 3 - calc = prsfunc.calculate_sliding_slopes(data, window_size) + rate_coverage = 1.0 + calc = prsfunc.calculate_sliding_slopes(data, window_size, rate_coverage) np.testing.assert_allclose(calc, xpctd, rtol=0.0, atol=1.e-12) + """ + Additional test, RDesiderio, May 12, 2017 + Verify nan placement after routine has been altered to accept a 75% fill criterion. + """ + # set up test data with slope=2 + data = 2 * np.arange(150.0) + # nan out middle section + data[50:100] = np.nan + + # when the sliding slope routine is applied, the result will be sets of contiguous nans + # and non-Nan values whose indices depend on the window size and on the location of nans + # in the test data. + # 29 point window (4weekrate case) + xpctd_w29 = np.ones(150) + 1.0 + xpctd_w29[0:21] = np.nan + xpctd_w29[57:121] = np.nan + window_size = 29 + rate_coverage = 0.75 + calc = prsfunc.calculate_sliding_slopes(data, window_size, rate_coverage) + np.testing.assert_array_equal(calc, xpctd_w29) + # 57 point window (8weekrate case) + xpctd_w57 = np.ones(150) + 1.0 + xpctd_w57[0:42] = np.nan + xpctd_w57[64:142] = np.nan + window_size = 57 + rate_coverage = 0.75 + calc = prsfunc.calculate_sliding_slopes(data, window_size, rate_coverage) + np.testing.assert_array_equal(calc, xpctd_w57) def test_prs_botsflu_auxiliary_timestamps(self): """ @@ -479,11 +515,33 @@ def test_prs_botsflu_daydepth(self): Test the calculation of DAYDEPTH. Implemented by Russell Desiderio, January 2015. + RDesiderio. May 12, 2015. Added 0.9 and default coverage test """ + # case(1) January 2015 unit test values xpctd_daydepth = np.array([-1510.8705, -1510.9055, -1510.9119, -1510.8491, -1510.8596, -1510.8354, -1510.8791, -1510.8807, -1510.8378, -1510.8279, -1510.8530, -1512.2859, -1513.2018, -1513.1660, -1513.1128, -1513.0478]) + dday_coverage = 0.0 # if a bin contains a non-Nan value, its datum will exist + daydepth = prsfunc.prs_botsflu_daydepth(self.ss1900, self.botpres, dday_coverage) + calc = daydepth[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_daydepth, rtol=0, atol=1.e-4) + + # case(2) May 2017 unit test values, coverage = 0.90 + xpctd_daydepth = np.array([np.nan, -1510.9055, np.nan, -1510.8491, + -1510.8596, -1510.8354, -1510.8791, -1510.8807, + -1510.8378, -1510.8279, -1510.8530, -1512.2859, + -1513.2018, -1513.1660, -1513.1128, np.nan]) + dday_coverage = 0.9 + daydepth = prsfunc.prs_botsflu_daydepth(self.ss1900, self.botpres, dday_coverage) + calc = daydepth[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_daydepth, rtol=0, atol=1.e-4) + + # case(3) May 2017 default, no coverage specified in argument list + xpctd_daydepth = np.array([np.nan, -1510.9055, np.nan, -1510.8491, + -1510.8596, -1510.8354, -1510.8791, -1510.8807, + -1510.8378, -1510.8279, -1510.8530, -1512.2859, + -1513.2018, -1513.1660, -1513.1128, np.nan]) daydepth = prsfunc.prs_botsflu_daydepth(self.ss1900, self.botpres) calc = daydepth[self.idx_24h] np.testing.assert_allclose(calc, xpctd_daydepth, rtol=0, atol=1.e-4) @@ -493,11 +551,53 @@ def test_prs_botsflu_4wkrate(self): Test the calculation of 4WKRATE. Implemented by Russell Desiderio, January 2015. - """ + RDesiderio. May 12, 2015. 3 additional tests for coverage parameters; 4 total + (1) dday_coverage = 0.0; rate_coverage = 1.0 (January 2015 unit test values) + (2) dday_coverage = 0.0; rate_coverage = 0.75 + (3) dday_coverage = 0.9; rate_coverage = 0.75 + (4) default (no argument) default (no argument) (should be same as (3)) + """ + # case (1) + dday_coverage = 0.0 + rate_coverage = 1.0 xpctd_4wkrate = np.array([np.nan, np.nan, np.nan, np.nan, np.nan, 70.5000, 28.8000, -23.6372, 7.8316, 38.7091, 57.6564, -407.2349, -3024.1664, -4291.1577, -3982.1142, -2118.0963]) + b_4wkrate = prsfunc.prs_botsflu_4wkrate( + self.ss1900, self.botpres, dday_coverage, rate_coverage) + calc = b_4wkrate[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_4wkrate, rtol=0, atol=1.e-4) + + # case (2) + dday_coverage = 0.0 + rate_coverage = 0.75 + xpctd_4wkrate = np.array([np.nan, np.nan, np.nan, np.nan, + 64.1039, 70.5000, 28.8000, -23.6372, + 7.8316, 38.7091, 57.6564, -407.2349, + -3024.1664, -4291.1577, -3982.1142, -2118.0963]) + b_4wkrate = prsfunc.prs_botsflu_4wkrate( + self.ss1900, self.botpres, dday_coverage, rate_coverage) + calc = b_4wkrate[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_4wkrate, rtol=0, atol=1.e-4) + + # case (3) + dday_coverage = 0.9 + rate_coverage = 0.75 + xpctd_4wkrate = np.array([np.nan, np.nan, np.nan, np.nan, + np.nan, 63.2291, 0.4971, -23.6372, + 7.8316, 38.7091, 57.6564, -407.2349, + -3024.1664, -4291.1577, -3982.1142, -2276.7628]) + b_4wkrate = prsfunc.prs_botsflu_4wkrate( + self.ss1900, self.botpres, dday_coverage, rate_coverage) + calc = b_4wkrate[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_4wkrate, rtol=0, atol=1.e-4) + + # case (4) + xpctd_4wkrate = np.array([np.nan, np.nan, np.nan, np.nan, + np.nan, 63.2291, 0.4971, -23.6372, + 7.8316, 38.7091, 57.6564, -407.2349, + -3024.1664, -4291.1577, -3982.1142, -2276.7628]) b_4wkrate = prsfunc.prs_botsflu_4wkrate(self.ss1900, self.botpres) calc = b_4wkrate[self.idx_24h] np.testing.assert_allclose(calc, xpctd_4wkrate, rtol=0, atol=1.e-4) @@ -507,11 +607,53 @@ def test_prs_botsflu_8wkrate(self): Test the calculation of 8WKRATE. Implemented by Russell Desiderio, January 2015. - """ + RDesiderio. May 12, 2015. 3 additional tests for coverage parameters; 4 total + (1) dday_coverage = 0.0; rate_coverage = 1.0 (January 2015 unit test values) + (2) dday_coverage = 0.0; rate_coverage = 0.75 + (3) dday_coverage = 0.9; rate_coverage = 0.75 + (4) default (no argument) default (no argument) (should be same as (3)) + """ + # case (1) + dday_coverage = 0.0 + rate_coverage = 1.0 xpctd_8wkrate = np.array([np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, 35.7686, -90.1254, -882.7773, -1506.2385, -1921.8774, -2120.7657]) + b_8wkrate = prsfunc.prs_botsflu_8wkrate( + self.ss1900, self.botpres, dday_coverage, rate_coverage) + calc = b_8wkrate[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_8wkrate, rtol=0, atol=1.e-4) + + # case (2) + dday_coverage = 0.0 + rate_coverage = 0.75 + xpctd_8wkrate = np.array([np.nan, np.nan, np.nan, np.nan, + np.nan, np.nan, np.nan, 24.5993, + 28.7563, 32.6131, 35.7686, -90.1254, + -882.7773, -1506.2385, -1921.8774, -2120.7657]) + b_8wkrate = prsfunc.prs_botsflu_8wkrate( + self.ss1900, self.botpres, dday_coverage, rate_coverage) + calc = b_8wkrate[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_8wkrate, rtol=0, atol=1.e-4) + + # case (3) + dday_coverage = 0.9 + rate_coverage = 0.75 + xpctd_8wkrate = np.array([np.nan, np.nan, np.nan, np.nan, + np.nan, np.nan, np.nan, np.nan, + 26.1114, 30.4322, 31.6141, -113.8871, + -882.7773, -1506.2385, -1921.8774, -2150.3933]) + b_8wkrate = prsfunc.prs_botsflu_8wkrate( + self.ss1900, self.botpres, dday_coverage, rate_coverage) + calc = b_8wkrate[self.idx_24h] + np.testing.assert_allclose(calc, xpctd_8wkrate, rtol=0, atol=1.e-4) + + # case (4) + xpctd_8wkrate = np.array([np.nan, np.nan, np.nan, np.nan, + np.nan, np.nan, np.nan, np.nan, + 26.1114, 30.4322, 31.6141, -113.8871, + -882.7773, -1506.2385, -1921.8774, -2150.3933]) b_8wkrate = prsfunc.prs_botsflu_8wkrate(self.ss1900, self.botpres) calc = b_8wkrate[self.idx_24h] np.testing.assert_allclose(calc, xpctd_8wkrate, rtol=0, atol=1.e-4)