Fixes for ERAD 2022 short course (#293)

* Implement support for different x- and y-pixel sizes in aggregate_fields_space

* Add tests for different x- and y-window sizes in aggregate_fields_space

* Fix possible numerical instability when checking the spatial window size

* Another fix for numerical instability in window size checks

* Replace int with round

* Update pysteps rc to correct RMI radar path names

* Fix typo in pysteps rc

* Add check for non-finite values

* Add the VET paper by Laroche and Zawadzki to bibliography

* Fix typo

* Remove unnecessary log base argument (default=10, throws exception with older matplotlib versions)

Co-authored-by: Ruben Imhoff <31476760+RubenImhoff@users.noreply.github.com>

doc/source/references.bib

@@ -147,6 +147,17 @@ @article{Hwang2015
   DOI = "10.1175/WAF-D-15-0057.1"
 }
 
+@ARTICLE{LZ1995,
+  AUTHOR = "S. Laroche and I. Zawadzki",
+  TITLE = "Retrievals of Horizontal Winds from Single-Doppler Clear-Air Data by Methods of Cross Correlation and Variational Analysis",
+  JOURNAL = "Journal of Atmospheric and Oceanic Technology",
+  VOLUME = 12,
+  NUMBER = 4,
+  PAGES = "721--738",
+  YEAR = 1995,
+  DOI = "10.1175/1520-0426(1995)012<0721:ROHWFS>2.0.CO;2",
+}
+
 @ARTICLE{NBSG2017,
   AUTHOR = "D. Nerini and N. Besic and I. Sideris and U. Germann and L. Foresti",
   TITLE = "A non-stationary stochastic ensemble generator for radar rainfall fields based on the short-space {F}ourier transform",

pysteps/extrapolation/semilagrangian.py

@@ -80,7 +80,7 @@ def extrapolate(
         Default: False
     vel_timestep: float
         The time step of the velocity field. It is assumed to have the same
-        unit as the timesteps argument. Applicable if timeseps is a list.
+        unit as the timesteps argument. Applicable if timesteps is a list.
         Default: 1.
     interp_order: int
         The order of interpolation to use. Default: 1 (linear). Setting this

pysteps/motion/darts.py

@@ -95,6 +95,9 @@ def DARTS(input_images, **kwargs):
             "invalid output_type=%s, must be 'spatial' or 'spectral'" % output_type
         )
 
+    if np.any(~np.isfinite(input_images)):
+        raise ValueError("the input images contain non-finite values")
+
     if verbose:
         print("Computing the motion field with the DARTS method.")
         t0 = time.time()

pysteps/pystepsrc

@@ -91,7 +91,7 @@
 		"rmi": {
             "root_path": "./radar/rmi/radqpe",
             "path_fmt": "%Y%m%d",
-            "fn_pattern": "%Y%m%d%H%M00.rad.bhbjbwdnfa.comp.rate.qpe2",
+            "fn_pattern": "%Y%m%d%H%M00.rad.best.comp.rate.qpe",
             "fn_ext": "hdf",
             "importer": "odim_hdf5",
             "timestep": 5,

pysteps/tests/test_utils_dimension.py

@@ -152,6 +152,13 @@ def test_aggregate_fields_time(R, metadata, time_window_min, ignore_nan, expecte
         False,
         4 * np.ones((1, 5, 5)),
     ),
+    (
+        np.ones((1, 10, 10)),
+        {"unit": "mm/h", "xpixelsize": 1, "ypixelsize": 2},
+        (2, 4),
+        False,
+        np.ones((1, 5, 5)),
+    ),
 ]
 
 

pysteps/utils/dimension.py

@@ -130,23 +130,23 @@ def aggregate_fields_space(R, metadata, space_window, ignore_nan=False):
         Metadata dictionary containing the xpixelsize, ypixelsize and unit
         attributes as described in the documentation of
         :py:mod:`pysteps.io.importers`.
-    space_window: float or None
+    space_window: float, tuple or None
         The length of the space window that is used to upscale the fields.
-        The space_window unit is the same used in the geographical projection
-        of R
-        and hence the same as for the xpixelsize and ypixelsize attributes.
-        The space spanned by the m and n dimensions of R must be a multiple of
-        space_window.
-        If set to None, it returns a copy of the original R and metadata.
+        If a float is given, the same window size is used for the x- and
+        y-directions. Separate window sizes are used for x- and y-directions if
+        a two-element tuple is given. The space_window unit is the same used in
+        the geographical projection of R and hence the same as for the xpixelsize
+        and ypixelsize attributes. The space spanned by the n- and m-dimensions
+        of R must be a multiple of space_window. If set to None, the function
+        returns a copy of the original R and metadata.
     ignore_nan: bool, optional
         If True, ignore nan values.
 
     Returns
     -------
     outputarray: array-like
-        The new array of aggregated fields of shape (k,j), (t,k,j)
-        or (l,t,k,j),
-        where k = m*ypixelsize/space_window and j = n*xpixelsize/space_window.
+        The new array of aggregated fields of shape (k,j), (t,k,j) or (l,t,k,j),
+        where k = m*ypixelsize/space_window[1] and j = n*xpixelsize/space_window[0].
     metadata: dict
         The metadata with updated attributes.
 
@@ -177,15 +177,23 @@ def aggregate_fields_space(R, metadata, space_window, ignore_nan=False):
     if len(R.shape) > 4:
         raise ValueError("The number of dimensions must be <= 4")
 
+    if np.isscalar(space_window):
+        space_window = (space_window, space_window)
+
     # assumes that frames are evenly spaced
-    if ypixelsize == space_window and xpixelsize == space_window:
+    if ypixelsize == space_window[1] and xpixelsize == space_window[0]:
         return R, metadata
-    if (R.shape[axes[0]] * ypixelsize) % space_window or (
-        R.shape[axes[1]] * xpixelsize
-    ) % space_window:
+
+    ysize = R.shape[axes[0]] * ypixelsize
+    xsize = R.shape[axes[1]] * xpixelsize
+
+    if (
+        abs(ysize / space_window[1] - round(ysize / space_window[1])) > 1e-10
+        or abs(xsize / space_window[0] - round(xsize / space_window[0])) > 1e-10
+    ):
         raise ValueError("space_window does not equally split R")
 
-    nframes = [int(space_window / ypixelsize), int(space_window / xpixelsize)]
+    nframes = [int(space_window[1] / ypixelsize), int(space_window[0] / xpixelsize)]
 
     # specify the operator to be used to aggregate the values
     # within the space window
@@ -204,8 +212,8 @@ def aggregate_fields_space(R, metadata, space_window, ignore_nan=False):
     R = aggregate_fields(R, nframes[0], axis=axes[0], method=method)
     R = aggregate_fields(R, nframes[1], axis=axes[1], method=method)
 
-    metadata["ypixelsize"] = space_window
-    metadata["xpixelsize"] = space_window
+    metadata["ypixelsize"] = space_window[1]
+    metadata["xpixelsize"] = space_window[0]
 
     return R, metadata
 

pysteps/verification/plots.py

@@ -160,7 +160,7 @@ def plot_reldiag(reldiag, ax=None):
         color="gray",
         edgecolor="black",
     )
-    iax.set_yscale("log", base=10)
+    iax.set_yscale("log")
     iax.set_xticks(reldiag["bin_edges"])
     iax.set_xticklabels(["%.1f" % max(v, 1e-6) for v in reldiag["bin_edges"]])
     yt_min = int(max(np.floor(np.log10(min(reldiag["sample_size"][:-1]))), 1))