Allow for semi dynamic incremental mask in blending setup

[RubenImhoff]

This PR fixes issue #381.

In the extrapolation based nowcasting, the mask is dilated iteratively at every timestep. This is easy to do in Lagrangian space for a pure extrapolation-based nowcast. However, in the case of STEPS blending with NWP, this is less straightforward because we are no longer in Lagrangian space (NWP is not advected). Right now, it's dilated only once but this does not grow further per timestep.

I've made a first start to make the incremental mask computation more dynamic (i.e., more incremental) over time. I've done so by including a max_mask_rim next to the mask_rim in the mask_kwargs. With increasing lead time, the mask_rim is then allowed to increase with the time step index on top of the mask_rim up to the point where max_mask_rim is reached. This to ensure that we do not increase run times too much, while a very large mask_rim is not doing much as it gets cut off by the rainfall threshold at some point.

Another option would be to make the iterate_structure setup below dependent on the time step:

# iterate it to expand it nxn
n = mask_f * self.__config.timestep / self.__config.kmperpixel
self.__params.struct = iterate_structure(struct, int((n - 1) / 2.0))

However, this should be stopped after a certain amount of lead times, as the struct parameter becomes way to big (and slow) otherwise.

Or, if any of you has other suggestions, please let me know!

To do's:

• Implement first approach to improve the current incremental mask setup.
• Decide on final approach.
• Add tests.

[RubenImhoff]

To show what typically happens, after some time steps (often the longer lead times during the blending), the rainfall fields that predominantly follow from the noise generation tend to get cut off too strongly by the mask:

With the suggested changes, this will change to (which I think looks more natural):

By the way, the 'hard' edges in the north and east are the edge of the domain (this is zoomed in).

[sidekock]

To show what typically happens, after some time steps (often the longer lead times during the blending), the rainfall fields that predominantly follow from the noise generation tend to get cut off too strongly by the mask:

With the suggested changes, this will change to (which I think looks more natural):

By the way, the 'hard' edges in the north and east are the edge of the domain (this is zoomed in).

Hi Ruben,
I think I understand what the above comments are talking about but just to be sure, would it be possible to draw exactly which parts change? For instance, would it be possible to plot the mask with some sort of transparency?
Thanks!

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 84.55%. Comparing base (5fd4dc8) to head (d44cb1f).

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #459   +/-   ##
=======================================
  Coverage   84.55%   84.55%           
=======================================
  Files         162      162           
  Lines       13458    13460    +2     
=======================================
+ Hits        11379    11381    +2     
  Misses       2079     2079           

Flag	Coverage Δ
unit_tests	84.55% <100.00%> (+<0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[RubenImhoff]

Yes, good suggestion, @sidekock. I went for the quick and dirty way (otherwise it becomes quite some work to plot this nicely including the final precip fields):

The figures below show 1 = within mask and 0 = outside mask (so masked out in the post-processing). I've chosen a max_mask_rim of 40 (which maybe could even be lower for a similar effect, saves runtime) and mask_rim the default of 10.

At first lead time (both masks should still be the same), left is proposed approach, right is 'old' approach:

Lead time 10 (15-min time step, so after ~2.5 hours):

Lead time 20 (so after ~5 hours):

Lead time 52 (so after 13 hours):

The rainfall fields (forecast) for these lead times were as follows:
Lead time 1:

Lead time 10:

Lead time 20:

Lead time 52:

I still see some strong edge effect, @sidekock, while I'm using a smooth_radar_mask_range of 100. Any ideas on that?

[RubenImhoff]

Effect of a smaller max_mask_rim = 20 at the last lead time:

Slightly different from max_mask_rim = 40.

The runtimes were hard to measure in a clean way on my laptop, due to side processes, but runtimes are about 6% higher for max_mask_rim = 40 and about 2% higher for max_mask_rim = 20.

[mats-knmi]

I don't know a lot about this part of the code, but the proposed solution looks good to me. A simple way to achieve the goal. I agree that this results in more realistic looking precipitation fields.

Is there a downside to this approach vs the other option or any other possible solution?

[sidekock]

@RubenImhoff, thanks for these plots. Could I ask for some more :)

• I was wondering if I could see the difference in output of the final radar fields between 10 and 40 max_mask_rim. That way I can better understand the effects of the rim... And maybe give you an answer on the edge effects...

Currently, I have no clue what is going on here. Maybe just a check but what is the data just pat that edge on timestep 52? are those nans? I assume so and that is indeed what my implementation is using.

[RubenImhoff]

@RubenImhoff, thanks for these plots. Could I ask for some more :)

• I was wondering if I could see the difference in output of the final radar fields between 10 and 40 max_mask_rim. That way I can better understand the effects of the rim... And maybe give you an answer on the edge effects...

Currently, I have no clue what is going on here. Maybe just a check but what is the data just pat that edge on timestep 52? are those nans? I assume so and that is indeed what my implementation is using.

You mean the forecast rainfall fields, then the one of 40 mask_max_rim is already there (lead time 52 in previous comments), but I can add 10 (if that is what you were looking for).

The edge effects actually seem strongest at timestep 1, which is indeed where there are nans just at the edges of the domain.

[sidekock]

@RubenImhoff, thanks for these plots. Could I ask for some more :)

• I was wondering if I could see the difference in output of the final radar fields between 10 and 40 max_mask_rim. That way I can better understand the effects of the rim... And maybe give you an answer on the edge effects...

Currently, I have no clue what is going on here. Maybe just a check but what is the data just pat that edge on timestep 52? are those nans? I assume so and that is indeed what my implementation is using.

You mean the forecast rainfall fields, then the one of 40 mask_max_rim is already there (lead time 52 in previous comments), but I can add 10 (if that is what you were looking for).

The edge effects actually seem strongest at timestep 1, which is indeed where there are nans just at the edges of the domain.

Indeed, I want to compare the 10 to the 40 on the actual rain values and not just the mask. I'll dive into the code to check what is going on with this mask. I don't fully know what is going on, but I will run some tests and plot the results to see what is going on...

[RubenImhoff]

I don't know a lot about this part of the code, but the proposed solution looks good to me. A simple way to achieve the goal. I agree that this results in more realistic looking precipitation fields.

Is there a downside to this approach vs the other option or any other possible solution?

Not really a downside, I think. It increases computation time slightly and the buffer cannot grow indefinitely (but maybe that is okay?), because the max_rim also results in diluted values towards the edges of the buffer. That means that at some point, the values become so low that they are cut off by the used rainfall threshold. Although I find it hard to judge if this is the right approach, the advantage of this approach is that it uses already existing functionality in pysteps.

[RubenImhoff]

This is for max_mask_rim = 10 at lead time 52; 40 is already plotted above.

As you can see, the differences with 40 are small, but become visible for the smaller rainfall fields/cells that pop up. So especially for convective days, this will make a difference (and that is where I also used to see the cut off that was too strict).

[ladc]

Thanks a lot for your work on this, I was also looking for the source of these undesirable "diamonds" and a possible solution. The effect was especially noticeable on the probability plots. Did you have a look at these in the new version?

[sidekock]

@RubenImhoff would it be possible to give me the exact configuration of you blended pysteps run? That makes it easyer to run some tests

[ladc]

Thanks a lot for your work on this, I was also looking for the source of these undesirable "diamonds" and a possible solution. The effect was especially noticeable on the probability plots. Did you have a look at these in the new version?

PS: actually the maps you gave sent are clear enough and probably the same pattern can be expected in the probability maps

[sidekock]

@RubenImhoff I took another look at the plots you sent over and I think the issue is just the data you are using. Correct me if I am wrong but as far as I can see, this is what I observe:

• This dataset is not in the same projection as the map
• The edge of the domain is visible on this plot.
• Am I correct that the following is roughly the study domain:
  
  If my previous assumptions where right, the visual blend of radar mask will not have an effect. I implemented a blend for a smooth transition between the radar mask (which usually covers a smaller domain than the NWP). In this case since both domains outer coverage would be the same, there is nothing to be done. If some of the assumptions are wrong. Please let me know

Hi @sidekock, you are right. That might be it then. If I see any other strange effect for later timesteps, I'll let you know. But this explains it for now, I think.

[RubenImhoff]

Shall we continue with this method? Then I'll add some tests. :)

[dnerini]

Shall we continue with this method? Then I'll add some tests. :)

why the question? is there any argument against it? :)

[RubenImhoff]

Shall we continue with this method? Then I'll add some tests. :)

why the question? is there any argument against it? :)

None so far - just checking whether we agree to indeed go for it. :)

[ladc]

Looks good, Ruben! I'll just push a tiny simplification and run a test and get back to you.

[ladc]

Sorry, it seems I accidentally marked this ready for review. Anyway, I checked the output with a max_mask_rim of 40 and it does look more physical with less artifacts. These changes are all good for me. Though I suppose an extra test wouldn't hurt.

[ladc]

Lookin' good.

[dnerini]

@RubenImhoff Is this ready to
Be merged or were you planning additional work ?