.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/linda_nowcasts.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_linda_nowcasts.py: LINDA nowcasts ============== This example shows how to compute and plot a deterministic and ensemble LINDA nowcasts using Swiss radar data. .. GENERATED FROM PYTHON SOURCE LINES 10-24 .. code-block:: Python from datetime import datetime import warnings warnings.simplefilter("ignore") import matplotlib.pyplot as plt from pysteps import io, rcparams from pysteps.motion.lucaskanade import dense_lucaskanade from pysteps.nowcasts import linda, sprog, steps from pysteps.utils import conversion, dimension, transformation from pysteps.visualization import plot_precip_field .. GENERATED FROM PYTHON SOURCE LINES 25-27 Read the input rain rate fields ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 27-62 .. code-block:: Python date = datetime.strptime("201701311200", "%Y%m%d%H%M") data_source = "mch" # Read the data source information from rcparams datasource_params = rcparams.data_sources[data_source] # Find the radar files in the archive fns = io.find_by_date( date, datasource_params["root_path"], datasource_params["path_fmt"], datasource_params["fn_pattern"], datasource_params["fn_ext"], datasource_params["timestep"], num_prev_files=2, ) # Read the data from the archive importer = io.get_method(datasource_params["importer"], "importer") reflectivity, _, metadata = io.read_timeseries( fns, importer, **datasource_params["importer_kwargs"] ) # Convert reflectivity to rain rate rainrate, metadata = conversion.to_rainrate(reflectivity, metadata) # Upscale data to 2 km to reduce computation time rainrate, metadata = dimension.aggregate_fields_space(rainrate, metadata, 2000) # Plot the most recent rain rate field plt.figure() plot_precip_field(rainrate[-1, :, :]) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_linda_nowcasts_001.png :alt: linda nowcasts :srcset: /auto_examples/images/sphx_glr_linda_nowcasts_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 63-65 Estimate the advection field ---------------------------- .. GENERATED FROM PYTHON SOURCE LINES 65-69 .. code-block:: Python # The advection field is estimated using the Lucas-Kanade optical flow advection = dense_lucaskanade(rainrate, verbose=True) .. rst-class:: sphx-glr-script-out .. code-block:: none Computing the motion field with the Lucas-Kanade method. --- 5 outliers detected --- --- LK found 176 sparse vectors --- --- 63 samples left after declustering --- --- total time: 0.33 seconds --- .. GENERATED FROM PYTHON SOURCE LINES 70-72 Deterministic nowcast --------------------- .. GENERATED FROM PYTHON SOURCE LINES 72-118 .. code-block:: Python # Compute 30-minute LINDA nowcast with 8 parallel workers # Restrict the number of features to 15 to reduce computation time nowcast_linda = linda.forecast( rainrate, advection, 6, max_num_features=15, add_perturbations=False, num_workers=8, measure_time=True, )[0] # Compute S-PROG nowcast for comparison rainrate_db, _ = transformation.dB_transform( rainrate, metadata, threshold=0.1, zerovalue=-15.0 ) nowcast_sprog = sprog.forecast( rainrate_db[-3:, :, :], advection, 6, n_cascade_levels=6, R_thr=-10.0, ) # Convert reflectivity nowcast to rain rate nowcast_sprog = transformation.dB_transform( nowcast_sprog, threshold=-10.0, inverse=True )[0] # Plot the nowcasts fig = plt.figure(figsize=(9, 4)) ax = fig.add_subplot(1, 2, 1) plot_precip_field( nowcast_linda[-1, :, :], title="LINDA (+ 30 min)", ) ax = fig.add_subplot(1, 2, 2) plot_precip_field( nowcast_sprog[-1, :, :], title="S-PROG (+ 30 min)", ) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_linda_nowcasts_002.png :alt: LINDA (+ 30 min), S-PROG (+ 30 min) :srcset: /auto_examples/images/sphx_glr_linda_nowcasts_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Computing LINDA nowcast ----------------------- Inputs ------ dimensions: 320x355 number of time steps: 3 Methods ------- nowcast type: deterministic feature detector: blob extrapolator: semilagrangian kernel type: anisotropic Parameters ---------- number of time steps: 6 ARI model order: 1 localization window radius: 64.0 Detecting features... found 15 blobs in 0.75 seconds. Transforming to Lagrangian coordinates... 0.11 seconds. Estimating the first convolution kernel... 12.03 seconds. Estimating the ARI(p,1) parameters... 0.18 seconds. Estimating the second convolution kernel... 9.53 seconds. Computing nowcast for time step 1... 3.18 seconds. Computing nowcast for time step 2... 3.31 seconds. Computing nowcast for time step 3... 3.30 seconds. Computing nowcast for time step 4... 3.31 seconds. Computing nowcast for time step 5... 3.33 seconds. Computing nowcast for time step 6... 3.30 seconds. Computing S-PROG nowcast ------------------------ Inputs ------ input dimensions: 320x355 Methods ------- extrapolation: semilagrangian bandpass filter: gaussian decomposition: fft conditional statistics: no probability matching: cdf FFT method: numpy domain: spatial Parameters ---------- number of time steps: 6 parallel threads: 1 number of cascade levels: 6 order of the AR(p) model: 2 precip. intensity threshold: -10.0 ************************************************ * Correlation coefficients for cascade levels: * ************************************************ ----------------------------------------- | Level | Lag-1 | Lag-2 | ----------------------------------------- | 1 | 0.999004 | 0.996024 | ----------------------------------------- | 2 | 0.996659 | 0.984401 | ----------------------------------------- | 3 | 0.988483 | 0.960524 | ----------------------------------------- | 4 | 0.947887 | 0.851055 | ----------------------------------------- | 5 | 0.795798 | 0.597699 | ----------------------------------------- | 6 | 0.339528 | 0.180669 | ----------------------------------------- **************************************** * AR(p) parameters for cascade levels: * **************************************** ------------------------------------------------------ | Level | Phi-1 | Phi-2 | Phi-0 | ------------------------------------------------------ | 1 | 1.912659 | -0.914566 | 0.018047 | ------------------------------------------------------ | 2 | 1.842805 | -0.848982 | 0.043159 | ------------------------------------------------------ | 3 | 1.703894 | -0.723747 | 0.104429 | ------------------------------------------------------ | 4 | 1.390818 | -0.467283 | 0.281683 | ------------------------------------------------------ | 5 | 0.873046 | -0.097069 | 0.602702 | ------------------------------------------------------ | 6 | 0.314433 | 0.073910 | 0.938023 | ------------------------------------------------------ Starting nowcast computation. Computing nowcast for time step 1... done. Computing nowcast for time step 2... done. Computing nowcast for time step 3... done. Computing nowcast for time step 4... done. Computing nowcast for time step 5... done. Computing nowcast for time step 6... done. .. GENERATED FROM PYTHON SOURCE LINES 119-124 The above figure shows that the filtering scheme implemented in LINDA preserves small-scale and band-shaped features better than S-PROG. This is because the former uses a localized elliptical convolution kernel instead of the cascade-based autoregressive process, where the parameters are estimated over the whole domain. .. GENERATED FROM PYTHON SOURCE LINES 126-128 Probabilistic nowcast --------------------- .. GENERATED FROM PYTHON SOURCE LINES 128-177 .. code-block:: Python # Compute 30-minute LINDA nowcast ensemble with 40 members and 8 parallel workers nowcast_linda = linda.forecast( rainrate, advection, 6, max_num_features=15, add_perturbations=True, vel_pert_method=None, num_ens_members=40, num_workers=8, measure_time=True, )[0] # Compute 40-member STEPS nowcast for comparison nowcast_steps = steps.forecast( rainrate_db[-3:, :, :], advection, 6, 40, n_cascade_levels=6, R_thr=-10.0, mask_method="incremental", kmperpixel=2.0, timestep=datasource_params["timestep"], vel_pert_method=None, ) # Convert reflectivity nowcast to rain rate nowcast_steps = transformation.dB_transform( nowcast_steps, threshold=-10.0, inverse=True )[0] # Plot two ensemble members of both nowcasts fig = plt.figure() for i in range(2): ax = fig.add_subplot(2, 2, i + 1) ax = plot_precip_field( nowcast_linda[i, -1, :, :], geodata=metadata, colorbar=False, axis="off" ) ax.set_title(f"LINDA Member {i+1}") for i in range(2): ax = fig.add_subplot(2, 2, 3 + i) ax = plot_precip_field( nowcast_steps[i, -1, :, :], geodata=metadata, colorbar=False, axis="off" ) ax.set_title(f"STEPS Member {i+1}") .. image-sg:: /auto_examples/images/sphx_glr_linda_nowcasts_003.png :alt: LINDA Member 1, LINDA Member 2, STEPS Member 1, STEPS Member 2 :srcset: /auto_examples/images/sphx_glr_linda_nowcasts_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Computing LINDA nowcast ----------------------- Inputs ------ dimensions: 320x355 number of time steps: 3 Methods ------- nowcast type: ensemble feature detector: blob extrapolator: semilagrangian kernel type: anisotropic Parameters ---------- number of time steps: 6 ARI model order: 1 localization window radius: 64.0 error dist. window radius: 48.0 error ACF window radius: 80.0 ensemble size: 40 parallel workers: 8 seed: None Detecting features... found 15 blobs in 0.74 seconds. Transforming to Lagrangian coordinates... 0.09 seconds. Estimating the first convolution kernel... 11.63 seconds. Estimating the ARI(p,1) parameters... 0.18 seconds. Estimating the second convolution kernel... 9.27 seconds. Estimating forecast errors... Computing nowcast for time step 1... done. 0.79 seconds. Estimating perturbation parameters... 40.97 seconds. Computing nowcast for time step 1... 19.03 seconds. Computing nowcast for time step 2... 19.51 seconds. Computing nowcast for time step 3... 19.44 seconds. Computing nowcast for time step 4... 19.50 seconds. Computing nowcast for time step 5... 19.51 seconds. Computing nowcast for time step 6... 19.60 seconds. Computing STEPS nowcast ----------------------- Inputs ------ input dimensions: 320x355 km/pixel: 2.0 time step: 5 minutes Methods ------- extrapolation: semilagrangian bandpass filter: gaussian decomposition: fft noise generator: nonparametric noise adjustment: no velocity perturbator: None conditional statistics: no precip. mask method: incremental probability matching: cdf FFT method: numpy domain: spatial Parameters ---------- number of time steps: 6 ensemble size: 40 parallel threads: 1 number of cascade levels: 6 order of the AR(p) model: 2 precip. intensity threshold: -10.0 ************************************************ * Correlation coefficients for cascade levels: * ************************************************ ----------------------------------------- | Level | Lag-1 | Lag-2 | ----------------------------------------- | 1 | 0.999004 | 0.996024 | ----------------------------------------- | 2 | 0.996659 | 0.984401 | ----------------------------------------- | 3 | 0.988483 | 0.960524 | ----------------------------------------- | 4 | 0.947887 | 0.851055 | ----------------------------------------- | 5 | 0.795798 | 0.597699 | ----------------------------------------- | 6 | 0.339528 | 0.180669 | ----------------------------------------- **************************************** * AR(p) parameters for cascade levels: * **************************************** ------------------------------------------------------ | Level | Phi-1 | Phi-2 | Phi-0 | ------------------------------------------------------ | 1 | 1.912659 | -0.914566 | 0.018047 | ------------------------------------------------------ | 2 | 1.842805 | -0.848982 | 0.043159 | ------------------------------------------------------ | 3 | 1.703894 | -0.723747 | 0.104429 | ------------------------------------------------------ | 4 | 1.390818 | -0.467283 | 0.281683 | ------------------------------------------------------ | 5 | 0.873046 | -0.097069 | 0.602702 | ------------------------------------------------------ | 6 | 0.314433 | 0.073910 | 0.938023 | ------------------------------------------------------ Starting nowcast computation. Computing nowcast for time step 1... done. Computing nowcast for time step 2... done. Computing nowcast for time step 3... done. Computing nowcast for time step 4... done. Computing nowcast for time step 5... done. Computing nowcast for time step 6... done. .. GENERATED FROM PYTHON SOURCE LINES 178-184 The above figure shows the main difference between LINDA and STEPS. In addition to the convolution kernel, another improvement in LINDA is a localized perturbation generator using the short-space Fourier transform (SSFT) and a spatially variable marginal distribution. As a result, the LINDA ensemble members preserve the anisotropic and small-scale structures considerably better than STEPS. .. GENERATED FROM PYTHON SOURCE LINES 184-187 .. code-block:: Python plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_linda_nowcasts_004.png :alt: linda nowcasts :srcset: /auto_examples/images/sphx_glr_linda_nowcasts_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (4 minutes 34.966 seconds) .. _sphx_glr_download_auto_examples_linda_nowcasts.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: linda_nowcasts.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: linda_nowcasts.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: linda_nowcasts.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_