---
# System prepended metadata

title: pysteps blending

---

# pysteps blending

Steps structure:

```
check_inputs()

# FFT, extrapolation, noise and initialize
get_methods()

# Get perturbation generator and compute stddev adjustment
adj=get_noise_adj()

# Decompose input precipitation fields, normalize and rearrange for AR
decompose_cascade()

# Compute autocorrelation  GAMMA and AR coefficients PHI
get_autocorr()

# Apply velocity perturbations
perturb_velocities()

compute_mask()

#CALCULATE NOWCAST
for t in timesteps:
    for i in ensemble:
        create_cascade(generate_noise())*adj
        iterate_AR()
        
        recompose_cascade()
        
        # Adjust CDF of forecast to match most recently observed precip field
        probability_matching()
        
        # Advect recomposed precipitation field
        extrapolate()
        
        
```

STEPS blended structure:

According to Seed 2013: 

Radar:
- Track radar (not needed explicitly since we provide velocity field)
- Decompose cascade
- Estimate parameters

(NWP: update skill)

Per ensemble member:
- Initialize noise cascade
- Perturb radar cascade
- Perturb radar track

Per forecast:
- NWP:
    - Track
    - Decompose cascade (will be done ahead of time)
    - Estimate parameters
- Blend
    - Tracks
    - Parameters (alpha, beta)
- Update
    - Radar cascade
    - Noise cascade
    - Blending weights
- Combine


```
check_inputs()

# FFT, extrapolation, noise and initialize
get_methods()

# Decompose input precipitation fields, normalize and rearrange for AR
decompose_cascade(R)

# Update radar parameters
Update_Radar() -- FROM CPP

# If we have NWP data, update the skill:
Update_Skill() -- FROM CPP

 Perform backward advection to compute QPE / accumulation 
Make_QPE()

### See code from merged_cascade.cpp:Make_Forecast()

# Read file with the recent and climatological skill and regress to mean
determine_skill()

# Get perturbation generator and compute stddev adjustment
# NB: changes might be required to allow alpha/beta to blend linearly to the 
# NWP values
adjR=get_noise_adj(R)

# Apply velocity perturbations
perturb_velocities()

# Note that NODATA are replaced by NWP?
# compute_mask()

#CALCULATE NOWCAST
for t in timesteps:
    decompose_cascade(NWP)
    for i in ensemble:
        # 1. Regress correlations towards climatology and update AR parameters
        regress_correlations()
        # 2. Calculate AR parameters using the Yule Walker equations
        update_phi()
        
        # 3. Calculate current weights for the forecast
        calculate_weights(correlations) # implemented in steps_blending
        
        # 4. Merge velocity fields (with correct weights)
        merge_velocities() #including smooth_velocities()
        # 5. Merge the mean and stddev (for normalized fields!)
        blend_means_sigmas()  # almost implemented in steps_blending
        
        # 6. Set parameters for stoch cascade and generate the next noise cascade
        # 7. Add model, radar and stoch cascades for the merged forecast.
        # NB calculate separate weights for missing and non-missing data.
        blend_cascades()  # implemented in steps_blending
        
        # 8. Renormalize cascade levels to have mean 0 and std 1
        normalize_cascade()
        
        # 9. Recompose the cascade back into a single rainfall field and
        # convert back to rainfall intensity, mm/h.
        recompose_cascade() # implemented in steps_blending
        
        # 10. Adjust CDF of forecast to match most recently observed precip field
        probability_matching()
        
        # 11. Add back orographic component (not done in pysteps)
        
        # 12. Normalize rainfall field mean (not done in pysteps)
        
        # 13. Advect recomposed precipitation field
        extrapolate()
        
        
```

NOTES:

* The parametric noise parameters alpha/beta need to be blended too: from the observed values to the one fitted to the NWP forecast data. NB: time frequency of NWP data will need to be high enough...
* 

