Weather Radar Predictions Using Images with pySTEPS

Recent advances in computing and data processing have made it possible to predict short-term weather — nowcasting — by leveraging visual radar imagery. In this article, we will explore how to generate weather radar predictions using Python and the pySTEPS library. We’ll go through every step: from setting up your environment and understanding key techniques like Lucas–Kanade motion estimation and S-PROG nowcasting, to building an iterative script that updates its forecast using previous outputs.

In our examples, we use precipitation data from weather radars. Currently, the Lithuanian Hydrometeorological Service under the Ministry of Environment operates two weather radars in Lithuania: one in Laukuva (Šilalė district, western Lithuania) and another in Trakų Vokė (on the outskirts of Vilnius).

These radars provide information every 5 minutes and cover an area of approximately 150 kilometers each. Throughout the article, you’ll find comparisons of the observed data with forecasts at 5 minutes, 30 minutes, 1 hour, and 2 hours into the future. Few example sets are included to demonstrate these comparisons.

Why?

Imagine being able to predict rainfall for the next 30 minutes or even several hours by analyzing just a series of radar images. This is the promise of nowcasting — a method to forecast very short-term weather by extrapolating current observations. Our approach combines two essential techniques:

• Lucas–Kanade Motion Estimation: A classic optical flow method that calculates the movement of features (in our case, precipitation patterns) between consecutive radar images.
• S-PROG Nowcasting: A statistical technique that decomposes the precipitation field into multiple spatial scales (cascade levels) and uses an autoregressive (AR) model to generate forecasts.

Python script that we will analyse in this article is designed to work iteratively: each forecast is generated using previous observations, and in some configurations, even previous forecasts, allowing the model to “learn” from its own output. While the basic approach primarily warps (or shifts) existing features, thoughtful parameter tuning can promote a more generative output — where new details emerge (or old details disappear).

In this article, we break down the code, explain the underlying principles, and provide detailed instructions for setting up and running your own nowcasting system. Whether you’re a seasoned Python programmer or a newcomer to meteorological data processing, you’ll find our explanations and code examples accessible and practical.

While our current approach relies on classical methods, it also offers a solid foundation for integrating advanced techniques later that could further enhance forecast detail and accuracy.

Prerequisites

Before you begin, ensure your environment is ready with the following tools and libraries:

• Python 3.x
• NumPy: For numerical operations.
• Pillow (PIL): For image processing.
• Requests: For fetching data from remote APIs.
• Matplotlib: For visualization and color mapping.
• pySTEPS: A specialized library for precipitation nowcasting.

For the dataset, it is recommended to start with at least 20 images for nowcasting.

Installing Dependencies using pip

• Ensure Python is Installed:
  Make sure you have Python 3.x installed on your system. You can verify by running:

python --version

• Install Required Packages:
  Once you know that Python is present on your system, install the packages with pip:

pip install numpy pillow requests matplotlib pysteps

• Verifying the Installation:
  To check that all packages were installed correctly, you can run:

pip list

This will display a list of installed packages and their versions.

Installing Dependencies using Conda

As an alternative, you can use Conda package manager. Below is an explanation on how to install the required dependencies using Conda.

• Setting Up a Conda Environment:
  It’s a good practice to create a new Conda environment for your project. Managing dependencies in an isolated environment reduces the chance of conflicts with other projects on your system. This is especially useful when working on scientific projects where package versions can be critical. For example, you can create one named weather_env using:

conda create --name weather_env python=3.9
conda activate weather_env

Feel free to change the Python version if needed.

• Installing Dependencies with Conda:
  Once you have activated your new environment, you can install most of the required packages from the conda-forge channel. Run:

conda install -c conda-forge numpy pillow requests matplotlib pysteps

If, for any reason, a package isn’t available on conda-forge, you can install it via pip after installing the rest with conda. For example:

conda install -c conda-forge numpy pillow requests matplotlib
pip install pysteps

This approach combines Conda’s robust environment management with pip’s flexibility.

Access to Images

Additionally, you’ll need access to radar images. If you have your own dataset, place your images in an input folder. When using the script, remember to update the folder paths in the code from: "/your/path/to/input/images" to your local directory.

Otherwise, if you choose remote way you’ll need to set an URL to fetch data from. Do not forget to update the decoding part of script too. It will be outlined later in the code snippet.

An example of one of the weather radar images (out of 20+ sequenced imagery) for my use case:

Understanding Key Concepts

Lucas–Kanade Motion Estimation

One of the main points of our approach is the Lucas–Kanade method — a popular algorithm in computer vision for estimating optical flow. In simple terms, it determines how different parts of an image shift from one frame to the next. For weather radar, this means calculating how rainfall patterns move over time. The output of Lucas–Kanade is a velocity field that provides the speed and direction of movement at each pixel. This information is crucial for nowcasting because it allows us to extrapolate (or “warp”) the current precipitation field into the near future.

S-PROG Nowcasting

S-PROG stands for Statistical Precipitation Nowcasting by a Stochastic Realization Of a Generative model. It is a specialized method for short-term precipitation forecasting. The process involves:

• Extrapolation: Moving the current precipitation pattern forward in time using the velocity field obtained from Lucas–Kanade.
• Autoregressive Generation: Decomposing the precipitation field into multiple cascade levels and applying an AR model to generate new details.

An AR model (autoregressive model) is a type of statistical model used to describe and forecast time series data. In an AR model, the current value of a variable is expressed as a linear combination of its previous values plus a random error term. In the context of precipitation nowcasting with S-PROG, the AR model helps capture the temporal evolution of the precipitation field by using past radar images (or their decomposed components) to predict future states.

For example, in an AR(1) model (first-order autoregressive), the current forecast is based on the immediately preceding value; an AR(2) model would use the two previous values, and so on. This autoregressive process allows the model to “learn” patterns and trends from past observations, which is critical for generating plausible forecasts.

By integrating an AR model into S-PROG, the method not only shifts the observed precipitation field (using the velocity field from Lucas–Kanade) but also introduces variability based on the historical behavior of the precipitation, potentially generating new details in the forecast.

This combination of motion-based extrapolation and autoregressive generation is what gives S-PROG its ability to produce forecasts that evolve over time, blending both the warping of existing patterns and the generation of new details.

S-PROG primarily uses the observed patterns and their motion to forecast future precipitation. However, by iteratively updating its inputs, the method can begin to generate new structures rather than merely shifting the old ones.

Code Overview

Our code is structured in a logical sequence that mirrors the nowcasting process:

1. Output Folder Cleanup:
   The script starts by clearing the output folder to ensure no imagery of previous runs interfere with the current forecast.
2. Data Acquisition (Remote vs. Local):
   The code can either download radar images from a remote API (in our example, the images are encoded in Base64) or load images from a local folder.
3. Reading and Processing Input Images:
   Each image is read and converted into two formats: RGBA for final visualization and colorization and grayscale (L) for computing motion and feeding into the nowcasting model.
4. Motion Estimation:
   Using a selected number of recent grayscale images, the script computes a velocity field with the Lucas–Kanade algorithm.
5. Preparing Input for S-PROG:
   The last few frames (determined by the autoregressive order) are used as input to the S-PROG nowcasting model, providing the necessary temporal context.
6. Running the S-PROG Forecast:
   The S-PROG method is applied to generate forecast frames based on the input images and the velocity field.
7. Fixed Normalization and Custom Color Mapping:
   The forecast frames are normalized using a global intensity range, and a custom colormap is applied to convert numerical intensities into colorful, visually interpretable images.
8. Iterative Forecast Updating:
   The script updates its input iteratively by blending the new forecast with the previous input, allowing the model to generate more new detail instead of just warping existing patterns.
9. Saving Forecast Frames:
   Finally, each forecast frame is saved as a PNG file, and debug information is printed to help you understand the processing at each step.

Throughout the article, there will be detailed explanations of each step included, and I’ve also provided images comparing observed data with forecasts at different time intervals.

Detailed Step-by-Step Explanation

Let’s analyse how our script is going to look like in detail:

1. Output Folder Cleanup
Before starting any processing, it’s important to clear the output folder. The function clear_output_folder(folder)deletes all files in the specified directory. This ensures that any forecasts from previous runs are removed, and you begin with a clean start.

2. Data Acquisition: Remote vs. Local
The script supports two modes:

• Remote Mode:
  When you run the script in remote mode, it contacts a specified API endpoint. In our example, the API returns radar data as a JSON object with Base64-encoded images. The function fetch_radar_json() retrieves this data, and decode_and_save_radar_images() decodes and saves each image as a PNG in the input folder. Make sure to update the decode_and_save_radar_images() function for the remote fetch to correctly decode the source.
• Local Mode:
  If you already have radar images stored locally (e.g., downloaded from the Internet), simply place them in your specified input folder. The script will then use these images directly.

3. Reading and Processing Input Images
Once your images are in the input folder, the script reads them one by one. Each image is processed into two different formats:

• RGBA Format:
  This format retains color and transparency information and is used for the final output visualization.
• Grayscale (L) Format:
  The grayscale version reduces the complexity of the data, making it easier to compute the velocity field and run the nowcasting algorithm.

4. Motion Estimation
Using the last NUM_FRAMES_FOR_MOTION images (typically around 20), the script computes a velocity field with the Lucas–Kanade algorithm via lk.dense_lucaskanade(frames_for_motion). This velocity field represents how the precipitation features move across frames and is essential for predicting the future state of the precipitation field.

5. Preparing Input for S-PROG
The S-PROG nowcasting model requires a sequence of recent images to provide temporal context. The script selects the last AR_ORDER+1 frames from the grayscale image stack. For example, if AR_ORDER is set to 1, it uses the last 2 frames as input.

6. Running the S-PROG Forecast
The core function run_sprog_forecast() takes the input frames and the velocity field, and runs the S-PROG model to generate a forecast for a single lead time. The S-PROG model uses parameters such as:

• PRECIP_THR: The minimum precipitation intensity threshold.
• EXTRAP_METHOD: The method used for extrapolation (set to "semilagrangian").
• SPROG_KWARGS: A dictionary of additional parameters controlling the cascade decomposition and AR model.

Any invalid values in the forecast output (e.g., NaNs) are replaced with zeros to ensure clean data.

7. Fixed Normalization and Custom Color Mapping
After forecasting, the script calculates the global minimum and maximum intensity values across all forecast frames. These values define a fixed normalization range so that every forecast frame is mapped consistently. A custom colormap (defined by the list custom_colors) then maps these normalized values to specific colors. The colorize_forecast_field() function applies normalization, gamma correction (controlled by GAMMA_VALUE), and then uses the custom colormap to generate an 8-bit RGBA image.

8. Iterative Forecast Updating
One of the key features of our approach is its iterative nature. Instead of forecasting each frame independently, the script updates its input sequence for each lead time by blending the new forecast with the previous last frame. This is done using a weighted update:

blended_last = GEN_WEIGHT * forecast_field + (1 - GEN_WEIGHT) * current_input[-1]

If you set GEN_WEIGHT to 1, it means that the new forecast fully replaces the last frame. You could reduce have this value reduced (e.g., to 0.7) if you’d like to retain some information in next image from the previous one. This iterative blending helps the system generate new features rather than simply shifting the observed pattern.

9. Saving Forecast Frames
Once all forecast frames have been generated, the script normalizes, colorizes, and saves each one as a PNG file in the output folder. For each frame, it also creates a dynamic alpha mask: pixels with an intensity above ALPHA_THRESH are rendered fully opaque, while the rest become transparent.

Debug messages print useful information about the intensity ranges, which can be helpful for further parameter tuning.

Example Visual Comparisons

Here, you can find visual comparisons that illustrate the forecasting evolution. It’ll help understand how the nowcasting model extrapolates precipitation patterns and how well the forecasts match observed conditions (accuracy).

Dataset #1

• Last fact — The most recent observed radar imagery.

• 5 minutes into Forecast — A radar image captured 5 minutes after the last observation compared with the first forecasted image.

• 30 minutes into Forecast — Observed and forecasted images 30 minutes into.

• 1 hour into Forecast — Observed and forecasted images 1 hour into.

• 2 hours into Forecast — Observed and forecasted images 2 hours into.

Dataset #2

• Last fact — The most recent observed radar imagery.

• 5 minutes into Forecast — A radar image captured 5 minutes after the last observation compared with the first forecasted image.

• 30 minutes into Forecast — Observed and forecasted images 30 minutes into.

• 1 hour into Forecast — Observed and forecasted images 1 hour into.

• 2 hours into Forecast — Observed and forecasted images 2 hours into.

Dataset #3

• Last fact — The most recent observed radar imagery.

• 5 minutes into Forecast — A radar image captured 5 minutes after the last observation compared with the first forecasted image.

• 30 minutes into Forecast — Observed and forecasted images 30 minutes into.

• 1 hour into Forecast — Observed and forecasted images 1 hour into.

• 2 hours into Forecast — Observed and forecasted images 2 hours into.

We can notice that it sometimes fails to retain the correct color and puts more intense spectrum where the precipitation is light and vice versa. While this is not perfect, the problems can be solved by further improving the algorithm. I’ll leave it to your creative mind! :)

Solution

Modifiable Variables and Their Limitations

Before you use the full code, it’s important to understand the key modifiable variables and their practical limitations:

• INPUT_FOLDER and OUTPUT_FOLDER — Paths to directories where input images are stored and output images will be saved.
  Limitation: Ensure the folders exist and that you have proper read/write permissions.
• NUM_FRAMES_FOR_MOTION — The number of recent frames used for motion estimation (e.g., 4).
  Limitation: Too few frames may lead to noisy velocity estimates; too many may smooth out fine-scale motion.
• N_LEAD_TIMES — The total number of forecast frames to generate.
  Limitation: A higher number can provide a longer forecast but may increase error accumulation.
• GEN_WEIGHT — A weight between 0 and 1 determining how much of the new forecast replaces the previous frame in the input update.
  Limitation: A value of 1 means full replacement (pure warping). Lower values (e.g., 0.7) blend in part of the previous frame, which may help introduce new details but might also retain outdated information.
• AR_ORDER — The order of the autoregressive model used in S-PROG (e.g., 1 means two input frames are used).
  Limitation: Higher values (like 6 or 8) are not practical (and possibly throw an error), as the model is typically designed for low AR orders (often 1–3). Too high an order might lead to overfitting or unstable forecasts.
• PRECIP_THR — The precipitation threshold; values below this are treated as no-rain.
  Limitation: If set too high, subtle precipitation signals may be ignored; if too low, noise may be amplified.
• EXTRAP_METHOD — The method for extrapolation, set here to "semilagrangian".
  Limitation: This choice impacts how the forecast field is shifted forward in time. Other methods may be available but could require additional tuning.
• SPROG_KWARGS — A dictionary containing additional parameters for S-PROG, such as the number of cascade levels (n_cascade_levels), the bandpass filter method, and the domain (spatial or spectral).
  Limitation: For n_cascade_levels a value between 2 and 8 should be used. Too few levels may not capture all spatial details, while too many levels can result in over-smoothing.
• domain: The spatial domain preserves more detail but may be more sensitive to noise. The spectral domain can yield smoother forecasts but may lose fine detail.
• MARGIN — An optional margin used during normalization.
  Limitation: Usually set to 0, but adjusting it can affect how aggressively the data is clipped during visualization.
• ALPHA_THRESH — The intensity threshold used to create a dynamic alpha mask for visualization.
  Limitation: This parameter controls which pixels are considered significant (opaque) in the final image. Adjusting it too high may result in fewer visible details, while too low may produce overly bright images.
• GAMMA_VALUE — Gamma correction factor applied during colorization.
  Limitation: Gamma values significantly alter the color mapping — a value of 1 means no gamma correction. Values greater than 1 shift lower intensities further down the color scale (making them appear darker or more “bluish” if your colormap starts with blue). Values lower than 1 brighten lower intensities, possibly pushing them into mid-range colors.
• n_lead — The number of forecast frames produced in a single call to S-PROG.
  Limitation: Setting n_lead to a higher value (e.g., 24) can generate a longer forecast horizon at once, but may also accumulate errors more quickly. Using a smaller value (like 1) and calling the function iteratively can yield more frequent updates or blending, yet risks AR-model instability if the frames become too similar over time.
• Custom Colormap (CUSTOM_CMAP) — A colormap created using a list of hex color codes. In this example, the colors progress from bluish to yellow.
  Limitation: You can modify these hex codes to suit your visual preference. Ensure that the chosen colors provide a meaningful gradient for your data.
  You can also use the built-in colormap from Matplotlib instead of a custom one. To do this, you can simply set your colormap variable to the nipy_spectral, turbo, rainbow or another colormap. For example, you can replace the custom colormap definition with:

import matplotlib.cm as cm

CUSTOM_CMAP = cm.get_cmap("nipy_spectral")

Source code

Below is the complete code used in this article. Make sure to adjust the folder paths, parameters and remote data decoding algorithm to suit your needs.

⚠️ A basic understanding of tools mentioned in the article is helpful when modifying or troubleshooting the code as it is probably not going to work from the first time. Even better, ChatGPT can help a lot and solve the issues you may face in seconds. Really.

#!/usr/bin/env python3
import os
import shutil
import json
import requests
import base64
import io
import argparse
import numpy as np
from PIL import Image

# PySTEPS imports: used for motion estimation and nowcasting.
import pysteps
import pysteps.motion.lucaskanade as lk
import pysteps.nowcasts.sprog as sprog

# Matplotlib for color mapping.
import matplotlib
import matplotlib.colors as mcolors
import matplotlib.colors as colors

##############################################################################
# CONFIGURATION
##############################################################################
# Radar API endpoint.
RADAR_ENDPOINT = "https://api.weather.com/radar-imagery"

# Set your own paths for input and output directories.
INPUT_FOLDER = "/your/path/to/input/images"     # Replace with your input folder path.
OUTPUT_FOLDER = "/your/path/to/output/images"     # Replace with your output folder path.

def clear_output_folder(folder):
    """Deletes all files in the specified folder to ensure a clean run."""
    if os.path.isdir(folder):
        for f in os.listdir(folder):
            path = os.path.join(folder, f)
            if os.path.isfile(path):
                os.remove(path)

os.makedirs(INPUT_FOLDER, exist_ok=True)
os.makedirs(OUTPUT_FOLDER, exist_ok=True)
clear_output_folder(OUTPUT_FOLDER)

# Forecast parameters.
NUM_FRAMES_FOR_MOTION = 4        # Number of frames used for motion estimation.
N_LEAD_TIMES = 24                # Total number of forecast frames to generate.

# Generation weight: the fraction of the new forecast that is used when updating input.
# GEN_WEIGHT = 1 means the new forecast fully replaces the last frame.
GEN_WEIGHT = 1

# AR order and precipitation threshold.
AR_ORDER = 1                    # Order of the autoregressive model (practically, 1 to 3 are common)
PRECIP_THR = 80                 # Minimum precipitation threshold for forecast.
EXTRAP_METHOD = "semilagrangian"

# S-PROG configuration: parameters for cascade decomposition and AR model.
SPROG_KWARGS = {
    "conditional": False,
    "probmatching_method": None,
    "n_cascade_levels": 3,      # Number of cascade levels (tune between 3 and 8; higher may oversmooth).
    "bandpass_filter_method": "gaussian",
    "domain": "spatial",        # Choose "spatial" or "spectral".
    "norain_thr": 0.0           # Minimum rain intensity below which the data is considered as 'no rain'
}

MARGIN = 0                     # Additional margin for normalization.
ALPHA_THRESH = 10              # Alpha threshold for visualization.

# Custom colormap: define colors from low to high intensity.
custom_colors = [
    "#5AC4C5",  # bluish (low intensity)
    "#429443",  # green
    "#6AE446",  # brighter green
    "#79FB54",  # greenish-yellow
    "#FFFF5A"   # yellow (high intensity)
]
CUSTOM_CMAP = mcolors.LinearSegmentedColormap.from_list("custom_cmap", custom_colors)

GAMMA_VALUE = 9  # Gamma correction factor for colorization.

##############################################################################
def fetch_radar_json():
    """Fetches radar data from the API endpoint."""
    resp = requests.get(RADAR_ENDPOINT)
    resp.raise_for_status()
    return resp.json()

# Update decoding algorithm appropariate to your own data decoding
def decode_and_save_radar_images(radar_list):
    """
    Decodes base64-encoded radar images from JSON and saves them as PNG files in INPUT_FOLDER.
    
    Returns:
        A list of saved filenames.
    """
    saved_filenames = []
    for entry in radar_list:
        date_str = entry["date"]
        b64_str = entry["base64"]
        raw_bytes = base64.b64decode(b64_str)
        pil_img = Image.open(io.BytesIO(raw_bytes))
        safe_date = date_str.replace(":", "").replace(" ", "_")
        filename = f"radar_{safe_date}.png"
        path = os.path.join(INPUT_FOLDER, filename)
        pil_img.save(path)
        saved_filenames.append(filename)
    return saved_filenames

##############################################################################
# Forecast function: generates one forecast frame using S-PROG.
##############################################################################
def run_sprog_forecast(precip_input, velocity_field, n_lead=1):
    """
    Generates a forecast using S-PROG.
    
    Parameters:
        precip_input: np.array, input precipitation frames.
        velocity_field: Motion field computed via Lucas–Kanade.
        n_lead: Number of lead times to forecast (typically 1 for iterative updates).
    
    Returns:
        Forecast fields as a numpy array of shape (n_lead, H, W).
    """
    precip_used = np.nan_to_num(precip_input, nan=0.0, posinf=0.0, neginf=0.0)
    forecast_raw = sprog.forecast(
        precip_used,
        velocity_field,
        n_lead,
        ar_order=AR_ORDER,
        precip_thr=PRECIP_THR,
        extrap_method=EXTRAP_METHOD,
        **SPROG_KWARGS
    )
    forecast_fields = np.nan_to_num(forecast_raw, nan=0.0, posinf=0.0, neginf=0.0)
    return forecast_fields

def colorize_forecast_field(frame, cmap=CUSTOM_CMAP, fixed_vmin=None, fixed_vmax=None, gamma=1.0):
    """
    Converts a forecast frame (float array) into an RGBA image using a custom colormap.
    
    Parameters:
        frame: np.array, the forecast field.
        cmap: Custom colormap.
        fixed_vmin, fixed_vmax: Normalization bounds (if None, computed from the frame).
        gamma: Gamma correction factor.
    
    Returns:
        An RGBA image as a numpy array (dtype uint8).
    """
    if fixed_vmin is None:
        fixed_vmin = np.min(frame)
    if fixed_vmax is None:
        fixed_vmax = np.max(frame)
    clipped = np.clip(frame, fixed_vmin, fixed_vmax)
    norm = colors.Normalize(vmin=fixed_vmin, vmax=fixed_vmax)
    normalized = norm(clipped)
    normalized = normalized ** gamma
    rgba = cmap(normalized)
    return (rgba * 255).astype(np.uint8)

##############################################################################
def main(mode):
    # Data acquisition: remote or local.
    if mode == "remote":
        print("Mode: REMOTE => Downloading images from API.")
        data_list = fetch_radar_json()
        if not data_list:
            print("No data from API. Exiting.")
            return
        decode_and_save_radar_images(data_list)
    else:
        print("Mode: LOCAL => Using images in input folder.")

    # Read input images from INPUT_FOLDER.
    all_files = sorted([f for f in os.listdir(INPUT_FOLDER) if f.lower().endswith(".png")])
    if len(all_files) < 2:
        print("Need >= 2 frames for motion. Exiting.")
        return

    color_frames = []
    gray_frames = []
    for fname in all_files:
        path = os.path.join(INPUT_FOLDER, fname)
        with Image.open(path) as img:
            # Convert image to RGBA for final visualization.
            rgba_arr = np.array(img.convert("RGBA"), dtype=np.uint8)
            color_frames.append(rgba_arr)
            # Convert image to grayscale (L mode) for processing.
            gray_arr = np.array(img.convert("L"), dtype=float)
            gray_frames.append(gray_arr)

    color_frames = np.stack(color_frames, axis=0)
    gray_frames = np.stack(gray_frames, axis=0)

    # Estimate motion using the last NUM_FRAMES_FOR_MOTION frames.
    use_from = max(0, len(gray_frames) - NUM_FRAMES_FOR_MOTION)
    frames_for_motion = gray_frames[use_from:]
    velocity_field = lk.dense_lucaskanade(frames_for_motion)

    # Prepare initial input for S-PROG using the last AR_ORDER+1 frames.
    needed = AR_ORDER + 1
    if len(gray_frames) >= needed:
        current_input = gray_frames[-needed:].copy()
    else:
        current_input = np.repeat(gray_frames[-1][None, ...], needed, axis=0)

    forecast_frames = []
    # Iteratively forecast each lead time frame using a blended update.
    for i in range(N_LEAD_TIMES):
        # Generate forecast for one lead time (n_lead=1).
        forecast_field = run_sprog_forecast(current_input, velocity_field, n_lead=1)[0]
        forecast_frames.append(forecast_field)
        # Update input: blend the new forecast with the previous last frame.
        blended_last = GEN_WEIGHT * forecast_field + (1 - GEN_WEIGHT) * current_input[-1]
        current_input = np.concatenate([current_input[1:], blended_last[None, ...]], axis=0)

    forecast_frames = np.array(forecast_frames)  # Shape: (N_LEAD_TIMES, H, W)

    # Compute global intensity range for normalization.
    global_min = float(np.min(forecast_frames))
    global_max = float(np.max(forecast_frames))
    global_min = max(0, global_min - MARGIN)
    global_max = global_max + MARGIN
    print(f"[DEBUG] Global forecast range: {global_min:.1f} .. {global_max:.1f}")

    # Colorize and save each forecast frame.
    for i in range(N_LEAD_TIMES):
        intensities = forecast_frames[i]
        print(f"[DEBUG] Frame {i+1} raw: min={intensities.min()}, max={intensities.max()}")

        colored = colorize_forecast_field(
            intensities,
            cmap=CUSTOM_CMAP,
            fixed_vmin=global_min,
            fixed_vmax=global_max,
            gamma=GAMMA_VALUE
        )

        # Create a dynamic alpha mask: pixels above ALPHA_THRESH are fully opaque.
        alpha_mask = np.where(intensities > ALPHA_THRESH, 255, 0).astype(np.uint8)
        colored[..., 3] = alpha_mask
        zero_mask = (alpha_mask == 0)
        colored[zero_mask] = (0, 0, 0, 0)

        out_name = f"forecast_{i+1}.png"
        out_path = os.path.join(OUTPUT_FOLDER, out_name)
        Image.fromarray(colored, mode="RGBA").save(out_path)
        print(f"Saved {out_path}")

        # Debug information for the colorized output.
        mask = colored[..., 3] > 0
        if np.any(mask):
            rgb = colored[..., :3].astype(float)
            intensity_zone = np.mean(rgb, axis=2)
            zone_min = intensity_zone[mask].min()
            zone_max = intensity_zone[mask].max()
            print(f"[DEBUG] Forecast frame {i+1} colored zone: min intensity={zone_min:.1f}, max intensity={zone_max:.1f}")
        else:
            print(f"[DEBUG] Forecast frame {i+1} colored zone: no non-transparent pixels.")

if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description="Weather radar nowcasting using S-PROG with iterative, blended forecast updating to balance motion extrapolation and generative detail."
    )
    parser.add_argument("--mode", choices=["local", "remote"], default="local", help="Use local input folder or fetch from API.")
    args = parser.parse_args()
    main(args.mode)

You can run this Python code in several ways. Here are a few common methods:

• Command Line (Terminal):
  Open your terminal (or Command Prompt on Windows) and navigate to the folder where your script is saved. Then run:

python your_script_name.py --mode local

Replace your_script_name.py with the name of your file in the directory. You can also use --mode remote if you want to fetch data from the API. Again, for the remote fetch you will need to update your decoding script to match the source and decode it correctly.

• Integrated Development Environment (IDE):
  You can run the code from within an IDE such as VS Code, PyCharm, or Spyder. Open your script in the IDE and run it using the built-in run command.

Risks, Limitations, and Opportunities

• Data Quality:
  The reliability of motion estimation and nowcasting is highly dependent on the quality and resolution of your radar images. Noisy or low-resolution data can lead to inaccurate forecasts.
• Iterative Forecasting:
  The iterative update mechanism (where each forecast is blended into the input) can introduce error accumulation over longer lead times. Fine-tuning parameters like GEN_WEIGHT and N_LEAD_TIMES is crucial to maintain forecast stability.
• Parameter Tuning:
  Key parameters such as AR_ORDER, PRECIP_THR, GEN_WEIGHT, and GAMMA_VALUE require careful tuning based on your specific dataset and forecasting goals. Small adjustments can significantly alter the forecast output.
• Prerequisites and Library Dependencies:
  This script uses several Python libraries — NumPy, Pillow, Requests, Matplotlib, and pySTEPS. A basic understanding of these tools is helpful when modifying or troubleshooting the code. Even better, ChatGPT can help a lot and solve the issues you may face in seconds.
• Opportunities for AI Integration:
  While this tutorial is based on classical methods (Lucas–Kanade and S-PROG), recent advances in deep learning (e.g., convolutional LSTMs, GANs, transformers) offer exciting possibilities for nowcasting. AI models can learn spatiotemporal patterns directly from large radar datasets, potentially providing more generative and detailed forecasts. However, these approaches require substantial training data and computational resources.

Why Choose This Solution?

This nowcasting solution is attractive because it is:

• Cost-Effective: Built entirely with open-source libraries, it does not require expensive proprietary software or hardware.
• Fast: The combination of classical methods ensures that forecasts can be generated quickly, making it suitable for real-time applications.
• Modular and Customizable: With numerous parameters available for tuning, you can tailor the system to different datasets and forecasting requirements.
• Accessible: The solution leverages widely used Python libraries, making it approachable for researchers, meteorologists, and enthusiasts alike.
• Foundation for Further Development: This system provides a robust baseline that you can enhance with AI-based techniques in the future.

Where to Go Next?

If you’ve enjoyed this exploration of weather radar nowcasting using pySTEPS, here are some directions for further study and development:

• MWNet:
  It is a deep learning-based nowcasting model that has shown promising results in precipitation forecasting. Exploring MWNet can help you integrate AI into your nowcasting pipeline for potentially more generative and detailed forecasts.
• COTREC:
  There is also another model designed for nowcasting using a combination of physical and statistical methods. Investigating COTREC could offer insights into alternative approaches that balance extrapolation with new feature generation.
• rainymotion Dense:
  This model focuses on dense optical flow techniques to generate forecasts. It might be interesting to compare its performance against the Lucas–Kanade method used in this tutorial.
• pySTEPS ANVIL:
  ANVIL is a component within pySTEPS that further refines the nowcasting process. Experimenting with ANVIL might yield improved forecasts with fewer artifacts.
• AI-Based Nowcasting:
  The rapidly evolving field of AI-based nowcasting offers a host of possibilities. Techniques such as convolutional LSTMs, GANs, and transformers can be trained on large datasets to learn complex spatiotemporal patterns. Integrating these models into your pipeline could lead to more accurate and creative forecasts.

Precipitation prediction for 60 minutes from 24. June 2021 18:40 CEST. Images from left by rows: last observation (18:40), ground truth (19:40), MWNet v1.0, MWNet v0.2, COTREC, rainymotion Dense, pySTEPS ANVIL, pySTEPS SPROG. Imagery from Meteopress.

Each of these solutions requires additional learning and experimentation. You may need to gather more data, train deep learning models, or adapt existing research code. However, exploring these advanced methods can significantly enhance the accuracy and detail of your weather forecasts.

Final Words

In this extensive guide, we explored how to generate weather radar nowcasts using Python and the pySTEPS library. We began by discussing the underlying principles of Lucas–Kanade motion estimation and S-PROG nowcasting, and why these methods are effective for short-term weather forecasting. We then walked through a detailed code example, explaining each step — from data acquisition and image processing to motion estimation, forecast generation, normalization, color mapping, and iterative updating.

This solution is not only cheap and fast but also highly customizable, making it an excellent starting point for anyone interested in nowcasting. Although our current approach is based on classical methods, it lays a really strong foundation for integrating advanced techniques in the future to further improve forecast detail and accuracy.

Whether you’re interested in enhancing the generative detail of your forecasts or exploring cutting-edge AI approaches, there is much to discover in this field.

Thank you for reading and good luck forecasting weather! :)

More about me and my experience here.