Linear blending

[wdewettin]

This addresses issue #226. I added a nowcasting method called linear blending. This works for all nowcast methods (deterministic and ensemble) and is also compatible with NWP ensembles. Furthermore, I added an example which uses some dummy NWP data to blend with and which plots the results. Lastly, I also added a test for the linear blending function.

[RubenImhoff]

Great work, @wdewettin! I was (actually, we were - that is including @wdewettin and @ladc) wondering whether we should place this method in a/the blending folder or whether it fits in the nowcasts folder as currently done by @wdewettin. Any thoughts?

[RubenImhoff]

Great, @wdewettin! I left some comments and questions. Some of them I could implement myself as well. Let's just wait for the comments of the others (and feel free to comment yourself on them), then we'll continue with the final changes.
Exciting to have this option in pysteps soon!

[RubenImhoff]

In case of R_nwp = None, the script is similar to steps.py, right? Do we want to use the script as a postprocessing tool, i.e. the user makes a nowcast with steps.py and passes the output on to linear_blending.py, or do we want to keep it this way and have some duplicate code with steps.py? As steps.py may also change to allow for other blending methods, I wonder what everyone thinks about the best way forward.

[cvelascof]

Here the object returned should be equal to the output of the nowcasting method used if R_nwp=None.
I understand that if we are going to create a function the ingest both sources of precipitation, then few lines of code must be duplicated. Perhaps in the future a higher level blending script should be created in the way that accept multiple sources of precipitation, and weights for each source (or a method to determine those weigths) and return blended fields. This approach may be easier to maintain, but perhaps harder to understand by users.

[RubenImhoff]

I believe the threshold is present in the metadata (the transformer at least returns the metadata). We could add this to the function, then we don't have to hard code it here.

[cvelascof]

I understand why this is done here, but this is assuming that R_nwp is in mm/hr and that the precip is in dB. If this is the case, then we need to indicate those limitations in the descrition of the function. Alternatively, both sources should be in mm/hr and the conversion of precip to dB be done within this function as well.

[cvelascof]

I understand why this is done here, but this is assuming that R_nwp is in mm/hr and that the precip is in dB. If this is the case, then we need to indicate those limitations in the description of the function. Alternatively, both sources could be in mm/hr and the conversion of precip to dB be done within this function as well.

Said that the returned output R_blended will be always in mm/hr

[RubenImhoff]

After consulting @wdewettin and @ladc, I figured out what this piece of code does. So if you e.g. have 10 ensemble nowcasts members and 3 NWP members, the output will be an ensemble of 10 members. Hence, the three NWP members are blended with the first three members of the nowcast (member one with member one, two with two, etc.), subsequently, the same NWP members are blended with the next three members (NWP member one with member 4, NWP member 2 with member 5, etc.), until 10 is reached.

Is this the way to go (works for me) or do we want an output of 30 (10 x 3 members here) in this case?

[RubenImhoff]

Oh and @wdewettin, let's add some comments to the code here so that this procedure is clear. :) (B.t.w., I can add that too, but I'll wait for the responses - and thus possible changes - of the others)

[wdewettin]

Yes! I will await the responses from the others, and then add some comments

[cvelascof]

I agree. It wuld be nice to add additional comments in this part of the script as @RubenImhoff suggested.

[cvelascof]

I agree. It would be nice to add additional comments in this part of the script as @RubenImhoff suggested.

[RubenImhoff]

With the way it is implemented at this moment: n_ens_members equals the maximum number of ensemble members in the input, so either the number of the nowcast or that of the NWP. Perhaps good to add that here.

[RubenImhoff]

(B.t.w., I can add that too, but I'll wait for the responses - and thus possible changes - of the others)

[RubenImhoff]

If not used, we could remove these lines.

[wdewettin]

True, in the plot_steps_nowcast.py example, this is how the missing data was handled, but here we use the (dummy) NWP data

[RubenImhoff]

Perhaps plot the radar rainfall and NWP rainfall forecast field somewhere, to give the reader (as this example may end up in the pysteps gallery) and idea of what the example fields look like.

[wdewettin]

Do you mean plotting them for every timestep? Or only the first one(s)?

[RubenImhoff]

First ones should be fine. :)

[dnerini]

Great work, @wdewettin! I was (actually, we were - that is including @wdewettin and @ladc) wondering whether we should place this method in a/the blending folder or whether it fits in the nowcasts folder as currently done by @wdewettin. Any thoughts?

Because of the way linear blending works, I personally think it could be provided as a simple postprocessing method (therefore within the postprocessing module) that takes as inputs two separate nowcasts and the relatives weights.

This said, I like the way @wdewettin implemented it, where one has to specify the the nowcasting method and pass the corresponding NWP forecast. But because of this dependency on the actual nowcasting methods, it could make more sense to move into a new blending module. We could say that as a general (and not necessarily very strict) rule, the methods included in a module should be independent one to the other. In this sense, the nowcating module should only include alternative nowcasting methods.

[dnerini]

Very nice contribution @wdewettin ! I highly appreciate the fact that you already included all the documentation (including a gallery example!) and tests, well done!

I plan to go through your PR in detail as soon as possible. For now, I have a more general request. Would it be possible to make your method work with an xarray DataArray object instead of numpy arrays?

[wdewettin]

Yes, @dnerini, I will take a look at making my method work with xarrays!

[wdewettin]

@dnerini Sorry for my late response! I looked into making it work with xarrays instead of numpy arrays and I found that it would be better to wait with that for the time being. All the functions on which the linear blending function is based still output numpy arrays, and it will be a lot easier to make the linear blending function work with xarrays if those functions are changed first.

[codecov]

Codecov Report

❗ No coverage uploaded for pull request base (pysteps-v2@3802f8d). Click here to learn what that means.
The diff coverage is n/a.

@@              Coverage Diff              @@
##             pysteps-v2     #229   +/-   ##
=============================================
  Coverage              ?   73.83%           
=============================================
  Files                 ?      152           
  Lines                 ?    11213           
  Branches              ?        0           
=============================================
  Hits                  ?     8279           
  Misses                ?     2934           
  Partials              ?        0           

Flag	Coverage Δ
unit_tests	73.83% <0.00%> (?)

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

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3802f8d...758f88b. Read the comment docs.

[cvelascof]

@RubenImhoff @wdewettin thanks for implementing this module.

[cvelascof]

Here the object returned should be equal to the output of the nowcasting method used if R_nwp=None.
I understand that if we are going to create a function the ingest both sources of precipitation, then few lines of code must be duplicated. Perhaps in the future a higher level blending script should be created in the way that accept multiple sources of precipitation, and weights for each source (or a method to determine those weigths) and return blended fields. This approach may be easier to maintain, but perhaps harder to understand by users.

[cvelascof]

I understand why this is done here, but this is assuming that R_nwp is in mm/hr and that the precip is in dB. If this is the case, then we need to indicate those limitations in the descrition of the function. Alternatively, both sources should be in mm/hr and the conversion of precip to dB be done within this function as well.

[cvelascof]

Great job @wdewettin and @RubenImhoff

[cvelascof]

I agree. It wuld be nice to add additional comments in this part of the script as @RubenImhoff suggested.

[cvelascof]

I agree. It would be nice to add additional comments in this part of the script as @RubenImhoff suggested.

[cvelascof]

I understand why this is done here, but this is assuming that R_nwp is in mm/hr and that the precip is in dB. If this is the case, then we need to indicate those limitations in the description of the function. Alternatively, both sources could be in mm/hr and the conversion of precip to dB be done within this function as well.

Said that the returned output R_blended will be always in mm/hr

[cvelascof]

@RubenImhoff @wdewettin @dnerini apologies for the mess in my comments ... not sure what did it go wrong?

[RubenImhoff]

All requested changes implemented. Main task left is to make the code xarray ready for pysteps-v2. As discussed before, we may also do that in a separate PR.