Merge pull request #195 from NREL/gb/precip_skill

Gb/precip skill

sup3r/bias/bias_calc.py

@@ -44,7 +44,8 @@ def __init__(self,
                  bias_handler='DataHandlerNCforCC',
                  base_handler_kwargs=None,
                  bias_handler_kwargs=None,
-                 decimals=None):
+                 decimals=None,
+                 match_zero_rate=False):
         """
         Parameters
         ----------
@@ -94,6 +95,15 @@ class to be retrieved from the rex/sup3r library. If a
             decimals, this gets passed to np.around(). If decimals
             is negative, it specifies the number of positions to
             the left of the decimal point.
+        match_zero_rate : bool
+            Option to fix the frequency of zero values in the biased data. The
+            lowest percentile of values in the biased data will be set to zero
+            to match the percentile of zeros in the base data. If
+            SkillAssessment is being run and this is True, the distributions
+            will not be mean-centered. This helps resolve the issue where
+            global climate models produce too many days with small
+            precipitation totals e.g., the "drizzle problem".
+            Ref: Polade et al., 2014 https://doi.org/10.1038/srep04364
         """
 
         logger.info('Initializing DataRetrievalBase for base dset "{}" '
@@ -110,6 +120,7 @@ class to be retrieved from the rex/sup3r library. If a
         self.bias_handler_kwargs = bias_handler_kwargs or {}
         self.bad_bias_gids = []
         self._distance_upper_bound = distance_upper_bound
+        self.match_zero_rate = match_zero_rate
 
         if isinstance(self.base_fps, str):
             self.base_fps = sorted(glob(self.base_fps))
@@ -367,6 +378,7 @@ def get_bias_data(self, bias_gid):
         bias_data : np.ndarray
             1D array of temporal data at the requested gid.
         """
+
         idx = np.where(self.bias_gid_raster == bias_gid)
         if self.bias_dh.data is None:
             self.bias_dh.load_cached_data()
@@ -486,6 +498,46 @@ def get_base_data(cls,
 
         return out_data, out_ti
 
+    @staticmethod
+    def _match_zero_rate(bias_data, base_data):
+        """The lowest percentile of values in the biased data will be set to
+        zero to match the percentile of zeros in the base data. This helps
+        resolve the issue where global climate models produce too many days
+        with small precipitation totals e.g., the "drizzle problem".
+        Ref: Polade et al., 2014 https://doi.org/10.1038/srep04364
+
+        Parameters
+        ----------
+        bias_data : np.ndarray
+            1D array of biased data observations.
+        base_data : np.ndarray
+            1D array of base data observations.
+
+        Returns
+        -------
+        bias_data : np.ndarray
+            1D array of biased data observations. Values below the quantile
+            associated with zeros in base_data will be set to zero
+        """
+
+        q_zero_base_in = np.nanmean(base_data == 0)
+        q_zero_bias_in = np.nanmean(bias_data == 0)
+
+        q_bias = np.linspace(0, 1, len(bias_data))
+        min_value_bias = np.interp(q_zero_base_in, q_bias, sorted(bias_data))
+
+        bias_data[bias_data < min_value_bias] = 0
+
+        q_zero_base_out = np.nanmean(base_data == 0)
+        q_zero_bias_out = np.nanmean(bias_data == 0)
+
+        logger.debug('Input bias/base zero rate is {:.3e}/{:.3e}, '
+                     'output is {:.3e}/{:.3e}'
+                     .format(q_zero_bias_in, q_zero_base_in,
+                             q_zero_bias_out, q_zero_base_out))
+
+        return bias_data
+
     @staticmethod
     def _read_base_sup3r_data(dh, base_dset, base_gid):
         """Read baseline data from a sup3r DataHandler
@@ -727,7 +779,8 @@ def _run_single(cls,
                     daily_reduction,
                     bias_ti,
                     decimals,
-                    base_dh_inst=None):
+                    base_dh_inst=None,
+                    match_zero_rate=False):
         """Find the nominal scalar + adder combination to bias correct data
         at a single site"""
 
@@ -739,6 +792,9 @@ def _run_single(cls,
                                          decimals=decimals,
                                          base_dh_inst=base_dh_inst)
 
+        if match_zero_rate:
+            bias_data = cls._match_zero_rate(bias_data, base_data)
+
         out = cls.get_linear_correction(bias_data, base_data, bias_feature,
                                         base_dset)
         return out
@@ -931,6 +987,7 @@ def run(self,
                         self.bias_ti,
                         self.decimals,
                         base_dh_inst=self.base_dh,
+                        match_zero_rate=self.match_zero_rate,
                     )
                     for key, arr in single_out.items():
                         self.out[key][raster_loc] = arr
@@ -963,6 +1020,7 @@ def run(self,
                             daily_reduction,
                             self.bias_ti,
                             self.decimals,
+                            match_zero_rate=self.match_zero_rate,
                         )
                         futures[future] = raster_loc
 
@@ -1006,7 +1064,8 @@ def _run_single(cls,
                     daily_reduction,
                     bias_ti,
                     decimals,
-                    base_dh_inst=None):
+                    base_dh_inst=None,
+                    match_zero_rate=False):
         """Find the nominal scalar + adder combination to bias correct data
         at a single site"""
 
@@ -1018,6 +1077,9 @@ def _run_single(cls,
                                                decimals=decimals,
                                                base_dh_inst=base_dh_inst)
 
+        if match_zero_rate:
+            bias_data = cls._match_zero_rate(bias_data, base_data)
+
         base_arr = np.full(cls.NT, np.nan, dtype=np.float32)
         out = {}
 
@@ -1117,10 +1179,12 @@ def _init_out(self):
                        f'bias_{self.bias_feature}_std',
                        f'bias_{self.bias_feature}_skew',
                        f'bias_{self.bias_feature}_kurtosis',
+                       f'bias_{self.bias_feature}_zero_rate',
                        f'base_{self.base_dset}_mean',
                        f'base_{self.base_dset}_std',
                        f'base_{self.base_dset}_skew',
                        f'base_{self.base_dset}_kurtosis',
+                       f'base_{self.base_dset}_zero_rate',
                        ]
 
         self.out = {k: np.full((*self.bias_gid_raster.shape, self.NT),
@@ -1138,13 +1202,17 @@ def _init_out(self):
             self.out[bias_k] = arr.copy()
 
     @classmethod
-    def _run_skill_eval(cls, bias_data, base_data, bias_feature, base_dset):
+    def _run_skill_eval(cls, bias_data, base_data, bias_feature, base_dset,
+                        match_zero_rate=False):
         """Run skill assessment metrics on 1D datasets at a single site.
 
         Note we run the KS test on the mean=0 distributions as per:
         S. Brands et al. 2013 https://doi.org/10.1007/s00382-013-1742-8
         """
 
+        if match_zero_rate:
+            bias_data = cls._match_zero_rate(bias_data, base_data)
+
         out = {}
         bias_mean = np.nanmean(bias_data)
         base_mean = np.nanmean(base_data)
@@ -1154,13 +1222,20 @@ def _run_skill_eval(cls, bias_data, base_data, bias_feature, base_dset):
         out[f'bias_{bias_feature}_std'] = np.nanstd(bias_data)
         out[f'bias_{bias_feature}_skew'] = stats.skew(bias_data)
         out[f'bias_{bias_feature}_kurtosis'] = stats.kurtosis(bias_data)
+        out[f'bias_{bias_feature}_zero_rate'] = np.nanmean(bias_data == 0)
 
         out[f'base_{base_dset}_mean'] = base_mean
         out[f'base_{base_dset}_std'] = np.nanstd(base_data)
         out[f'base_{base_dset}_skew'] = stats.skew(base_data)
         out[f'base_{base_dset}_kurtosis'] = stats.kurtosis(base_data)
+        out[f'base_{base_dset}_zero_rate'] = np.nanmean(base_data == 0)
+
+        if match_zero_rate:
+            ks_out = stats.ks_2samp(base_data, bias_data)
+        else:
+            ks_out = stats.ks_2samp(base_data - base_mean,
+                                    bias_data - bias_mean)
 
-        ks_out = stats.ks_2samp(base_data - base_mean, bias_data - bias_mean)
         out[f'{bias_feature}_ks_stat'] = ks_out.statistic
         out[f'{bias_feature}_ks_p'] = ks_out.pvalue
 
@@ -1175,7 +1250,7 @@ def _run_skill_eval(cls, bias_data, base_data, bias_feature, base_dset):
     @classmethod
     def _run_single(cls, bias_data, base_fps, bias_feature, base_dset,
                     base_gid, base_handler, daily_reduction, bias_ti,
-                    decimals, base_dh_inst=None):
+                    decimals, base_dh_inst=None, match_zero_rate=False):
         """Do a skill assessment at a single site"""
 
         base_data, base_ti = cls.get_base_data(base_fps, base_dset,
@@ -1189,12 +1264,11 @@ def _run_single(cls, bias_data, base_fps, bias_feature, base_dset,
                f'bias_{bias_feature}_std_monthly': arr.copy(),
                f'base_{base_dset}_mean_monthly': arr.copy(),
                f'base_{base_dset}_std_monthly': arr.copy(),
-               f'{bias_feature}_scalar': arr.copy(),
-               f'{bias_feature}_adder': arr.copy(),
                }
 
         out.update(cls._run_skill_eval(bias_data, base_data,
-                                       bias_feature, base_dset))
+                                       bias_feature, base_dset,
+                                       match_zero_rate=match_zero_rate))
 
         for month in range(1, 13):
             bias_mask = bias_ti.month == month
@@ -1208,6 +1282,6 @@ def _run_single(cls, bias_data, base_fps, bias_feature, base_dset,
                 for k, v in mout.items():
                     if not k.endswith(('_scalar', '_adder')):
                         k += '_monthly'
-                    out[k][month - 1] = v
+                        out[k][month - 1] = v
 
         return out

tests/bias/test_bias_correction.py

@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 """pytests bias correction calculations"""
 import os
+import shutil
 import tempfile
 
 import h5py
@@ -526,7 +527,46 @@ def test_nc_base_file():
 
     out = calc.run(fill_extend=True, max_workers=1)
 
-    assert (out['rsds_scalar'] == 1).all()
-    assert (out['rsds_adder'] == 0).all()
+    assert np.allclose(out['base_rsds_mean_monthly'],
+                       out['bias_rsds_mean_monthly'])
     assert np.allclose(out['base_rsds_mean'], out['bias_rsds_mean'])
     assert np.allclose(out['base_rsds_std'], out['bias_rsds_std'])
+
+
+def test_match_zero_rate():
+    """Test feature to match the rate of zeros in the bias data based on the
+    zero rate in the base data. Useful for precip where GCMs have a low-precip
+    "drizzle" problem."""
+    bias_data = np.random.uniform(0, 1, 1000)
+    base_data = np.random.uniform(0, 1, 1000)
+    base_data[base_data < 0.1] = 0
+
+    skill = SkillAssessment._run_skill_eval(bias_data, base_data, 'f1', 'f1')
+    assert skill['bias_f1_zero_rate'] != skill['base_f1_zero_rate']
+    assert (bias_data == 0).mean() != (base_data == 0).mean()
+
+    skill = SkillAssessment._run_skill_eval(bias_data, base_data, 'f1', 'f1',
+                                            match_zero_rate=True)
+    assert (bias_data == 0).mean() == (base_data == 0).mean()
+    assert skill['bias_f1_zero_rate'] == skill['base_f1_zero_rate']
+    for p in (1, 5, 25, 50, 75, 95, 99):
+        assert np.allclose(skill[f'base_f1_percentile_{p}'],
+                           np.percentile(base_data, p))
+
+    with tempfile.TemporaryDirectory() as td:
+        fp_nsrdb_temp = os.path.join(td, os.path.basename(FP_NSRDB))
+        shutil.copy(FP_NSRDB, fp_nsrdb_temp)
+        with h5py.File(fp_nsrdb_temp, 'a') as nsrdb_temp:
+            ghi = nsrdb_temp['ghi'][...]
+            ghi[:1000, :] = 0
+            nsrdb_temp['ghi'][...] = ghi
+        calc = SkillAssessment(fp_nsrdb_temp, FP_CC, 'ghi', 'rsds',
+                               target=TARGET, shape=SHAPE,
+                               distance_upper_bound=0.7,
+                               bias_handler='DataHandlerNCforCC',
+                               match_zero_rate=True)
+        out = calc.run(fill_extend=True, max_workers=1)
+
+    bias_rate = out['bias_rsds_zero_rate']
+    base_rate = out['base_ghi_zero_rate']
+    assert np.allclose(bias_rate, base_rate, rtol=0.005)