Main module blending

[RubenImhoff]

First 'steps' towards a main module for a blended nowcasts-NWP forecast following the original STEPS approach as described in Bowler et al. (2006) and Seed et al. (2013).

• Working main module
• Tests
• Examples
• Documentation
• Probability matching following a different method
• Implementation of blending weights as in Bowler et al. (2006)
• Implementation of blending weights as in Seed et al. (2013)
• Fixes to prability matching method, masking method and noisy output
• blending.utils tests

Part of issue #212

[RubenImhoff]

pystepsrc

You have included a new path_workdir. I'm not against, but what's its purpose? Perhaps you could include a short comment to help understand its scope?

Good point. We implemented this as a place to store the decomposed NWP fields and NWP velocity fields, but seeing how it used now, we can reach the same with 'just' the path_outputs in pystepsrc. @ladc, what do you think?

[RubenImhoff]

In the blended forecast, I have the impression that something is a bit off with the NWP values, as I observed a strong increase of the intensities with leadtime. I anyway still need to have a closer look at your implementation of the masking and probability matching.

I completely forgot to have a new look at the examples forecast after the latest changes. The minor code changes will follow soon. Indeed, something seems to go wrong with either the masking or prob_matching.. I don't see the same effect (at least it is not as clear if it is there) in the test forecasts with the RMI data, so I wonder what goes wrong.

[RubenImhoff]

In the blended forecast, I have the impression that something is a bit off with the NWP values, as I observed a strong increase of the intensities with leadtime. I anyway still need to have a closer look at your implementation of the masking and probability matching.

I completely forgot to have a new look at the examples forecast after the latest changes. The minor code changes will follow soon. Indeed, something seems to go wrong with either the masking or prob_matching.. I don't see the same effect (at least it is not as clear if it is there) in the test forecasts with the RMI data, so I wonder what goes wrong.

While working on the tests for pysteps.blending.utils, I came across an error in pysteps.blending.utils. In the function stack_cascades, the mean and sigma could be set to 0.0 and 1.0 (I don't know why), which eventually gives wrong values after blending with NWP. As this only happened with the NWP component(s), the prob_matching was compensating for this and that is why we saw the effect with increasing time in the examples scripts. However, the changes above fix the problems (nevertheless, masking and prob_matching still require a critical look)! Also for the RMI test scripts the results look better:

[dnerini]

Very nice @RubenImhoff! The animation is very promising!

I'm having a closer look at the probability matching. Do I understand it right that you basically perform a linear blending of radar extrapolation and NWP based on the weights from the second cascade level in order to build the reference field for the probability matching? This is now a bit of a detail, but the problem with this is that you would loose the peak values, especially when both weights are around 0.5.

An alternative approach would be a simple resampling: for each grid point you sample either the radar or the NWP value (and this could be based on the same weights you are already using). This way you would build a new target distribution that is not smoother than the original radar and NWP fields.

I have some code for this (I used it in the EnKF approach). Should I give it a try to implement it in your code?

[RubenImhoff]

An alternative approach would be a simple resampling: for each grid point you sample either the radar or the NWP value (and this could be based on the same weights you are already using). This way you would build a new target distribution that is not smoother than the original radar and NWP fields.

I have some code for this (I used it in the EnKF approach). Should I give it a try to implement it in your code?

Great idea, @dnerini, and thanks if you want to implement that!

[RubenImhoff]

Yesterday, @ladc, @cvelascof and me had a nice talk with Alan Seed about the blending code and the weights implementation. He was happy to see that we're making this open source in pysteps and thought that the results so far look promising. Based on our talk, some main suggestions from Alan (and things to work on probably in the coming days):

• Have a look at default skill values for cases where the weights implementation either becomes unstable or for when there is hardly any rain in the domain. This should ensure that the noise weights always increases with time and the nowcast weight decreases.
• Have a look at the combined weight term of the nowcast and NWP component (prior to the noise weight estimation) and make sure this always decreases with lead time.
• For the transformations, perhaps apply a log sine transformation which can include zeroes, thereby making the variance of the individual components more representative than with a threshold set (Alan said it is better not to mask the radar and NWP data with a given threshold prior to the blending). This is something we could implement. I wonder what the optical flow methods (and perhaps also the FFT) need, maybe we should at least keep a threshold for the motion estimations, but not for the remaining part of the blending?

[ladc]

pystepsrc

You have included a new path_workdir. I'm not against, but what's its purpose? Perhaps you could include a short comment to help understand its scope?

Good point. We implemented this as a place to store the decomposed NWP fields and NWP velocity fields, but seeing how it used now, we can reach the same with 'just' the path_outputs in pystepsrc. @ladc, what do you think?

Hi guys, sorry for the super late reply. I think we made the distinction because the one path is aimed at files for internal use, including input and output (also e.g. the weights file, pre-decomposed NWP, ...) and the other is where we store only the output products for users. But I also understand if you would prefer to reduce the number of different directories or parameters in the config file, and also store these in path_outputs.

[ladc]

That is not a bad idea, Daniele. How much time do you think we still need for v2? I'm thinking of our idea (at least @ladc and mine - but not limited to) to work towards a blending paper, including the open-source implementation of the STEPS code in pySTEPS and a hydrometeorological analysis. As long as we have a release ready by that time, it should suffice.

@ladc, what do you think?

The xarray overhaul (among other things) is indeed not a small task, and it takes some manpower to work on pysteps v2. If you think it's doable to integrate the blending into pysteps v1.6, then it might be better to go for that option indeed. I hope this does not create too much duplicate work?

Only changes I can think of, is that the importers and subsequent testing schemes are (partly) based on the xarray implementation.

[dnerini]

The xarray overhaul (among other things) is indeed not a small task, and it takes some manpower to work on pysteps v2. If you think it's doable to integrate the blending into pysteps v1.6, then it might be better to go for that option indeed. I hope this does not create too much duplicate work?

It should be doable in my opinion. To simplify the task, I suggest to exclude the nwp importers from the branch and publish them as separate pysteps plugins: https://pysteps.readthedocs.io/en/latest/developer_guide/importer_plugins.html

The rest of the code is already very much in line with v1 I believe? I'd also suggest a minimal refactoring to reorganize the main method in few components (I'd be happy to give it a try).

If we all agree, I'd start refactoring this PR so that it can be merged into master rather than pysteps-v2.

@RubenImhoff very cool that you could exchange with Alan about your work! One comment from my side on your list of improvements: at this point I'd try to limit adding complexity to this implementation. We want a robust base blending method that works, possibly as close as possible to the reference publication. Further improvements can come incrementally later on (with your publication? ;-)), and as separate modules, to avoid adding too much complexity to the base method. For example: climatological skill. This is an important component for an operational system (it has to run in any situation), but I would not include it in the base method. Instead, we can make sure that the method can accept is as an optional argument, and leave to the user to compute it in his/her implementation. Same idea with the log sine transformation: would be nice to experiment with it, but this can be included in a separate PR.

[RubenImhoff]

If we all agree, I'd start refactoring this PR so that it can be merged into master rather than pysteps-v2.

@RubenImhoff very cool that you could exchange with Alan about your work! One comment from my side on your list of improvements: at this point I'd try to limit adding complexity to this implementation. We want a robust base blending method that works, possibly as close as possible to the reference publication. Further improvements can come incrementally later on (with your publication? ;-)), and as separate modules, to avoid adding too much complexity to the base method. For example: climatological skill. This is an important component for an operational system (it has to run in any situation), but I would not include it in the base method. Instead, we can make sure that the method can accept is as an optional argument, and leave to the user to compute it in his/her implementation. Same idea with the log sine transformation: would be nice to experiment with it, but this can be included in a separate PR.

Sounds good, Daniele, and I agree! I'll upload a couple of minor changes and implement your idea to make the decomposed NWP input optional (so, you could also start with a non-decomposed NWP forecast) and then we're ready to go.

[RubenImhoff]

I've just merged Daniele's great refactoring work. Should we continue the review here, or should I merge this branch in the steps_blending branch and renew that for pysteps 1.6 (and review from there)?

[dnerini]

I've just merged Daniele's great refactoring work. Should we continue the review here, or should I merge this branch in the steps_blending branch and renew that for pysteps 1.6 (and review from there)?

Very good. Yes, please merge into steps_blending and we'll continue from there for the integration into master

[dnerini]

I've just merged Daniele's great refactoring work. Should we continue the review here, or should I merge this branch in the steps_blending branch and renew that for pysteps 1.6 (and review from there)?

Very good. Yes, please merge into steps_blending and we'll continue from there for the integration into master

mmmh maybe that's not a good idea if we want to merge into master after that. Let's instead merge into a new branch from master, I'll do it now.