RainFARM : Terzago et al. 2018

[simonbeylat]

I wanted to ask you if you could add the version of RainFARM improved by Terzago et al.(https://doi.org/10.5194/nhess-18-2825-2018).

I would like to use it using Python. By the way, I am open to collaborate and help you to do it.

[dnerini]

Hello @simonbeylat

I'm not familiar with the work of Terzago et al, can you please briefly summarize what they improved with respect to the original RainFARM?

In any case, we would of course welcome your contribution! You can submit a draft of this new RainFARM version as a pull request following our developer guidelines. We will then support you during the review of your code, by suggesting code changes and contributing with code to make sure that your contribution will be properly integrated in pysteps.

Before starting your implementation, it might be worth discussing some higher level details here. For example, should this become a separate method, an option in the existing RainFARM method, or simply replace some existing code? For this, I come back to my initial question: please try to describe in more detail what you'd like to do and how you propose to implement them.

[simonbeylat]

Hello @dnerini,
Thank you for your reply.

Terzago has added 2 elements to RainFARM. The first is a smoothing operator at the end of the algorithm.

The second (and most important) is a weight applied at the fine scale of the field, before the last step.
These weights are calculated from the climatology. It allows RainFARM to reproduce the climatology of fine-scale precipitation.
The aim is to provide RainFARM with complex orographic information in order to reproduce a realistic fine-scale precipitation field.

I think we can follow what von Hardenberg (one of the authors of the RainFARM articles) has done in R and Julia and add an option to the current function. (https://github.com/jhardenberg/rainfarmr).

[dnerini]

I like your proposition, @simonbeylat! If possible, I'd keep the two changes separate, hence work on two pull requests. The smoothing operator looks like a fairly straightforward addition, while the weighting is somewhat more involved, as it requires a climatological analysis. In this sense, I assume we would need to implement the methods to compute such weights based on a data archive?

Also, keep in mind that currently the rainfarm implementation in pysteps (courtesy of @jleinonen) doesn't include the temporal component, would you be willing to work on this missing part, too?

Finally, I was thinking we should get in touch with @jhardenberg, perhaps he's interested in contributing to a python version of his code within pysteps?

[simonbeylat]

Great, it's a good idea to see if von Hardenberg is interested.

So I think we can do a first pull request for the smoothing operator where we create a little internal function to do that, and add a boolean parameter called Smoothing to the RainFARM function (downscale) and add this smoothing option to the main function.

A second pull request with a function to compute the weight that has as parameter the climatology we want to use, the output of this function would be the weights.
Then add a parameter called Weight to the RainFARM function with None as the default value, in case we don't want to use the weight, otherwise you can just set this parameter to the output of the weight function.

Do you like it this way?

For the time component, it will be interesting to work on, but maybe we need to think a bit more about how we can implement it.

I assume that you used this environment to develop this package. Do you think I can use the xarray library, which I like, to read NetCDF files? Otherwise I would adapt

[jhardenberg]

Hi folks, I have actually seen this development and it would be really cool to have a python implementation of the more recent developments. The Terzago et al. 2018 version actually differs from the original Rebora et al. 2006 version also in another detail: we introduced (actually already in D'Onofrio et al 2014 https://doi.org/10.1175/JHM-D-13-096.1 ) a 'spectral merging' feature which greatly increases the realism of the resulting fields. Indeed this occured at the cost of dropping the downscaling in time. Actually, I will give a look at the current implementation in pysteps and the come back with some comments.

[jhardenberg]

The time downscaling was easy to implement in the original Rebora et al. 2006 by using 3D FFTs, but actually in D'Onofrio et al. 2014 we dropped it for climate downscaling when we introduced the spectral merging. Indeed the entire concept of downscaling in time made a lot sense at the weather event scale, but is less well defined for climate fields (should we do a fft in time over all the length of the simulation ? Only over pieces , which pieces? ). Also, as long as daily precipitation is considered, correlations in time are guaranteed by the conservation of the large-scale field by RainFARM (we do not expect significant subgrid correlations in time at the subdaily scales if we downscale for example from 50 to 1 km ....) so it is not needed.

[simonbeylat]

Thank you for your reply @jhardenberg

In this case, it might be better to start by implementing this "spectral fusion" of D'Onofrio et al 2014. To follow the same historical development.

[dnerini]

Very good @simonbeylat, I very much like your plan and look forward to review and contribute to your PRs!

I assume that you used this environment to develop this package.

For development you can use the dev environment. You should be able to follow these instructions, let me know if you encounter any problem.

Do you think I can use the xarray library, which I like, to read NetCDF files? Otherwise I would adapt

The introduction of xarray has been a long-standing issue, but as of today, we still haven't adapted pysteps (it's work in progress although somewhat stale unfortunately...) . So this means that we would prefer not to add xarray as a dependency (not for now), but this doesn't mean that you shouldn't use it for developing and testing your implementation locally. In the end, the methods should accept a numpy array, which you can always easily get from xarray, if you see what I mean.

[dnerini]

Also (and apologies for the slow answer), hi @jhardenberg and thanks for your inputs! Very nice that you're interested in contributing to this effort.

It's very interesting what you mention about the downscaling in time, it sure does make a lot of sense. I wonder whether we should still try to implement that component for weather applications, though. Something to keep in mind...

@simonbeylat how do you plan to proceed?

[simonbeylat]

Hello @dnerini,

I have some important things to finish (with deadline) for the next few weeks and then I'll have a few days off.
I think I'll focus fully on this in September, if that's okay for you ?

[dnerini]

Absolutely, there is no hurry, just wanted to have an idea ;-)

[simonbeylat]

Hi @dnerini



I am about to start the implementation, if I remember correctly we have 3 different implementations in 3 different pull requests (respecting what we said in July): 







1) Spectral fusion from D'Onofrio et al 2014. Directly integrated with the main function. 







2) The smoothing operator added with a boolean parameter called Smoothing to the RainFARM function (downscale) and a function that does it.







3) The weight function (which gives the weights using a climatology as input) and the weight parameter (where we can put the weights computed by the function).





Is it okey for you ?







[jhardenberg]

Hi @simonbeylat, thank you again for the initiative, that is great.
If I may suggest I would indeed start from the smoothing operator, which looks like a quite easy change. A reference implementation can be found both in my R and in my julia implementations. Notice that in both cases the smoothing is implemented as a convolution using the trick of transforming it into a product in spectral space. A tricky thing to pay attention to is the treatment of missing values, as you will see in both cases.

[dnerini]

Hey @simonbeylat very nice to hear from you ! The plan sounds very good, I look forward to see and discuss your pull requests!
Also, the suggestion from @jhardenberg sounds very reasonable, what do you think? Shall we start with the smoothing operator?
Final note, do not hesitate to draft a pull request early on with a preliminary version of your code: this often helps to get feedback and get things right from the beginning (and avoid dead ends).

[simonbeylat]

Hi @dnerini

I have created a first pull request with the smoothing operator.

I used @jhardenberg's R script to do it! I tested it with the same input and got the same result (with R and Python).

I'm waiting for your comments and starting the next pull request which should be the "spectral fusion of D'Onofrio et al 2014. Directly integrated with the main function. "

[simonbeylat]

Hi @dnerini @jhardenberg,

I have looked at the articles by D'Onofrio et al.(2014), if I understand correctly the change is related to section 3 i) and iv). 
I was wondering if we should incorporate it into the main function or add it as an option (with a boolean argument).



I think the boolean argument option might be interesting, especially if we want to focus on time downscaling.



As I understand it, D'onofrio et al. have removed the time part. So I would say that we need to be able to switch off this change in case we want to add the time downscaling latter. 






What is your opinion ?

[dnerini]

Hi @simonbeylat, I agree, the spectral fusion can be added as an option to the main downscaling function, and, as you mention it, such option will be set to false when downscaling in time.