pySTEPS / pysteps

Python framework for short-term ensemble prediction systems.
https://pysteps.github.io/
BSD 3-Clause "New" or "Revised" License
441 stars 160 forks source link

STEPS-blending: Losing small scale features in the first time steps and other issues #332

Closed mpvginde closed 3 months ago

mpvginde commented 1 year ago

Hi everyone,

In setting up pySTEPS for a blended nowcast with STEPS, I have noticed that during the first time steps some of the radar image small-scale features are immediately lost. This becomes particularly apparent when the blended STEPS nowcast is compared to the non-blended version of the STEPS nowcast.

An example over the Belgian radar domain can be seen here: blendVSnoblend_t0 Fig. 1: Identical starting conditions (i.e. radar field) for the blenden nowcast (left) and non-blended nowcast (right)

blendVSnoblend_t1 Fig. 2: The blended (left) and non-blended (right) nowcast after 1 time step of 5 minutes

It's clear that even after the first time step (Fig, 2), the blended nowcast field has lost some fine-scale structures w.r.t. the non-blended nowcast (i.e. the blended fields seem smoother). In addition, there also seem to be more broad structures with more intense rainfall in the blended nowcast.

When investigating this phenomenon I came across several issues/points of attention in the blending code.\ I already discussed some of them with @RubenImhoff, but I will try to explain them here too so everything is well documented. \ \ In order to explore the problem I plotted the time evolution of the Power (displayed as standard deviation [St. Dev.]) of the different cascade levels for the different fields.\ For blending nowcast there are 2 relevant fields per level:

Both are evolved separately using an AR(2) model\ In addition I also plotted the time evolution of their respective blending weights with which they are combined. I have deliberately omitted the NWP-component as its weight is typically very small during the first time steps.

For the non-blending nowcast there is only 1 relevant field as the noise is directly incorporated in the AR(2) process.

Finally, also included is a figure of the actual fields themselves at a specific time step. For the blended nowcast:

For the non-blended nowcast:

overview_l6_t1 Fig. 3: AR(2)-evolution overview for cascade level 6 at time step 1

Fig 3. shows this AR(2)-evolution overview for cascade level 6 (2nd smallest level). We notice some problems:

  1. The AR(2) evolution of the precipitation cascade without the noise term EPS is no longer stationary in the variance (something that is necessary according to Bowler et al. 2006: \ "The noise term has been omitted from this formulation and the cascade is renormalized at each step to ensure no loss of power") In practice, this means a dampening of the amplitudes (= loss of power) for all cascade levels
  2. The noise starts from 0 and needs to spin-up. For the largest levels (due to the large temporal autocorrelations) this may take a long time.

This results in the following issues:

What does this mean: For the smallest levels, due to the low temporal autocorrelation, the precipitation field decays very fast and also the blending weight immediately is very low. This mean that the original precipitation signal is almost immediately lost. This can be seen by comparing the blended field from the blended nowcast with the field from the non-blended nowcast (the two right-most panels in the bottom row of Fig. 3). Additionally also the noise has a reduced amplitude (due to its initialization to zero) which results in a blended field with a too small power (= St. Dev.) as seen in the right-most upper panels in Fig. 3.

For the largest scales, there is also a loss of power in the blended field. This is mainly due to the spin-up of the noise field. Due to the large temporal-autocorrelation they stay close to zero for a long time. Keeping in mind the blending strategy, I guess the St. Dev. of the noise field should be (close to) 1 and the blending weights should handle their contribution to the blended field. But since this is not the case the overall power of the blended fields is too low for the largest scales (as seen by comparing the two right-most upper row panels in Fig. 4). overview_l2_t10 Fig. 4: AR(2)-evolution overview for cascade level 2 at time step 10

These issues lead to two effects: The first one being the loss of small scale features in the first time steps of the nowcast (which triggered this investigation). \ The second one, which is more subtle, being the large scale noise having not enough power in the beginning of the nowcast. This syndrome can be seen when there is precipitation from the NWP-forecast outside the radar-domain (e.g. in Fig. 2). In this case the noisy fields (since the NWP blending weight is typically rather small) appear with rather small rainfall amounts (which actually gives them a quite realistic feel, but theoretically their amplitudes are too small.)

Unfortunately there is another problem. Namely with the semi-Lagrangian advection of the small scales. In the blending nowcast all the different cascade levels need to be advected separately because they need to be blended with the NWP cascade levels. But for the smallest levels this also leads to a loss in power (as seen in Figs. 5 and 6) noise_adv_l7 Fig. 5: Standard deviation (~ power) of the noise field before and after semi-Lagrangian advection for level 6. noise_adv_l6 Fig. 6: Standard deviation (~ power) of the noise field before and after semi-Lagrangian advection for level 7.

The reason for this is less clear to me, but I suspect this has to do with the interpolation needed to find the starting points of the advection. Due to the high spatial variance of the smallest levels, potentially a large error is introduced when interpolating.\ The non-blending nowcast does not suffer from this because their the levels are recomposed to the total precipitation field before the advection takes place.

Some possible solution:

  1. Renormalize the precipitation cascade after each AR(2) iteration
  2. Initialize the noise cascade with the appropriate noise

I did a first implementation if these ideas which resulted in the following: overview_l6_t1_adapted Fig. 7: AR(2)-evolution overview for cascade level 6 at time step 1 after adaptation

overview_l2_t10_adapted Fig. 8: AR(2)-evolution overview for cascade level 2 at time step 10 after adaptation

As Figs. 7 and 8 show, the blended cascade levels now have again the correct St. Dev. (or power) after the AR(2) evolution. The resulting nowcast is shown in Fig. 9

blendVSnoblendVSadapted_t1 Fig. 9: The blended (left) and non-blended (middle) and adapted-blended (right) nowcast after 1 time step of 5 minutes

As you can see the adapted nowcast retains better the small-scale features in the first time steps. Also notice the precipitation outside the radar domain immediately coming in at full power (a result of the noise component of the large scales having the correct power).

However, these adaptations still have some drawbacks:

@ladc and myself have reserved some time in July to further investigate this issue. Comments and suggestions are much appreciated.

Kind regards, Michiel

RubenImhoff commented 1 year ago

Hi @mpvginde,

Many thanks for this thorough analysis! We already discussed it earlier, but I'm very happy to see you went through the code and these results in so much detail. And good to see this as a pysteps issue now, so that we can think along in coming up with a solution.

Regarding initializing the noise cascade with the appropriate noise, I think this is the way to go and what you have tested, could in my opinion directly be implemented in the code.

Regarding the renormalization, I remember that @ladc and I discussed this and I can't entirely remember why we decided not to implement that (directly). I think we were worried about any issues at the end of the forecast when the x times renormalized cascades had to be recomposed. Have you tested this by any chance? If there are no issues, I think this would be a great addition as well. The localized version of the AR(2) process would even be better, perhaps this is eventually the way to go for the blending procedure. Maybe we should consider this as step 2 in the improvements? I can imagine that the localized version also takes somewhat more implementation time.

Finally, regarding the semi-Lagrangian advection related issues, perhaps @pulkkins has ideas here?

Cheers,

Ruben

dnerini commented 1 year ago

Very interesting issue and a lot of detail to digest! Many thanks for the analysis @mpvginde!

The noise starts from 0 and needs to spin-up. For the largest levels (due to the large temporal autocorrelations) this may take a long time.

I agree this doesn't sound right. Should we consider it as a bug or does this follow the original formulation in the paper? Can you submit a PR with the fix you've tested above?

The AR(2) evolution of the precipitation cascade without the noise term EPS is no longer stationary in the variance (something that is necessary according to Bowler et al. 2006

Again, should we consider this a bug? From you answer @RubenImhoff , it seemed that it was rather a design choice, was it?

Finally, I'd like to ping @albertocarpentieri here, since he encountered a similar issue in his work (preprint ), although it was related to the extrapolation only case without blending. But he also added a renormalization step after each AR iteration, so maybe there is a more general solution that we can discuss with respect to this?

dnerini commented 1 year ago

@mpvginde I sent you an invitation to join the pysteps organization. Feel free to create new branches directly in the pysteps repo.

RubenImhoff commented 5 months ago

Hi @mpvginde, is it also getting time to merge our blending improvements to the main branch of pysteps, as that will solve this issue?