Add parametric nowcasting and calibrated stochastic rainfall generation for v2

MANIFEST.in

@@ -4,6 +4,5 @@ include pysteps/pystepsrc_schema.json
 include pysteps/io/mch_lut_8bit_Metranet_AZC_V104.txt
 include pysteps/io/mch_lut_8bit_Metranet_v103.txt
 recursive-include pysteps *.pyx
+recursive-include pysteps/tests/data *.json *.npz *.rst
 include pyproject.toml
-
-

doc/source/index.rst

@@ -34,6 +34,9 @@ __ https://github.com/pySTEPS/pysteps
    Gallery <../auto_examples/index>
    My first nowcast (Colab Notebook) <https://colab.research.google.com/github/pySTEPS/pysteps/blob/master/examples/my_first_nowcast.ipynb>
    API Reference <pysteps_reference/index>
+   Parametric rainfall estimators <user_guide/parametric_estimators>
+   Calibrated stochastic rainfall fields <user_guide/calibrated_generator>
+   Parametric nowcasting and simulation <user_guide/parametric_nowcast>
    Example data <user_guide/example_data>
    Configuration file (pystepsrc) <user_guide/set_pystepsrc>
    Machine learning applications <user_guide/machine_learning_pysteps>

doc/source/pysteps_reference/noise.rst

@@ -2,10 +2,11 @@
 pysteps.noise
 =============
 
-Implementation of deterministic and ensemble nowcasting methods.
+Methods for generating spatially correlated noise and velocity perturbations.
 
 
-.. automodule:: pysteps.noise.interface
 .. automodule:: pysteps.noise.fftgenerators
+.. automodule:: pysteps.noise.interface
 .. automodule:: pysteps.noise.motion
+.. automodule:: pysteps.noise.parametric
 .. automodule:: pysteps.noise.utils

doc/source/pysteps_reference/nowcasts.rst

@@ -8,8 +8,9 @@ Implementation of deterministic and ensemble nowcasting methods.
 .. automodule:: pysteps.nowcasts.interface
 .. automodule:: pysteps.nowcasts.anvil
 .. automodule:: pysteps.nowcasts.extrapolation
-.. automodule:: pysteps.nowcasts.linda
 .. automodule:: pysteps.nowcasts.lagrangian_probability
+.. automodule:: pysteps.nowcasts.linda
+.. automodule:: pysteps.nowcasts.parametric
 .. automodule:: pysteps.nowcasts.sprog
 .. automodule:: pysteps.nowcasts.sseps
 .. automodule:: pysteps.nowcasts.steps

doc/source/pysteps_reference/timeseries.rst

@@ -7,4 +7,5 @@ Methods and models for time series analysis.
 
 .. automodule:: pysteps.timeseries.autoregression
 .. automodule:: pysteps.timeseries.correlation
+.. automodule:: pysteps.timeseries.parametric
 

doc/source/pysteps_reference/utils.rst

@@ -16,5 +16,6 @@ Implementation of miscellaneous utility functions.
 .. automodule:: pysteps.utils.pca
 .. automodule:: pysteps.utils.reprojection
 .. automodule:: pysteps.utils.spectral
+.. automodule:: pysteps.utils.statistics
 .. automodule:: pysteps.utils.tapering
 .. automodule:: pysteps.utils.transformation

doc/source/user_guide/calibrated_generator.rst

@@ -0,0 +1,121 @@
+Calibrated stochastic rainfall fields
+=====================================
+
+The calibrated generator creates random rainfall patterns from a rain fraction,
+a rain-intensity distribution, and spectral slopes describing spatial structure.
+It returns a NumPy array in transformed rainfall units. In parametric nowcasting,
+these fields provide new spatial variability at each forecast step.
+
+Here, the marginal distribution means the distribution of values across pixels,
+without regard to their positions. Spectral slopes describe how variability is
+distributed across spatial scales. See :ref:`sphx_glr_auto_examples_plot_calibrated_noise_generator.py`
+for a radar-data example using rainfall statistics, parameter selection,
+calibration diagnostics, rainfall distributions, spectra, and normalized
+cascade innovations. It uses the same Swiss composite as the original
+stochastic-noise tutorial on a complete subdomain.
+
+Generating a field
+------------------
+
+Supply signed spectral slopes, rain fraction, and caller-selected Normal location
+and scale in transformed rainfall units. The example values illustrate the
+interface; the parametric engine estimates and selects these parameters from
+the observations::
+
+    import numpy as np
+    from pysteps.noise.parametric import generate_calibrated_field
+
+    parameters = {
+        "beta_1": -2.4,
+        "beta_2": -3.1,
+        "rain_fraction": 0.08,
+        "mean_db": 6.0,
+        "stdev_db": 3.5,
+    }
+    field, diagnostics = generate_calibrated_field(
+        (256, 256), parameters,
+        kmperpixel=1.0, scale_break=15.0, threshold=0.0,
+        seed=0, return_diagnostics=True,
+    )
+    assert field.shape == (256, 256) and np.isfinite(field).all()
+
+``shape`` is in ``(y, x)`` order; pixel spacing and break wavelength are in km.
+The threshold is in the field's transformed units, so 1 mm/h rainfall occurrence
+corresponds to a 0 dBR threshold. Location and scale are typically latent Normal
+parameters, but this helper consumes the active moments supplied by the caller
+and performs no latent-versus-ordinary moment selection or regional slope repair.
+The input mapping is not modified.
+
+The default result is the field alone. Optional diagnostics report estimated
+spectral slopes, distance, number of attempts, selected attempt, tolerance status,
+target and calibrated parameters, settings, and all candidate estimates and
+distances. A finite returned field does not imply that the requested spectral
+tolerance was met.
+
+Generation and spectral correction
+----------------------------------
+
+Independent ``U(0,1)`` samples are FFT-filtered using half the signed power slopes
+as amplitude exponents. The resulting correlated field is mapped to the target
+rain-intensity distribution. This mapping also changes its spectrum. To
+compensate, both filter slopes receive an additive correction, default -0.5,
+before generation; the Normal location and scale are used as supplied.
+
+The histogram marginal mapper uses 500 bins and clips probabilities to
+``[0.001, 0.999]``. It selects a wet-area cutoff using a bin's upper edge and maps
+the selected sample to Normal quantiles using bin-centre interpolation. The
+returned dry anomalies remain continuous. Mapping can change the realised wet
+fraction and spectrum, and the clipped tails constrain extremes; occurrence and
+marginal matching are approximate.
+
+Spectra are assessed on a copy with values at or below the transformed rain
+threshold replaced by the dry value. The distance is
+``sqrt(sum((estimated_beta-target_beta)**2/beta_variances))`` with default
+variances ``(0.067738, 0.076405)``. These are empirical weights. Fractions at or
+below 0.1 permit at most five attempts; higher fractions use one. The first
+candidate with distance at or below 1 is accepted. Otherwise the closest finite
+candidate is selected, retaining the first exact tie. If all spectral assessments
+are unavailable, the last field is returned with an explicit failure status.
+
+The correction, weights, acceptance threshold, retry settings, bin count, and
+quantile clip can all be overridden. The correction and error weights are
+empirical and may need calibration for other rainfall data.
+
+Random streams and dry cases
+----------------------------
+
+``seed`` creates a local NumPy random stream without changing global random state.
+Alternatively, pass a caller-owned ``rng``; it advances in place and must be
+independent for each ensemble member. Supply ``rng`` or ``seed``. Retry counts
+affect stream consumption, so reproducibility requires matching parameters,
+settings, and numerical implementation as well as the initial seed.
+
+An explicit zero rain fraction produces a constant dry field with zero attempts
+and no spectral-success claim. Positive-fraction mappings with too few distinct
+wet pixels raise an error. Missing/nonfinite input and invalid parameters also
+raise errors; they are not classified as dry rainfall.
+
+Cascade innovations and marginal adjustment
+-------------------------------------------
+
+The raw field retains its rainfall marginal. The parametric engine generates it
+first, then obtains normalized cascade innovations through decomposition::
+
+    from pysteps.cascade.bandpass_filters import filter_gaussian
+    from pysteps.cascade.decomposition import decomposition_fft
+
+    bandpass = filter_gaussian(field.shape, 6)
+    innovations = decomposition_fft(
+        field, bandpass, compute_stats=True, normalize=True
+    )
+    assert np.isfinite(innovations["cascade_levels"]).all()
+
+The raw field is not guaranteed to have global mean zero and variance one.
+Existing native STEPS noise-method registrations retain their perturbation
+interface; the calibrated parametric engine consumes this helper directly.
+
+The shared :func:`pysteps.postprocessing.probmatching.match_rainfall_distribution`
+also supports the post-AR marginal adjustment. Its default suppresses no positive
+rain fraction. The forecast passes ``min_rain_fraction=0.025`` to suppress
+very small target fractions after an AR update. Stochastic generation itself
+does not apply that cutoff.

doc/source/user_guide/parametric_estimators.rst

@@ -0,0 +1,115 @@
+Parametric rainfall estimators
+==============================
+
+The parametric estimators provide NumPy-based spatial statistics and temporal
+model utilities. They support the three-observation forecast interface described in :doc:`parametric_nowcast`
+and can also be used independently.
+
+Spatial estimation
+------------------
+
+Use a finite rainfall field and retain the threshold and dry-value convention
+when transforming it. Occurrence and marginal fitting use explicitly separate
+cutoffs: one identifies rainy pixels, and the other selects the intensity
+samples used for fitting. This example uses illustrative settings::
+
+    import numpy as np
+    import xarray as xr
+    from pysteps.utils.transformation import dB_transform
+    from pysteps.utils.statistics import compute_field_statistics
+    from pysteps.utils.spectral import rapsd_logarithmic, fit_power_law_spectrum
+
+    rng = np.random.default_rng(42)
+    rain = rng.lognormal(1.0, 0.8, size=(256, 256))
+    rain[rng.uniform(size=rain.shape) > 0.3] = 0.0
+    dataset = xr.Dataset(
+        {"rain": (("y", "x"), rain, {"units": "mm/h"})},
+        attrs={"precip_var": "rain"},
+    )
+    transformed_dataset = dB_transform(dataset, threshold=1.0, zerovalue=-0.1)
+    transformed = transformed_dataset["rain"].values
+    metadata = transformed_dataset["rain"].attrs
+    statistics = compute_field_statistics(
+        transformed, metadata["threshold"], lower_truncation=1.25
+    )
+    frequency, power = rapsd_logarithmic(transformed)
+    pixel_size_km = 1.0
+    scale_break_km = 15.0
+    spectrum_fit = fit_power_law_spectrum(
+        frequency, power, scale_break=scale_break_km / pixel_size_km
+    )
+
+``ordinary_mean`` and ``ordinary_stdev`` describe pixels at or above the
+occurrence threshold. ``latent_mean`` and ``latent_stdev`` estimate the parent
+Normal from samples at or above the fitting cutoff. They are different estimates;
+the fitting cutoff does not change the wet fraction. ``nonzero_fit`` reports
+sample count and fitting status, including finite bounded solutions whose
+parameters should not be interpreted as an unrestricted interior fit.
+
+The radial estimator subtracts the field mean, divides by its standard
+deviation, and averages unscaled FFT power in logarithmic frequency bins. It differs from
+:func:`pysteps.utils.spectral.rapsd`, which uses different binning and
+normalization. Its returned arrays contain linear frequencies in cycles per
+pixel and mean linear power per Fourier mode, not annulus-integrated energy.
+The fit uses signed slopes, an optional prescribed break, and a continuous
+piecewise-linear model in 10*log10 units. The two initial bins and frequencies
+at or above 10**(-0.4) cycles per pixel are excluded by default. The fixed-break
+model is linear in its coefficients and is solved by unweighted least squares.
+
+Insufficient variance produces empty spectrum arrays. Insufficient or
+unidentified fits produce explicit statuses and unavailable coefficients, not
+climatological substitutions. Field statistics and spectra reject missing input
+pixels; callers must apply a deliberate missing-data policy before using these
+complete-grid estimators. Real inputs are converted to float64 for estimation
+without being modified. Neither function clips fitted parameters.
+
+Temporal estimation and modelling
+---------------------------------
+
+:func:`pysteps.timeseries.parametric.lagrangian_autocorrelation` accepts prepared
+cascade histories with shape ``(time, level, y, x)`` and a fixed velocity array
+with shape ``(2, y, x)`` in pixels per timestep. It uses successive unit-step
+advection. Its default ``boundary="zero"`` fills values outside the grid with
+zero after each step and correlates the whole grid. Zero is the mean anomaly
+of a normalized cascade. The ``boundary="mask"`` alternative excludes pixels
+that are outside the valid support of any aligned history. These two boundary
+choices produce different correlation estimates.
+
+:func:`pysteps.timeseries.parametric.ar2_correlation_time` converts observed lag
+pairs into the stationary AR(2) one-sided discrete ACF sum including lag zero.
+Output time units equal the supplied timestep units. Invalid and nonstationary
+pairs remain unavailable. This quantity is not a forecast-skill horizon.
+
+Reference-time estimation rescales the selected level times to one reference
+wavelength and takes their geometric mean. By default it uses the first three
+levels; every contributing time must be positive and finite. Use ``levels=(0,)``
+to estimate from the largest level alone, as the forecast does by default.
+Supply the scaling exponent separately. Use nominal wavelengths from the
+cascade filter::
+
+    from pysteps.cascade.bandpass_filters import filter_gaussian
+    from pysteps.timeseries.parametric import (
+        estimate_reference_correlation_time, model_cascade_correlations
+    )
+
+    bandpass = filter_gaussian((256, 256), 6, d=1 / pixel_size_km)
+    wavelengths_km = 1 / bandpass["central_freqs"]
+    observed_times_minutes = np.array([180.0, 70.0, 40.0, 20.0, 10.0, 5.0])
+    exponent = 0.8  # illustrative temporal scaling exponent
+    reference_minutes = estimate_reference_correlation_time(
+        observed_times_minutes, wavelengths_km, exponent
+    )
+    correlations, diagnostics = model_cascade_correlations(
+        reference_minutes, wavelengths_km, exponent,
+        alpha=2.72, timestep=5.0, return_diagnostics=True
+    )
+
+The first band's nominal wavelength is a label for a broad low-pass response,
+not an energy-weighted effective wavelength. The temporal model applies the
+nonperiodic lag-2 adjustment when inverting the target timescale. Requested times
+at or below one timestep use white-noise lag pairs, achieving one timestep.
+Diagnostics retain requested and achieved times and identify sub-timestep limits.
+
+The estimators leave regional parameter selection and bounds to the caller.
+See :doc:`parametric_nowcast` for their use within a complete forecast and
+:doc:`calibrated_generator` for generating fields from the estimated parameters.