Spatial scale values in pySTEPS

[tkokkone]

Dear pySTEPS experts,

This may be a pySTEPS issue, but just as well it can be an issue with my maths :-) Either way, I would really appreciate if someone could help.

I am having a hard time in trying to understand the relationship between spatial scales and cascade level indices the way it has been implemented in pySTEPS. I try to explain with an example where my problem is.

According to Figure 1 in Pulkkinen et al. (2019) (https://gmd.copernicus.org/articles/12/4185/2019/) the spatial scales corresponding to cascade levels are as follows:

The last Ratio column shows the ratio between two successive spatial scales. This ratio is constant for all other pairs of scales but for the first two. I also included the central wavenumbers in the table above (_centralwn).

The relationship between the central wavenumbers and spatial scales is (except for level 0):

Dim / _centralwn / 2

where Dim is the maximum grid dimension. This makes some sense to me as the wavenumber is defined as the reciprocal of the wavelength and the spatial scale as half of the wavelength. But what should happen at level zero where the wavenuber is zero, and hence the wavelength and the spatial scale are infinite? Or not? Anyway, it is the zero cascade level that causes the ratio between subsequent spatial scales not to be constant, which poses me a problem as I will explain in the next paragraph.

My colleagues and myself at Aalto University have mainly used STEPS/pySTEPS in ‘simulation mode’ to create an ensemble of rain fields (see e.g. Niemi et al., 2016, DOI:10.1002/2015WR017521). Now the autocorrelative properties between consecutive rain fields cannot be estimated in the same manner as in nowcasting mode, but a relationship between the life-time and the spatial scale of rainfall structures needs to be established. Here we have followed equations 9-11 in Seed at al. (2014)

https://www.researchgate.net/publication/264943060_Stochastic_simulation_of_space-time_rainfall_patterns_for_the_Brisbane_river_catchment

From the relationship between the correlation length and the spatial scale (equation 9) one can observe that the logarithm of the spatial scale and the logarithm of the correlation length are linearly related. Parameters a and b (equation 9) can hence be estimated as the intercept and slope of the line that forms when ln(correlation length) is plotted against ln(spatial scale). I attach below such figure for one event.

The grey line represents the data where spatial scale values have been calculated as in Figure 1 in Pulkkinen et al. (2019) (see above). One can clearly see the “jump” as a longer horizontal distance between the two last points. Also, it seems that as the consequence of the “jump” the last point does not fall on the line. The orange line is otherwise exactly the same, except for the spatial scale value of the last point, which has been modified so that the ratio between spatial scale values for neighbouring cascades remains constant. In practice, in the orange line the second last spatial scale value of 44 (ln(44) = 3.78) has been multiplied with the constant scale ratio of 2.1 to yield for the largest scale a spatial scale value of 93.8 (ln(93.8) = 4.5) – instead of 264 (ln(264) = 5.6) given by the pySTEPS convention. This scale ratio (here 2.1), as discussed above, is in pySTEPS constant for all other cascade pairs but for the last two. Now all points lie neatly on the line.

In Alan Seed’s and Australian Bureau of Meteorology’s C++ STEPS version the spatial scales are computed as follows (abbreviated):

r4ScaleRatio = pow( 2.0/(double)i4CascadeSize, 1.0/(double)(i4NumberLevels-1)) ;

float r4Scale = i4CascadeSizer4PixelSize; for ( i4Level = 0; i4Level < i4NumberLevels; i4Level++ ) { r4Scale = r4ScaleRatio; }

So evidently the scale ratio is constant for all neighbouring cascade pairs. This can also be seen from Equation 1 in Seed at al. (2013) (https://doi.org/10.1002/wrcr.20536).

So, I guess my questions are:

Do the spatial scale values in pySTEPS, as it seems, have a “jump” in the scale ratio between the two largest cascade scales? If yes, should it be so? If still yes, how and where has this been implemented (presumably in bandpass_filters.py in filter_gaussian method) in pySTEPS? As pySTEPS does not explicitly compute/report spatial scale values (as far as I know…) we have reverse engineered them from the central wavenumbers (returned by filter_gaussian) and Figure 1 in Pulkkinen et al. (2019) (see above).

Any help is greatly appreciated!

[dnerini]

Hi @tkokkone thanks for the interesting and detailed questions! @pulkkins is our expert on the subject, so I'll leave it to him to provide a response. And as you probably already realised, sometimes this can take some time, as we mostly take care of pysteps on a "best effort" basis... but please keep asking good questions and sharing results from your interesting work!

[tkokkone]

Hi @dnerini, @pulkkins and others,

Thanks for the reply! I am not expecting immediate responses, no worries, but great to see that this repository is alive! I guess you guys are my last resort to resolve this but there is no panic.

[pulkkins]

@tkokkone you are absolutely right. It would make more sense to choose the central frequencies so that the ratio between the two largest scales is consistent with the smaller ones. This ratio is in fact governed by the parameter l_0 of cascade.bandpass_filters.filter_gaussian: https://github.com/pySTEPS/pysteps/blob/9c252b5d87653cfe3b4ea18b41ff0d78193a90ce/pysteps/cascade/bandpass_filters.py#L102 In the above case with domain size 1226 and 8 cascade levels, such a value for l_0 is approximately 1.380966. I think we should implement a formula for automatically choosing l_0 so that the ratio between all successive scales is constant and use that as the default choice. I'll take a look into that.

[tkokkone]

Thanks, @pulkkins ! I am happy to help if I can.