fix noise tapering and error handling

[mats-knmi]

Fixes #451

This applies the 2 solutions proposed in the issue. For the error handling I choose a more pragmatic approach than the one I proposed in the issue.

I define the no_precip booleans based on the zero mask of the tapering method provided. This should be exactly what we want, because these booleans are intended to represent whether we can use the specified field for the noise generation right?

@RubenImhoff: can you verify that this will not have unintented side effects?

[codecov]

Codecov Report

Attention: Patch coverage is 90.00000% with 4 lines in your changes missing coverage. Please review.

Project coverage is 84.24%. Comparing base (006fb49) to head (f71b2db).

Files with missing lines	Patch %	Lines
pysteps/blending/utils.py	50.00%	2 Missing ⚠️
pysteps/utils/tapering.py	66.66%	2 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #460   +/-   ##
=======================================
  Coverage   84.24%   84.24%           
=======================================
  Files         160      161    +1     
  Lines       13243    13256   +13     
=======================================
+ Hits        11156    11168   +12     
- Misses       2087     2088    +1     

Flag	Coverage Δ
unit_tests	84.24% <90.00%> (+<0.01%)	⬆️

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

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[RubenImhoff]

Hi @mats-knmi, fantastic, thanks a lot for picking this up! This seems to work, I think. Could you maybe also plot what the tapering looks like now? @bwalraven, would you like to have a look at it as well?

[mats-knmi]

Hi @mats-knmi, fantastic, thanks a lot for picking this up! This seems to work, I think. Could you maybe also plot what the tapering looks like now? @bwalraven, would you like to have a look at it as well?

This is what the tukey tapering looks like now for our domain (765x700):

And this is what it would look like for @bwalraven 's domain (199x113):

I tried to stretch the images so that x and y are the same scale.

[RubenImhoff]

Good and when win_fun is set to None, all are set to ones? (Is that recommended in this case?)

[RubenImhoff]

Besides the former questions, good to go from my side, @mats-knmi, great work!

[mats-knmi]

Just to be sure @RubenImhoff, did you check that the other locations in the blending code that use zero_precip_radar are not harmed by this change. I am not sure what these parts of the code do and if they would be negatively impacted:

pysteps/pysteps/blending/steps.py

Lines 1333 to 1340 in 76539bb

	# If there are values in the radar fields, compute the auto-correlations
	GAMMA = np.empty((self.__config.n_cascade_levels, self.__config.ar_order))
	if not self.__params.zero_precip_radar:
	# compute lag-l temporal auto-correlation coefficients for each cascade level
	for i in range(self.__config.n_cascade_levels):
	GAMMA[i, :] = correlation.temporal_autocorrelation(
	self.__state.precip_cascades[i], mask=self.__params.mask_threshold
	)

pysteps/pysteps/blending/steps.py

Lines 2118 to 2126 in 76539bb

	# Make sure we have values outside the mask
	if self.__params.zero_precip_radar:
	precip_forecast_recomp_subtimestep = np.nan_to_num(
	precip_forecast_recomp_subtimestep,
	copy=True,
	nan=self.__params.precip_zerovalue,
	posinf=self.__params.precip_zerovalue,
	neginf=self.__params.precip_zerovalue,
	)

pysteps/pysteps/blending/steps.py

Lines 2160 to 2168 in 76539bb

	# Make sure we have values outside the mask
	if self.__params.zero_precip_radar:
	precip_extrapolated_decomp = np.nan_to_num(
	precip_extrapolated_decomp,
	copy=True,
	nan=np.nanmin(precip_forecast_cascade_subtimestep),
	posinf=np.nanmin(precip_forecast_cascade_subtimestep),
	neginf=np.nanmin(precip_forecast_cascade_subtimestep),
	)

[bwalraven]

Thanks @mats-knmi! This would indeed solve 'my' issue in #451!

Perhaps a 'nice-to-have', with regards to @dnerini's comment about handling zero fields more graciously instead of raising an error I am wondering whether it would be useful to also generate a warning or more informative error message in the case that win_fun is not None, but precip_field * tapering does return a field of zeros?

This might not happen very often now that the adapted tukey function also covers the edges of a more rectangular grid, but in the case that it does happen, it would point the user more directly to the problem without having to dig through the code extensively. Alternatively this could also be done by a "no-rain check" like @mats-knmi proposed in point 2: error handling.

In any case, thanks for proposing a solution that solves the issue!

[RubenImhoff]

Just to be sure @RubenImhoff, did you check that the other locations in the blending code that use zero_precip_radar are not harmed by this change. I am not sure what these parts of the code do and if they would be negatively impacted:

pysteps/pysteps/blending/steps.py

Lines 1333 to 1340 in 76539bb

	# If there are values in the radar fields, compute the auto-correlations
	GAMMA = np.empty((self.__config.n_cascade_levels, self.__config.ar_order))
	if not self.__params.zero_precip_radar:
	# compute lag-l temporal auto-correlation coefficients for each cascade level
	for i in range(self.__config.n_cascade_levels):
	GAMMA[i, :] = correlation.temporal_autocorrelation(
	self.__state.precip_cascades[i], mask=self.__params.mask_threshold
	)

pysteps/pysteps/blending/steps.py

Lines 2118 to 2126 in 76539bb

	# Make sure we have values outside the mask
	if self.__params.zero_precip_radar:
	precip_forecast_recomp_subtimestep = np.nan_to_num(
	precip_forecast_recomp_subtimestep,
	copy=True,
	nan=self.__params.precip_zerovalue,
	posinf=self.__params.precip_zerovalue,
	neginf=self.__params.precip_zerovalue,
	)

pysteps/pysteps/blending/steps.py

Lines 2160 to 2168 in 76539bb

	# Make sure we have values outside the mask
	if self.__params.zero_precip_radar:
	precip_extrapolated_decomp = np.nan_to_num(
	precip_extrapolated_decomp,
	copy=True,
	nan=np.nanmin(precip_forecast_cascade_subtimestep),
	posinf=np.nanmin(precip_forecast_cascade_subtimestep),
	neginf=np.nanmin(precip_forecast_cascade_subtimestep),
	)

I have not tested it on real data yet (I would need to have a couple of good rain-no rain cases for that), but from what I can see, it will mainly change when we have or don't have no_rain, which should still give a stable result (and probably goes to a no rain alternative sooner, which is not bad in these cases, I think). Perhaps @bwalraven, could test it for a couple of his failing cases?

[mats-knmi]

Thanks @mats-knmi! This would indeed solve 'my' issue in #451!

Perhaps a 'nice-to-have', with regards to @dnerini's comment about handling zero fields more graciously instead of raising an error I am wondering whether it would be useful to also generate a warning or more informative error message in the case that win_fun is not None, but precip_field * tapering does return a field of zeros?

This might not happen very often now that the adapted tukey function also covers the edges of a more rectangular grid, but in the case that it does happen, it would point the user more directly to the problem without having to dig through the code extensively. Alternatively this could also be done by a "no-rain check" like @mats-knmi proposed in point 2: error handling.

In any case, thanks for proposing a solution that solves the issue!

The solution I proposed makes it so that if precip_field * tapering would be a field of zeros, then zero_precip_radar will be True. This means that the code that throws the error will not be executed. This should be enough for this error to not come back.

[bwalraven]

I have not tested it on real data yet (I would need to have a couple of good rain-no rain cases for that), but from what I can see, it will mainly change when we have or don't have no_rain, which should still give a stable result (and probably goes to a no rain alternative sooner, which is not bad in these cases, I think). Perhaps @bwalraven, could test it for a couple of his failing cases?

I am not familiar enough with the blending steps to say something about that, but I can have a look and test it for the simple nowcasting cases I have. But seeing @mats-knmi's resulting image of the adapted tukey function, and especially if the default tapering function becomes None i.e. np.ones() rather than tukey this should solve the issue I encountered when running the steps function.

[mats-knmi]

The default will still be tukey, you can set it to None yourself if you want, but that is not recommended as @dnerini said.

[dnerini]

hi @mats-knmi, thanks a lot for taking this up!

The changes to the tapering windows make a lot of sense to me and definitely improve over the current version, very nice!

For the approach with zero precip, I agree that the best way to handle this is to account for the effect of the tapering window on the norain check. I only left a very small comment concerning the deprecation warning, which I think it could be safely removed.

[dnerini]

I'd suggest to leave win_fun=None as default, so to maintain the same behaviour. The default can then be changed in the call in the blending method

[mats-knmi]

Hmmm, yes but in the fftgenerators the default method is tukey. So if a user uses an ensemble nowcast with default settings and then also the check_norain method with default settings, you would expect it to work. Which it won't if the default is different here as in the fftgenerators. It would be ideal if they could somehow be synced, but I don't have a good idea for that right now.

[dnerini]

agreed, but in principle check_norain is indipendent of the settings in the fftgenerators. So it'll be important to use the right win_fun argument in the nowcasting methods, but in itself I would still set the default to None here. Or am I missing on something?

[mats-knmi]

The reason I set the tapering method to tukey by default is because this check_norain method currently only exists for the purpose of running the fftgenerators. Since for the fftgenerators the tukey method is the default, I thought it would make sense to have it as the default here as well.

Ofcourse the check_norain method is a quite generic method that could be used for other purposes as well, in which case a default of None would make more sense.

I think that if we call check_norain from the nowcast methods, so that the user doesn't have to do that anymore, it should be fine to set the default to None. But then we would have to update the error message in the fftgenerators as well to indicate to a user like @SKB-7 (who uses the fftgenerators directly) that in order to successfully check for no rain in the context of the fftgenerators, you will have to call check_norain with the correct tapering method (which is tukey by default for the fftgenerators).

[mats-knmi]

I think the best approach is to do what I typed in the last paragraph of my previous message, so I will implement that today.

[dnerini]

@mats-knmi commit ee9a787 adds a test for the steps nowcast method with all zero inputs (it will fail), as a reminder that something needs to be still done for the use case, as discussed in #451.

[mats-knmi]

@dnerini Do you want me to add the check_norain to the beginning of the steps nowcast? This way it should not go wrong anymore for new users?

And what should I do with the other 2 nowcast methods that use noise generators? I see linda only uses pbs, but sseps could use fft generated noise as well right?

[dnerini]

@dnerini Do you want me to add the check_norain to the beginning of the steps nowcast? This way it should not go wrong anymore for new users?

yes, I think this could be a good approach. What do you do with the blending? return an empty, all zero output?

And what should I do with the other 2 nowcast methods that use noise generators? I see linda only uses pbs, but sseps could use fft generated noise as well right?

in principle yes, but maybe a better approach is to include a dedicated test for these methods too, then we'll need to include the check only for those methods that fail?

[mats-knmi]

@dnerini I saw your new tests, thanks. I will implement the no rain check in all the nowcast methods that fail these tests and that should fix them. I will return only zeros when the norain check fails indeed.

[dnerini]

Note that currently the failing test_plugins_support.py is outside the scope of this PR. I believe this is caused by changes in https://github.com/pySTEPS/cookiecutter-pysteps-plugin that will be addressed in #418. Maybe @ladc can quickly confirm that this is indeed the case?

[mats-knmi]

@dnerini I think I have fixed everything now, can you review the final version? If you are on board I would like to merge this and release this as a new version, since we want to start using this operationally. How do I do this? Do I just bump the version in the setup.py file and then merge this pull request?

[dnerini]

Looks very good to me! Left only two super minor comments, then I believe this is ready to be merged

[dnerini]

why remove this note? this is still a very experimental method ...

[mats-knmi]

I don't know how this happened. I did not mean to remove this part also I did not mean to add this line

pysteps.extrapolation.interface, pystepsemilagrangian

I am pretty confused, but I re-added this note.

[mats-knmi]

@dnerini I will merge this then, how do I turn this into a release? Or is there something I need to wait on before doing that?

[dnerini]

I would like to merge this and release this as a new version, since we want to start using this operationally. How do I do this? Do I just bump the version in the setup.py file and then merge this pull request?

So what I typically do:

1. merge PR
2. bump version directly on master (see e.g ee60fa6)
3. release new version (use automatic release notes). In this case, we probably want to release a new minor version, ie v1.15.0
4. check that the release github actions succesfully completes, eg https://github.com/pySTEPS/pysteps/actions/runs/13141012118

do you rely on the conda-forge installation? this typically takes few extra days (if everything goes well)

[mats-knmi]

Thanks, I'll try that!

do you rely on the conda-forge installation? this typically takes few extra days (if everything goes well)

No we install pysteps through pip

[mats-knmi]

@dnerini I cannot push directly to master (which is a good thing imo), but this means I cannot bump the version directly. So either you would have to do this for me or you could approve this PR: #462.

So maybe next time I bump the version directly in the feature PR or how do we do this?

[RubenImhoff]

Great work!!

[ladc]

Note that currently the failing test_plugins_support.py is outside the scope of this PR. I believe this is caused by changes in https://github.com/pySTEPS/cookiecutter-pysteps-plugin that will be addressed in #418. Maybe @ladc can quickly confirm that this is indeed the case?

Yes, the merge request fixing this is #405.