Snowfall and rain approximation from Dai 2008

[aulemahal]

Pull Request Checklist:

• This PR addresses an already opened issue (for bug fixes / features)
  • This PR fixes Add additionnal method to snowfall_approximation #1752
• Tests for the changes have been added (for bug fixes / features)
  • (If applicable) Documentation has been added / updated (for bug fixes / features)
• CHANGELOG.rst has been updated (with summary of main changes)
  • Link to issue (:issue:number) and pull request (:pull:number) has been added

What kind of change does this PR introduce?

• New methods "dai_annual" and "dai_seasonal" for snowfall and rain approximation.
• New argument landmask to control which formulation (land or ocean) we use.

Does this PR introduce a breaking change?

No. Only new arguments added, no default changed.

Other information:

@Zeitsperre : The landmask argument is polymorphic (my favorite, your nemesis (?)), it can either be a DataArray (which is assumed to be of boolean dtype, but no need to check) or a straight boolean. This is not a Quantified because boolean have no units. Right now, I put the annotation that resolved into the "optional variable" InputType, which means "bool" is not part of it. Should we create a new InputType for boolean mask ? It would be similar to Quantified in that it can be a scalar or a DataArray, but without the units check.

@tlogan2000 : I implemented the rain approximation as well. However, in opposition to the 3 other methods already implemented, the Dai rain approximation is not equal to pr - prsn because it considers some small fractions of the total precip to be of another phase (denoted "sleet" in the article). Thus, to imitate DonnéesClimatiques, ones need to compute prra as pr - prsn instead of using rain_approximation. Is this ok ?

Tests need to be added, I'll try to base my tests off of DonnéesClimatiques' data. At least it would compare xclim with another implementation, if not an "original" one.

[aulemahal]

Hi @kkkchow, I have a few questions about the prsn computed for CanDCS-M6 and available on pavics.

You may find my implementation of the snowfall fraction following Dai 2008 in this PR. However, when I computed prsn using pr and tas from the datasets linked above, I found some discrepancies with the version that's already there.

I used a randomly selected simulation : CanESM5/historical, over a domain including QC, On and part of the maritimes. I only used data for a single winter : 2000-12 to 2001-03 (but the analysis is the same over other winters).

I computed different flavours of the new method "Dai", with and without the clipping you described in your comment on the issue. I also tried different values for the clipping. The best fit between my implementation and yours is when I use 4°C (meaning snow fraction is 1 under -4°C and 0 above 4°C). I also computed prsn with the "binary" method for reference.

The next figure shows the total snow fall in mm of SWE over the period. The difference are shown as "CanDCS-M6 minus xclim", so red means my implementation underestimates prsn in comparison to CanDCS-M6. I also added the mean absolute error in the title.

This other figure looks at the distribution of this difference with temperature. You are looking at 2D histogram with a log norm on the colors of the difference vs temperature. I added more clipping values so we can see the behaviour.

As we can notice, the CanDCS-M6 prsn seems to be computed with a clipping of around 5°C ? I verified by computing the mean snowfall and mean rain (over all grid and all period) for 0.2°C bins of temperature from -10.1 to 10.1 °C and it gives the next figure:

The rain goes to 0 around -3°C and the snowfall goes to 0 around 5°C. (When I say "0" I mean 6.215e-5 kg m-2 d-1 because that's the smallest value of pr in the dataset. Probably some packing artefact.)

So with all that the question is: is this discrepancy problematic or does this sound non-significative. To me, an overall 1 mm of prsn error in the full winter sound ok, but still the differences are strange and that lowest error is for a different clipping value than what you said you used.

Question also applies to @tlogan2000 !

[aulemahal]

For the record, the issue denoted above was solved. The way the "clipping" was done was different from CanDCS and my first implementation. This has been fixed.

"clipping" is now a rescaling of the fraction function. The original tanh function from the paper is not actually bounded by 0 and 1, there's a small snow fraction even at very high temperatures (and vice-versa). The idea from the UVic and CCCS team is to rescale the function to 0 and 1 at a certain clip_temp. This idea has been implemented here.

[aulemahal]

The difference between this implementation and CanDCS-M6 is now insignificant. In my test, the total accumulation of snow in one winter differs (in the mean) by 0.01 mm (of SWE).

[Zeitsperre]

@Zeitsperre : The landmask argument is polymorphic (my favorite, your nemesis (?)), it can either be a DataArray (which is assumed to be of boolean dtype, but no need to check) or a straight boolean.

FTR, I don't hate polymorphism; I don't like when polymorphic variables aren't properly annotated or constrained. Live your best (polymorphic) life.

This is not a Quantified because boolean have no units. Right now, I put the annotation that resolved into the "optional variable" InputType, which means "bool" is not part of it. Should we create a new InputType for boolean mask ? It would be similar to Quantified in that it can be a scalar or a DataArray, but without the units check.

I think that makes the most sense. Do you see this being used elsewhere in the library? If so, try to make that annotation accessible to multiple modules.

[aulemahal]

Cool! I would argue that "xr.DataArray | bool" is simple enough it doesn't require a new type ?

Quickly browsing through the other indicators, I didn't find perfect candidate for the new type, but I found these potential ones:

• latitude_temperature_index : arg lat_factor might be on an array ?
• clausius_clapeyron_scaled_precipitation : arg cc_scale_factor, I don't know how an approximation this one is.
• potential_evapotranspiration: args peta and petb that seem to depend on spatially-sensitive calibration
  All these would be annotated xr.DataArray | float.

[Zeitsperre]

@aulemahal

The merge conflict is relatively annoying given that the conversion indices have been refactored. Let me know if you want me to take over that task.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[aulemahal]

@Zeitsperre I think your PR with the conversion refactor forgot to edit references to the submodules in the doc! The last commit adds convert where it makes sense, and rephrases a few sentences.

[Zeitsperre]

@Zeitsperre I think your PR with the conversion refactor forgot to edit references to the submodules in the doc! The last commit adds convert where it makes sense, and rephrases a few sentences.

Woops, yes, I did indeed forget. Thanks!