NA values for UTCI - issue with mean radiant temperature calculations?

[emsonali]

Setup Information

• Xclim version: 0.41.0
• Python version: 3.10.10
• Operating System: macOS 12.6.7

Context

I'm calculating UTCI for Southeast Asia and noticed that there are quite a lot of NA values, but this varies spatially and temporally, for parts of the year the northern hemisphere in Southeast Asia has missing values, and other parts of the year the southern hemisphere in Southeast Asia has missing values (see attached screenshots). Could this be due an issue with the mean radiant temperature, as it is influenced by latitude? There's a clear cutoff at the equator.

Code of Conduct

• I agree to follow this project's Code of Conduct

[emsonali]

We plotted csza, fdir_ratio and mrt. It could be that the bug is from fdir_ratio (one of the three conditions that set fdir_ratio to 0)

And here is the corresponding UTCI plot for this time step

[aulemahal]

Hi @emsonali !

This looks like the errors I tried to fix in #1399 this summer. The changes were included in xclim 0.44. However, I kinda focused on the arctic circle-related problems, I might have overlooked issues stemming form the equator... Would you be able to test your code again with the latest xclim version ?

[emsonali]

Hi @aulemahal I checked the environment that we used to run the UTCI calculations and we have xclim 0.45 there. xclim 0.41 was from another environment not used for the calculations, apologies for the mixup

[emsonali]

Hi @aulemahal I reran the calculation again with xclim 0.45 just to be sure and the issue still persists. Please let me know if I can provide further information that might help with your debugging

[aulemahal]

TL;DR Try passing stat='sunlit' to UTCI.

Hi @emsonali, I think I found why you're getting these errors, but if I'm right, that would mean UTCI has been broken for a while... Therefore, I'm invoking @ludwiglierhammer into this thread, as original author of the function. I see that the broken parts were my suggestion / addition, clearly there was something I didn't understand.

The issue has to do with the "stat" argument that UTCI passes to mrt and then to cosine_of_solar_zenith_angle. I don't understand what the default ("average") would mean. When we look at your plot for csza it is smooth, yes, but it contains values going below zero.

The solar zenith angle is the angle the sun makes with the vertical. Thus, at sunset the angle is 90° and at sunrise, 90° again. Thus the cosine can only be negative during the night, when the angle is greater than 90°. csza is the integral of that angle. What you plotted has negative values, so it includes the night period, which is what passing "average" and "sunlit=False" does.

But, fdir requires both the average of the cosine for the "sunlit" interval and for the full interval. The current code, with the default "average", skips that and gives the full interval only. I think stat should be sunlit.

I might have broken that when I fixed the cosine calculations for the arctic circles. Sorry.

@ludwiglierhammer , are you able to confirm my intuition? That "stat='average'" as it is coded now makes no sense and it should be "sunlit" for any requested interval ? The other choice ('instant') is for when we compute the instantaneous value of the MRT.

[ludwiglierhammer]

Hi @emsonali, could you please plot csza and fdir_ratio with both stat="instant" and stat="sunlit". Thus, we can see if @aulemahal's suggestion is right or there are some more bugs in the code. I think this is a good check. As I understand the integral of the solar zenith angle is calculated for both stat="average" and stat="sunlit". The difference is that with stat="sunlit" only the part during sunlit will be considered. Hence, during sunlit stat="average" should be equal to stat="sunlit" (e.g. for subdaily data) and I think we can delete "average" and set "sunlit" to default as @aulemahal did in #1501. Or do you think it makes sence to have the option to calculate an average integral of the angle of each interval considering the non-sunlit periods?

[ludwiglierhammer]

@aulemahal

What exactly does stat="integral" in https://github.com/Ouranosinc/xclim/blob/master/xclim/indices/helpers.py#L194?
I think this is redundant, isn't it?

In the docstrings of both mean_radiant_temperature and universal_thermal_climate_index If "instant", the instantaneous cosine of the solar zenith angle is calculated is mentioned twice.

[emsonali]

@aulemahal @ludwiglierhammer I plotted the czsa and fdir_ratio for stat = 'sunlit' and stat = 'instant' and here are the results - seems like the issue is with czsa_i for stat = 'sunlit'

stat = 'sunlit'

czsa_i

czsa_s

fdir_ratio

stat = 'instant'

czsa

fdir_ratio

[aulemahal]

@ludwiglierhammer I don't remember where I saw a usage of stat='integral' for the cosine... It is what it says, it returns the integral of the angle, the average being the integral divided by the interval length. I changed the default value to average, the one that makes most sense.

I went back in the sources and there's something I don't quite understand. Liljegren 2008 state the equations that we use, but when they speak of fdir, there is no distinction between the csza_i and csza_s. It seems to me that the article only talks about the instantaneous aspect, their equation is to be used to get WBGT at a single moment. Kong 2022 build upon that and there they use averages in a similar way to what we do. Their code (PyWGBT) introduces the distinction between csza_i and csza_s, but I haven't found any explanation of that. The notebook example uses a measured/observed fdir, so that specific function is not even used. (Same with thermofeel : fdir is a given, it is not estimated like here).

The average of the cosine of the solar zenith angle over a day will be smaller than 0 when the night is longer than the day. Haha, this is not a uncommon phenomenon, it happens half of the time ;). csza_i > 0 is what triggers the fdir problem because of the second test here :

xclim/xclim/indices/_conversion.py

Line 1912 in e4d6346

(fdir_ratio <= 0) | (csza_i <= np.cos(89.5 / 180 * np.pi)) | (rsds <= 0),

In Liljegren there is a $\theta < 0$ test , but that's because the computation is for an instant where the sun is set, it makes much sense. I think the test should be modified for the average case.

If I understand correctly, the idea is to mask with 0 if the sun is not over the horizon. In an average idea, that would be : if $max(\cos(\theta)) < 0$ for the interval.

Of we could use the fact that the sunlit average of the cosine is 0 if the sun does not shines in the interval. I thus suggest to simply drop csza_i. And use csza_s in the test.

[ludwiglierhammer]

@aulemahal In the notebook example they mention in block 20:

# the treatments below aim to avoid unrealistically high or low values of f which are also included in Liljegren's code.

Liliegren 2008 talks about instantaneous data and, as you mentioned, they do not differentiate between csza_i and csza_s. Therefor, we could simply drop csza_i and use csza_s in the test. I agree with you since we only consider csza_s in the calculation of fdir. If we wanna calculate mean_radiant_temperature for instananeous data (stat=="instant") csza_i and csza_s are set to the same value in the code. Thus, the test for fdir will stay the same.

stat="integral":
I do not see a different calculation for stat=="integral" and stat=="average" in cosine_of_solar_zentih_angle .
In both cases we call _sunlit_integral_of_cosine_of_solar_zenith_angle with:

   declination,
   lat,
   _wrap_radians(h_ss),
   _wrap_radians(h_s),
   _wrap_radians(h_e),
   stat == "average",

Don't we need a differentation before calling this function?

if stat == "average":
    average=True
elif stat == "integral":
    average=False

And don't we need to call this function with:

   declination,
   lat,
   _wrap_radians(h_ss),
   _wrap_radians(h_s),
   _wrap_radians(h_e),
   average=average,

since we have this:

   if average:
       out = out / denum
   return out

[aulemahal]

Nice!
Last commit to the PR made use of only the sunlit csza.

As for the "integral" case, I think the stat == 'average' (notice the 2 equal signs) does what your if-else proposition does. I'm passing that boolean to average as a positional arg.

@emsonali Are you able to test your code with xclim installed from the fix-mrt branch ? Should be as simple as :

pip install git+https://github.com/ouranosinc/xclim.git@fix-mrt

[ludwiglierhammer]

@aulemahal It was new to me that you can simply pass a boolean as a positional arg by using those 2 equal signs. Thanks for the explanation.

[emsonali]

Hi @aulemahal I tested with the fix-mrt branch and fdir seemed to be produced correctly now with stat="instant"and stat="sunlit" (i.e. where stat="average" and sunlit=False for csza_i and sunlit=True for csza_s). However, for stat="average", it's not working - not sure if this is because sunlit=False for that stat. I've used the stat settings as described in the documentation (see screenshot below), but of course, I have dropped csza_i for computing fdir_ratio as per your fix.

Here are the ouputs with the different stat settings. Please let me know if it's working as you expected, especially for stat="average".

stat="average"
csza

fdir_ratio

stat="instant"
csza_s

fdir_ratio

stat="sunlit"
csza_s

fdir_ratio

[aulemahal]

I'm not sure I fully understand your comment sorry. In the fix-mrt branch, stat='average' doesn't exist anymore for mean_radiant_temperature and utci. The code printscreen you showed is outdated. However, the 'sunlit' mode is the proper one to use (and now the default). The plots you show for that mode seem all right, so I think we can merge the PR!

Thanks!

[emsonali]

Thanks very much @aulemahal @ludwiglierhammer for your help with this!