Local Lagrangian probability forecasts

[dnerini]

I was wondering whether it could be possible to extend the extrapolation.semilagrangian.extrapolate method to include a stochastic component. This could be implemented as random deviation around the coordinates at which the precipitation field is evaluated. In practice, line 213

https://github.com/pySTEPS/pysteps/blob/e2c3484ef234e3ef1ac8f1f11b1e64bc6911f35c/pysteps/extrapolation/semilagrangian.py#L213

could become:

coords_warped = xy_coords + displacement + random_offset

where random_offset is a random displacement within an area which is proportional to the lead time (typically 1 km/min):

def sample_unit_circle(rows, cols):
    # source: https://stackoverflow.com/questions/46996866/sampling-uniformly-within-the-unit-circle
    length = np.sqrt(np.random.uniform(0, 1, size=(rows, cols)))
    angle = np.pi * np.random.uniform(0, 2, size=(rows, cols))
    dx = length * np.cos(angle)
    dy = length * np.sin(angle)
    return np.stack((dx, dy))

random_offset = sample_unit_circle(*xy_coords.shape) * td * 5

As a result, each extrapolation would sample the area around the original coordinate as in the Local Lagrangian approach of Germann and Zawadzki (2004). Worth noting that this approach would provide an efficient way to compute point probability forecasts (no spatial nor temporal correlation is considered).

There are of course a number of open issues, for example:

• How to efficiently account for the borders of the image?
• What's the best way to sample the offsets? Using a 2D Gaussian instead of uniform sampling could be a better idea.

I'll be happy to give it a try myself, but first it'd be good to hear your comments on this. Do you think this would be useful? Would the proposed implementation above work? Is is a good idea to try to implement it in extrapolation.semilagrangian or should we create a separate routine instead?

In particular, I suspect that @pulkkins and @loforest might have something to say in this regard ;-)

Reference

• Germann, U., & Zawadzki, I. (2004). Scale Dependence of the Predictability of Precipitation from Continental Radar Images. Part II: Probability Forecasts, Journal of Applied Meteorology, 43(1), 74-89. Retrieved Dec 19, 2020, from https://journals.ametsoc.org/view/journals/apme/43/1/1520-0450_2004_043_0074_sdotpo_2.0.co_2.xml

[pulkkins]

That would be a very nice addition to pysteps. In my opinion, we should to implement it in a separate routine and keep the basic method in extrapolation.semilagrangian.

[loforest]

The local Lagrangian approach would be a very good addition to pysteps, which is cheap to implement and to compute. It gives a simple probability nowcast when limited by computational resources (large domains, etc) or for less quantitative applications (smart-phone app, point probabilities are enough). From a research perspective, it may be worth performing a probabilistic verification to show the added value of generating stochastic ensembles.

The borders of the radar composite domain (or image) are an issue, but at the same time an intrinsic limitation of extrapolation nowcasts. Machine learning methods do not suffer from this issue and can continuously "create" rain close to the borders (although it is not clear to me the predictive value of this).

Whether to implement it in the existing function or not depends on how much code needs to be copy-pasted. If the change is small it may be worth keeping it in the same function. An additional argument would then set the rate of area growth per min. If set to None (default), we fall back to the deterministic extrapolation. Other questions: Is it better to pick random samples or space them regularly? Better circles or ellipses? See also: https://www.researchgate.net/publication/248403322_Precipitation_Nowcasting_Using_Radar-Derived_Atmospheric_Motion_Vectors

[dnerini]

Thanks @pulkkins and @loforest for the inputs. I drafted a PR #207 that provides a first implementation of what discussed above. It also includes an example. Please let me know what you think.