Merge branch 'master' of https://github.com/pySTEPS/pysteps

Author: dnerini

examples/ensemble_verification.py

@@ -62,9 +62,7 @@
     "n_cascade_levels"  : [6],
     "noise_adjustment"  : [True],
     "conditional"       : [False],
-    "precip_mask"       : [True],
     "mask_method"       : ["incremental"],      # obs, incremental, sprog
-    "prob_matching"     : [True],
 }
 
 # Conditional parameters
@@ -252,9 +250,7 @@ def export(X):
                             noise_method=p["noise_method"], 
                             noise_stddev_adj=p["noise_adjustment"],
                             ar_order=p["ar_order"],conditional=p["conditional"], 
-                            use_probmatching=p["prob_matching"], 
-                            mask_method=p["mask_method"], 
-                            use_precip_mask=p["precip_mask"], callback=export, 
+                            mask_method=p["mask_method"], callback=export, 
                             return_output=False, seed=p["seed"])
 
             ## save results

examples/run_ensemble_nowcast.py

@@ -53,8 +53,6 @@
 ar_order            = 2
 r_threshold         = 0.1               # rain/no-rain threshold [mm/h]
 adjust_noise        = True
-prob_matching       = True
-precip_mask         = True
 mask_method         = "incremental"     # sprog, obs or incremental
 conditional         = False
 unit                = "mm/h"            # mm/h or dBZ
@@ -117,8 +115,7 @@
                    bandpass_filter_method=bandpass_filter,
                    noise_method=noise_method, noise_stddev_adj=adjust_noise,
                    ar_order=ar_order, conditional=conditional,
-                   use_precip_mask=precip_mask, mask_method=mask_method,
-                   use_probmatching=prob_matching, seed=seed)
+                   mask_method=mask_method, seed=seed)
 
 ## if necessary, transform back all data
 R_fct, _    = transformer(R_fct, metadata, inverse=True)

pysteps/io/archive.py

@@ -15,8 +15,9 @@ def find_by_date(date, root_path, path_fmt, fn_pattern, fn_ext, timestep,
     root_path : str
         The root path to search the input files.
     path_fmt : str
-        Path format. It may consist of directory names separated by '/' and
-        date/time specifiers beginning with '%' (e.g. %Y/%m/%d).
+        Path format. It may consist of directory names separated by '/',
+        date/time specifiers beginning with '%' (e.g. %Y/%m/%d) and wildcards
+        (?) that match any single character.
     fn_pattern : str
         The name pattern of the input files without extension. The pattern can
         contain time specifiers (e.g. %H, %M and %S).

pysteps/io/importers.py

@@ -128,7 +128,10 @@ def import_bom_rf3(filename, **kwargs):
     metadata["unit"]        = "mm/h"
     metadata["transform"]   = None
     metadata["zerovalue"]   = np.nanmin(R)
-    metadata["threshold"]   = np.nanmin(R[R>np.nanmin(R)])
+    if np.any(np.isfinite(R)):
+        metadata["threshold"] = np.nanmin(R[R>np.nanmin(R)])
+    else:
+        metadata["threshold"] = np.nan
 
     return R, None, metadata
 
@@ -255,7 +258,10 @@ def import_fmi_pgm(filename, **kwargs):
     metadata["unit"]        = "dBZ"
     metadata["transform"]   = "dB"
     metadata["zerovalue"]   = np.nanmin(R)
-    metadata["threshold"]   = np.nanmin(R[R>np.nanmin(R)])
+    if np.any(np.isfinite(R)):
+        metadata["threshold"] = np.nanmin(R[R>np.nanmin(R)])
+    else:
+        metadata["threshold"] = np.nan
 
     return R,None,metadata
 
@@ -431,7 +437,7 @@ def import_mch_gif(filename, **kwargs):
     if np.any(R>np.nanmin(R)):
         metadata["threshold"]   = np.nanmin(R[R>np.nanmin(R)])
     else:
-        metadata["threshold"]   = None
+        metadata["threshold"]   = np.nan
     metadata["institution"] = "MeteoSwiss"
     metadata["product"] = product
 
@@ -538,14 +544,19 @@ def import_mch_hdf5(filename, **kwargs):
         unit = "mm/h"
         transform = None
 
+    if np.any(np.isfinite(R)):
+        thr = np.nanmin(R[R>np.nanmin(R)])
+    else:
+        thr = np.nan
+
     metadata.update({
                 "yorigin":"upper",
                 "institution":"MeteoSwiss",
                 "accutime":5.,
                 "unit":unit,
                 "transform":transform,
                 "zerovalue":np.nanmin(R),
-                "threshold":np.nanmin(R[R>np.nanmin(R)]) })
+                "threshold":thr })
 
     f.close()
 
@@ -778,6 +789,11 @@ def import_odim_hdf5(filename, **kwargs):
         unit = "mm/h"
         transform = None
 
+    if np.any(np.isfinite(R)):
+        thr = np.nanmin(R[R>np.nanmin(R)])
+    else:
+        thr = nan
+    
     metadata = {"projection":proj4str,
                 "ll_lon":LL_lon,
                 "ll_lat":LL_lat,
@@ -795,7 +811,7 @@ def import_odim_hdf5(filename, **kwargs):
                 "unit":unit,
                 "transform":transform,
                 "zerovalue":np.nanmin(R),
-                "threshold":np.nanmin(R[R>np.nanmin(R)])}
+                "threshold":thr}
 
     f.close()
 

pysteps/noise/motion.py

@@ -38,8 +38,8 @@ def initialize_bps(V, pixelsperkm, timestep, p_pert_par=None, p_pert_perp=None,
                    randstate=np.random, seed=None):
     """Initialize the motion field perturbator described in :cite:`BPS2006`.
     For simplicity, the bias adjustment procedure described there has not been
-    implemented. The perturbator generates a constant field whose magnitude
-    depends on lead time.
+    implemented. The perturbator generates a field whose magnitude increases
+    with respect to lead time.
 
     Parameters
     ----------

pysteps/nowcasts/steps.py

@@ -19,10 +19,10 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
              kmperpixel=None, timestep=None, extrap_method="semilagrangian",
              decomp_method="fft", bandpass_filter_method="gaussian",
              noise_method="nonparametric", noise_stddev_adj=False, ar_order=2,
-             vel_pert_method=None, conditional=False, use_precip_mask=True,
-             use_probmatching=True, mask_method="incremental", callback=None,
-             return_output=True, seed=None, num_workers=None, fft_method="numpy",
-             extrap_kwargs={}, filter_kwargs={}, noise_kwargs={}, vel_pert_kwargs={}):
+             vel_pert_method="bps", conditional=False, probmatching_method="cdf",
+             mask_method="incremental", callback=None, return_output=True,
+             seed=None, num_workers=None, fft_method="numpy", extrap_kwargs={},
+             filter_kwargs={}, noise_kwargs={}, vel_pert_kwargs={}):
     """Generate a nowcast ensemble by using the Short-Term Ensemble Prediction
     System (STEPS) method.
 
@@ -47,7 +47,7 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
     ----------------
     R_thr : float
       Specifies the threshold value for minimum observable precipitation
-      intensity. Must be set if use_probmatching is True or conditional is True.
+      intensity. Required if mask_method is not None or conditional is True.
     kmperpixel : float
       Spatial resolution of the input data (kilometers/pixel). Required if
       vel_pert_method is not None or mask_method is 'incremental'.
@@ -78,21 +78,20 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
       field is not perturbed.
     conditional : bool
       If set to True, compute the statistics of the precipitation field
-      conditionally by excluding the areas where the values are below the
-      threshold R_thr.
-    use_precip_mask : bool
-      If True, set pixels outside precipitation areas to the minimum value of
-      the observed field.
-    mask_method : {'obs', 'sprog', 'incremental'}
-      The precipitation/no precipitation method to use with mask: 'obs' = apply R_thr
-      to the most recently observed precipitation intensity field, 'sprog' = use the
-      smoothed forecast field from S-PROG, where the AR(p) model has been applied,
-      'incremental' = iteratively buffer the mask with a certain rate (currently
-      it is 1 km/min)
-    use_probmatching : bool
-      If True, apply probability matching to the forecast field in order to
-      preserve the distribution of the most recently observed precipitation
-      field.
+      conditionally by excluding pixels where the values are below the threshold
+      R_thr.
+    mask_method : {'obs','sprog','incremental',None}
+      The method to use for masking no precipitation areas in the forecast field. 
+      The masked pixels are set to the minimum value of the observations.
+      'obs' = apply R_thr to the most recently observed precipitation intensity
+      field, 'sprog' = use the smoothed forecast field from S-PROG, where the
+      AR(p) model has been applied, 'incremental' = iteratively buffer the mask
+      with a certain rate (currently it is 1 km/min), None=no masking.
+    probmatching_method : {'cdf','mean',None}
+      Method for matching the statistics of the forecast field with those of
+      the most recently observed one. Requires that mask_method is not None.
+      'cdf'=map the forecast CDF to the observed one, 'mean'=adjust only the
+      mean value of the forecast field, None=no matching applied.
     callback : function
       Optional function that is called after computation of each time step of
       the nowcast. The function takes one argument: a three-dimensional array
@@ -110,19 +109,19 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
       all available CPUs. Applicable if dask is enabled.
     fft_method : str or tuple
       A string or a (function,kwargs) tuple defining the FFT method to use
-      (see utils.fft.get_method). Defaults to "numpy".
+      (see utils.fft.get_method). Defaults to 'numpy'.
     extrap_kwargs : dict
-      Optional dictionary that is supplied as keyword arguments to the
-      extrapolation method.
+      Optional dictionary containing keyword arguments for the extrapolation
+      method. See the documentation of pysteps.extrapolation.
     filter_kwargs : dict
-      Optional dictionary that is supplied as keyword arguments to the
-      filter method.
+      Optional dictionary containing keyword arguments for the filter method.
+      See the documentation of pysteps.cascade.bandpass_filters.py.
     noise_kwargs : dict
-      Optional dictionary that is supplied as keyword arguments to the
-      initializer of the noise generator.
+      Optional dictionary containing keyword arguments for the initializer of
+      the noise generator. See the documentation of pysteps.noise.fftgenerators.
     vel_pert_kwargs : dict
-      Optional dictionary that is supplied as keyword arguments to the
-      initializer of the velocity perturbator.
+      Optional dictionary containing keyword arguments for the initializer of
+      the velocity perturbator. See the documentation of pysteps.noise.motion.
 
     Returns
     -------
@@ -139,9 +138,9 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
 
     Notes
     -----
-    If noise_method and vel_pert_method are set to None and n_ens_members is set
-    to 1, the produced nowcast is deterministic (i.e. the S-PROG nowcast, see
-    :cite:`Seed2003`).
+    If noise_method and vel_pert_method are None, n_ens_members is 1,
+    mask_method is 'sprog' and probmatching_method is 'mean', the deterministic
+    S-PROG nowcast is generated, see :cite:`Seed2003`.
 
     References
     ----------
@@ -156,14 +155,17 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
     if np.any(~np.isfinite(V)):
         raise ValueError("V contains non-finite values")
 
-    if mask_method not in ["obs", "sprog", "incremental"]:
-        raise ValueError("unknown mask method %s: must be 'obs', 'sprog' or 'incremental'" % mask_method)
+    if mask_method not in ["obs", "sprog", "incremental", None]:
+        raise ValueError("unknown mask method %s: must be 'obs', 'sprog' or 'incremental' or None" % mask_method)
 
     if conditional and R_thr is None:
         raise ValueError("conditional=True but R_thr is not set")
 
-    if use_probmatching and R_thr is None:
-        raise ValueError("use_probmatching=True but R_thr is not set")
+    if probmatching_method is not None and R_thr is None:
+        raise ValueError("probmatching_method!=None but R_thr is not set")
+
+    if probmatching_method is not None and mask_method is None:
+        raise ValueError("probmatching_method!=None but mask_method=None")
 
     if kmperpixel is None:
         if vel_pert_method is not None:
@@ -199,9 +201,8 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
     print("noise adjustment:       %s" % ("yes" if noise_stddev_adj else "no"))
     print("velocity perturbator:   %s" % vel_pert_method)
     print("conditional statistics: %s" % ("yes" if conditional else "no"))
-    print("precipitation mask:     %s" % ("yes" if use_precip_mask else "no"))
-    print("mask method:            %s" % mask_method)
-    print("probability matching:   %s" % ("yes" if use_probmatching else "no"))
+    print("precip. mask method:    %s" % mask_method)
+    print("probability matching:   %s" % probmatching_method)
     print("FFT method:             %s" % fft_method)
     print("")
 
@@ -219,14 +220,14 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
         print("velocity perturbations, perpendicular: %g,%g,%g" % \
             (vp_perp[0], vp_perp[1], vp_perp[2]))
 
-    if conditional or use_probmatching:
-        print("conditional precip. intensity threshold: %g" % R_thr)
+    if conditional or mask_method is not None:
+        print("precip. intensity threshold: %g" % R_thr)
 
     M,N = R.shape[1:]
     extrap_method = extrapolation.get_method(extrap_method)
     R = R[-(ar_order + 1):, :, :].copy()
 
-    if conditional or use_probmatching:
+    if conditional:
         MASK_thr = np.logical_and.reduce([R[i, :, :] >= R_thr for i in range(R.shape[0])])
     else:
         MASK_thr = None
@@ -342,29 +343,32 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
     D = [None for j in range(n_ens_members)]
     R_f = [[] for j in range(n_ens_members)]
 
-    if use_precip_mask:
+    if mask_method is not None:
+        MASK_prec = R[-1, :, :] >= R_thr
+
         if mask_method == "obs":
-            MASK_prec = R[-1, :, :] >= R_thr
+            pass
             # add a slight buffer to the mask
             # n=5
             # kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (n,n))
             # MASK_prec = MASK_prec.astype('uint8')
             # MASK_prec = cv2.dilate(MASK_prec,kernel).astype(bool)
         elif mask_method == "sprog":
             # compute the wet area ratio and the precipitation mask
-            MASK_prec = R[-1, :, :] >= R_thr
             war = 1.0*np.sum(MASK_prec) / (R.shape[1]*R.shape[2])
             R_m = R_c[0, :, :, :].copy()
         elif mask_method == "incremental":
             # initialize precip mask for each member
-            MASK_prec_ = R[-1, :, :] >= R_thr
-            MASK_prec = [MASK_prec_.copy() for j in range(n_ens_members)]
+            MASK_prec = [MASK_prec.copy() for j in range(n_ens_members)]
             # initialize the structuring element
             struct = scipy.ndimage.generate_binary_structure(2, 1)
             # iterate it to expand it nxn
             n = timestep/kmperpixel
             struct = scipy.ndimage.iterate_structure(struct, int((n - 1)/2.))
 
+        if probmatching_method == "mean":
+            mu_0 = np.mean(R[-1, :, :][MASK_prec])
+
     R = R[-1, :, :]
 
     print("Starting nowcast computation.")
@@ -375,7 +379,7 @@ def forecast(R, V, n_timesteps, n_ens_members=24, n_cascade_levels=6, R_thr=None
         sys.stdout.flush()
         starttime = time.time()
 
-        if use_precip_mask and mask_method == "sprog":
+        if mask_method == "sprog":
             for i in range(n_cascade_levels):
                 # use a separate AR(p) model for the non-perturbed forecast,
                 # from which the mask is obtained
@@ -434,24 +438,29 @@ def worker(j):
             # obtained from the AR(p) model(s)
             R_c_ = _recompose_cascade(R_c[j, :, :, :], mu, sigma)
 
-            if use_precip_mask:
+            if mask_method is not None:
                 # apply the precipitation mask to prevent generation of new
                 # precipitation into areas where it was not originally
                 # observed
                 if mask_method == "obs":
-                    R_c_[~MASK_prec] = R_c_.min()
+                    MASK_prec_ = ~MASK_prec
                 elif mask_method == "incremental":
-                    R_c_[~MASK_prec[j]] = R_c_.min()
+                    MASK_prec_ = ~MASK_prec[j]
                 elif mask_method == "sprog":
-                    R_c_[MASK_prec] = R_c_.min()
+                    MASK_prec_ = MASK_prec
 
-            if use_probmatching:
-                ## adjust the conditional CDF of the forecast (precipitation
-                ## intensity above the threshold R_thr) to match the most
-                ## recently observed precipitation field
+                R_c_[MASK_prec_] = R_c_.min()
+
+            if probmatching_method == "cdf":
+                # adjust the conditional CDF of the forecast (precipitation
+                # intensity above the threshold R_thr) to match the most
+                # recently observed precipitation field
                 R_c_ = probmatching.nonparam_match_empirical_cdf(R_c_, R)
+            elif probmatching_method == "mean":
+                mu_fct = np.mean(R_c_[~MASK_prec_])
+                R_c_[~MASK_prec_] = R_c_[~MASK_prec_] - mu_fct + mu_0
 
-            if use_precip_mask and mask_method == "incremental":
+            if mask_method == "incremental":
                 MASK_prec_ = R_c_ >= R_thr
                 MASK_prec_ = scipy.ndimage.morphology.binary_dilation(MASK_prec_, struct)
                 MASK_prec[j] = MASK_prec_
@@ -493,10 +502,7 @@ def worker(j):
                 R_f[j].append(R_f_[j])
 
     if return_output:
-        if n_ens_members == 1:
-            return np.stack(R_f[0])
-        else:
-            return np.stack([np.stack(R_f[j]) for j in range(n_ens_members)])
+        return np.stack([np.stack(R_f[j]) for j in range(n_ens_members)])
     else:
         return None