Perez_v3 output issue (SOLWEIG)

[bweeding]

I'm getting incorrect outputs from Perez_v3, which have showed up in kdown and Tmrt. I’m supplying my own diff and dir radiation (from the Australian Bureau of Meteorology’s BARRA reanalysis product). I’ve traced the errors all the way back to the Perez_v3 function in Solweig_2021a_calc:

lv = Perez_v3(zenDeg, azimuth, radD, radI, jday, patchchoice) # Relative luminance

The inputs for Perez_v3 that produce the bad results are:

zendeg: 85.7476675536682 azimuth: 67.36351431750992 radD: 73.1 radI: 100.1 jday: 229.0 patchchoice: 1

This produces an output with the 3rd column hundreds of times larger than for the other normal outputs:

Which then goes into the anisotropic luminosity:

aniLum = aniLum + diffsh[:,:,idx] * lv[0][idx][2]

dRad = aniLum * radD

Causing dRad to have extremely large values:

max: 3176.9642980328917 min: -2196.1060734137354 med: 1078.274174716783

compared to a normal dRad: max: 101.12 min: 0.15784201699530884 med: 50.55301209344854

Which is then fed into kdown. This results into large kdown and Tmrt values. I also traced an occurrence of low Tmrt and kdown values, and found that the same thing happened, though at 2pm in the afternoon!

zendeg: 22.940787488281416 azimuth: 355.04264645985876 radD: 110.58 radI: 87.59 jday: 21.0 patchchoice: 1 [ 6. 36. 9.76427633] [ 6. 36. 9.76427633] [ 6. 36. 9.76427633] [ 6. 36. 9.76427633] [ 6. 36. 9.76427633] [ 6. 36. 9.76427633] [ 6. 36. 9.76427633]

In both of these cases, extreme kdown led to extreme Tmrt.

[biglimp]

Thanks Ben. We should try to identify when this happens. It would be good if you can try to identify a pattern, e.g. when radI and radD is very similar or at high zenith angles etc. I would be good to also know radG at these occasions, as well as times of day.

Can you also proved an inputexample when this happens, e.g. all info needed (metfile, location etc.)

[bweeding]

I'll have a go at some deeper analysis first thing on Monday!

On Fri, 4 Feb 2022, 21:02 biglimp, @.***> wrote:

Thanks Ben. We should try to identify when this happens. It would be good if you can try to identify a pattern, e.g. when radI and radD is very similar or at high zenith angles etc. I would be good to also know radG at these occasions, as well as times of day.

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1029822946, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRBR44OG3WJUYSQOWA3UZOP3RANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[nilswallenberg]

Hi Ben,

Do you know if your input radI is for a horizontal surface or perpendicular to sun's beam? It should be perpendicular to the sun's beam. If the input is on a horizontal surface it is easy to change it to perpendicular with one line of code.

radI[i] = radI[i]/np.sin(altitude[0][i] * np.pi/180)

This line should be then be added to the solweig_algorithm.py just before so.Solweig_2021a_calc (approximately line 900 for UMEP processing and 300 for UMEP).

Cheers, Nils

[bweeding]

Hi Nils, in addressing the issue I think I found another one at my end which I have solved. I'll deal with that first.

I realised that because I'm using hourly mean data, I should be using timestamps on the half hour, so that I get an average sun position from the model as well. Therefore my data which is the average from 1300 to 1400 is stamped 1330 instead of 1300. You can see in the figure below that moving from stamping data on the hour compared to on the half hour improves things significantly:

I then tested out the modified radiation formula, as I'm fairly sure my radiation data is normal to the ground, as it comes from a surface level parameter in Wm^-2:

As you can see, the addition of the radiation formula introduces a range of strange results.

I'm going to run the same tests with some longer times periods and will send through the results when they're done.

Cheers Ben

[bweeding]

Investigation of strange results

I've produced two runs of a 50m by 50m tile from 2016 through 2018 using hourly mean reanalysis data.

As you can see, the seasonality in Tmrt is well represented, although there are two timestamps with extremely low average Tmrts that are clearly errors. Adding in the radiation modification as suggested by Nils reduces some of the extreme values, and adds an interesting peak in high but reasonable Tmrt values.

We can then look and the maximum and minimums for each run: While the minimums both appear sensible, the maximum values on the modified radiation run are quite strange, peaking in winter. The radiation modified run also contains far more extreme value spikes. Given this, I've focused the rest of this analysis on the default setting.

The next plot shows azimuth, altitude, radD, and radI plotted with Tmrt: This shows the extreme errors occurring at lower radiation values and sun altitudes, and wider azimuths. To clarify the times of year these errors occur at we can use a 2D histogram of Tmrt and day of year: We have non during Summer, but some occurring in all other seasons.

To check if a relationship appears with the ratio of altitude to azimuth and Tmrt: This shows extreme values only occur when the ratio is small and positive.

Doing the same with Rad I and Rad D: This shows the extreme values occur when Rad D is relatively large, in the context of both radiation values being small.

I haven't continued any deeper down the rabbit hole, as I'd like to hear if you have any guiding comments or ideas before I go any further?

Cheers Ben

[nilswallenberg]

There are a few things you could do or check:

1. Compare global radiation with I0 to check if you are off in time (does the reanalysis correct for daylight saving, etc, and how it is treated in your SOLWEIG run?). If you compare global radiation and I0 on a clear day (one in winter and one in summer) they should more or less match.

2. The example from 2PM on January 21st: zendeg: 22.940787488281416 azimuth: 355.04264645985876 radD: 110.58 radI: 87.59 jday: 21.0

Is it realistic to have only ~200 radG at this specific time? If it is overcast, radD should probably be higher if radI is 87. If it was fully overcast, radI should be zero. Could be something wrong with the input data here.

3. Another possible explanation is that the Perez model is not able to solve for some specific inputs. It is after all an empirical model.

4. Also found this (nr 8 in the following pdf): http://www.bom.gov.au/research/projects/reanalysis/Known_issues_in_BARRA_dataset_July_2019.pdf Maybe it helps, not sure which one of the BARRA products it relates to.

5. There are quite few errors, which you could filter out. :-)

Hope we can solve this! :-)

Cheers, Nils

[bweeding]

Hi Nils, I'm sure we can solve it!

I checked the I0 and global and there was indeed a mismatch occurring during winter! This was due to me using a local time conversion which included daylight savings. I fixed this by using a set offset of 11 hours:

I then ran the 2016 through 2018 period again: I only received one warning message during processing:

We now only have two timestamps with erroneous values, shown in the top two panels, with the following hour below them: The large white shape is a building where I've eliminated the Tmrt, while the speckled white areas must be NaNs produced by the Tmrt algorithm I'm assuming?

Given the uncertainty in the BARRA radiation and the empirical nature of the formula, 2 bad hours in 3 years of data seems pretty good! I'm going to now run the full 30 year time series over the weekend and will then produce plots on Monday.

Cheers Ben

[biglimp]

Great Ben! Looks like you sorted out almost all strange values. Is the last maps is from a time step where you have strange data? If everything is ok, no NaNs should be present.

[bweeding]

Thanks! The last maps are normal timestamps, but when I calculate UTCI using xarrays, I also set the Tmrt values of all the building locations to NaN so I'm only dealing with pedestrian accessible data for analysis.

[bweeding]

Hi Fredrik and Nils, after sorting out the time alignment issue, the data has cleared up significantly. However, there are still some "bad" timestamps. I've pulled out images of timestamps with Tmrt less than -20 and plotted them here. I'm going to do a bit more digging, but just wondering if there are any patterns or issues you recognise?

Cheers Ben

[biglimp]

Hm. I see some linear features in your output (about 526900). What are these and could they affect the output?

[bweeding]

I believe these are artifacts from the combining of the tiles. The whole 150m x 150m image is made from nine 50m x 50m tiles. Each of the 50m x 50m tiles is processes as the centre of its own 150m x 150m tile to account for tall building shadows. The linear features only seem to appear at low sun angles (dusk/dawn), so I think they must be from differences in the buildings that are used to calculate shadows. Does that make sense?

[biglimp]

Me theory is that at certain, very specific times during dawn and/or dusk, the calculation of sun angle and the observed Kdown does not match and somehow creates these strange Tmrt values. That could be if e.g. have positive sun angle but still no positive Kdown or something like that. Could you take one of the strange time steps and see what the full uotput looks like by putting in a POI somewhere in the sunlit part of the model. Then we could see from which flux the large negative values come from (Kdown, Kup, Kside, Ldown, Lup or Lside).

[bweeding]

Ah ok makes sense, I'll do it first thing Monday. I can produce maps of all the outputs using the xarray modification I've made.

On Fri, 4 Mar 2022, 7:13 pm biglimp, @.***> wrote:

Me theory is that at certain, very specific times during dawn and/or dusk, the calculation of sun angle and the observed Kdown does not match and somehow creates these strange Tmrt values. That could be if e.g. have positive sun angle but still no positive Kdown or something like that. Could you take one of the strange time steps and see what the full uotput looks like by putting in a POI somewhere in the sunlit part of the model. Then we could see from which flux the large negative values come from (Kdown, Kup, Kside, Ldown, Lup or Lside).

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1058934524, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRHHWLR456ZWMEJCRDDU6HAZLANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[biglimp]

It would also be good to get the sun angles etc. for that time step. Hence the POI.

[bweeding]

Ah yep no worries, I'll put it all together. Have a good weekend!

On Fri, 4 Mar 2022, 7:26 pm biglimp, @.***> wrote:

It would also be good to get the sun angles etc. for that time step. Hence the POI.

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1058944082, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRGJGIZIK4M2263C7NLU6HCMXANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[bweeding]

Here's a POI from the model - it's at 526995, 5252145.

There's a negative value for kdown - could that be it?

POI_1_for_FL.txt l

[biglimp]

Well, that is just a consequence of something. I see that you are using the wrong Kdir. It should be direct beam perpedicular to the solar position. Now you are using Kdir for a horisontal surface. I see this since Kdir+Kdiff = Kglobal. This could also be the reason why you get negative Kdown.

Kdir, Kdiff and Kdown are input parameters whereas the other K-components are model outputs.

Recalculate Kdir direct-beam radiation = (global−diffuse) / sin(sun altitude) and try again.

[bweeding]

Ah yes, this is something Nils mentioned before but I tried and the results were odd.

So I modified kdir by inserting: radI[i] = (radG[i]-radD[i])/np.sin(altitude[0][i] * np.pi/180) into solweig_algorithm.py just before the use of so.Solweig_2021a_calc.

Looking first at the 1 day time period from the POI file, modifying the radiation makes the Tmrt output more reasonable, though the numbers are still high in some locations:

If we look at the spatial maximum for each hour we see modifying the radiation has quite an impact in increasing the Tmrt values:

and if we look at the radiation values from the POI file, we can see increases in kdir and kdown:

Checking to see what happens on a longer time period, we see that using the radiation modification results in much higher Tmrt values, including some very large spikes.

If we do an hourly boxplot of all the data (not spatial max/min), we can see these extreme values occur at the 7:30, 8:30, and 9:30 timestamps.

If I perform the same plot without the radiation modification, we see strange values for 8:30 and 15:30, though only low rather than high:

Does that suggest any particular issue other than extremely high values cased by dividing by sin(alt) when alt is very low?

Cheers Ben

[biglimp]

The values looks more reasonable now when you are using the correct Kdir.

I think you are correct that the division of sin(alt) is one of the issues here. The model is very sensitive during mornings and evenings and if you have small time biases between sun position and the observed forcing data you can end up with these spikes. It is always the relation between sun altitude and radiation that can cause this. I suggest that you filter out values when you have:

• larger Kdiff than Kglobal
• Negative Kdown
• Unreasonable high Kdir (larger than e.g. 1500)

Me and Nils will think on how to get rid of these unreasonable high/low Tmrt values by putting in restrictions in the model. We might come back to you to test later.

[bweeding]

Ok great I'll have a tinker tomorrow and show you the results. Should I stick with the radiation modification though?

On Tue, 8 Mar 2022, 7:05 pm biglimp, @.***> wrote:

The values looks more reasonable now when you are using the correct Kdir.

I think you are correct that the division of sin(alt) is one of the issues here. The model is very sensitive during mornings and evenings and if you have small time biases between sun position and the observed forcing data you can end up with these spikes. It is always the relation between sun altitude and radiation that can cause this. I suggest that you filter out values when you have:

• larger Kdiff than Kglobal
• Negative Kdown
• Unreasonable high Kdir (larger than e.g. 1500)

Me and Nils will think on how to get rid of these unreasonable high/low Tmrt values by putting in restrictions in the model. We might come back to you to test later.

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1061509305, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRERO6X2MF4FUFA67EDU64C57ANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[biglimp]

Yes, radiation perpendicular to the solar beam.

[bweeding]

Ok great, I'll send through some results tomorrow! Enjoying figuring this out!!

On Tue, 8 Mar 2022, 8:02 pm biglimp, @.***> wrote:

Yes, radiation perpendicular to the solar beam.

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1061552944, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRAXLWSDMFDXZTIETZDU64JTJANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[bweeding]

So I included the following conditions to filter out the "bad" values:

          if (radDout>radG[i]).any():
                Tmrt.fill(np.nan)

            if (Kdown < 0).any():
                Tmrt.fill(np.nan)

            if (radIout > 1500).any():
                Tmrt.fill(np.nan)

and used data from starting_date = '1994-01-01'; finishing_date = '1996-12-31'

This gives reasonable looking mins and maxes:

If we do an hourly boxplot, it seems like there are some very high values in the morning:

The filtering removed 847 timestamps out of 26281, ~3.2%. If we plot the filtered timestamps by hour, we see morning contributes more than evening:

I can translate the good Tmrt values over to UTCI:

and then compare Tmrt, Tair, and UTCI:

[bweeding]

In relation to the problems caused by offsets between sunrise in altitude and in the meteorology, will that be connected to the position of the target location within the timezone? Or is that accounted for within the model?

[biglimp]

That is accounted for. I am thinking of other things that could happen e.g.:

• A clock in a logger is showing wrong time when using observed data as input forcing. Very small difference can produce the error that you have discovered
• The data is not average up to the recorded timestep but instead momentaneous values for that timestep

Thanks for presenting your restrictions in the comment above. We will include these in the model in the future.

[bweeding]

Ah ok. My data is averaged for the surrounding hour ie. 0930 data is the mean from 0900 to 1000. Could that be it?

On Wed, 9 Mar 2022, 6:35 pm biglimp, @.***> wrote:

That is accounted for. I am thinking of other things that could happen e.g.:

• A clock in a logger is showing wrong time when using observed data as input forcing. Very small difference can produce the error that you have discovered
• The data is not average up to the recorded timestep but instead momentaneous values for that timestep

Thanks for presenting your restrictions in the comment above. We will include these in the model in the future.

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1062632842, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRFHVEY64MA6UWVYFCLU7BIEBANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[biglimp]

Well, it depends where your data is coming from. If you have observed averaged data for 0900 and 1000 etc., you should not move the data to 0930. AND, if you have momentaneous data you should use 0930. SOLWEIG calulates the sun position inbetween timesteps, i.e. if you have data on the hour e.g. 0900 and 1000 then the position is calculated for 0930.

With obsereved avarge data I mean that fro 1000, an average of samples between 0900 and 1000 is recorded for the 1000 timestep.

Hope this makes sense...

[bweeding]

Ah so if my data is an average from 0900 to 1000 (it's a reanalysis product) then I should associate it with a timestamp of 1000?

On Wed, 9 Mar 2022, 7:28 pm biglimp, @.***> wrote:

Well, it depends where your data is coming from. If you have observed averaged data for 0900 and 1000 etc., you should not move the data to 0930. AND, if you have momentaneous data you should use 0930. SOLWEIG calulates the sun position inbetween timesteps, i.e. if you have data on the hour e.g. 0900 and 1000 then the position is calculated for 0930.

With obsereved avarge data I mean that fro 1000, an average of samples between 0900 and 1000 is recorded for the 1000 timestep.

Hope this makes sense...

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1062670935, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRF2TI2YUZLAKZ3NK4TU7BOLPANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[biglimp]

Yes, if reanalysis data timestamp 1000 is an average between 0900 and 1000.

[bweeding]

Ok I'll restamp the times tomorrow and send through a comparison.

On Wed, 9 Mar 2022, 8:18 pm biglimp, @.***> wrote:

Yes, if reanalysis data timestamp 1000 is an average between 0900 and 1000.

— Reply to this email directly, view it on GitHub https://github.com/UMEP-dev/UMEP/issues/365#issuecomment-1062713096, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANYGKRFNDUFPYOXCPMIGCRTU7BUE7ANCNFSM5NRIFVIQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[bweeding]

Ok I've used 1994 data to compare the old (on the half hour) to the new (at the end of the hour) timestamping. This reduced the number of "bad" timestamps from 266 to 178 (from 3% to 2%).

We can see a shift in the Tmrt values and what looks like a reduction in extreme values.

Checking this with a histogram confirms that an actual change in values has occurred, rather than simply a shift of values from one hour to another:

[biglimp]

@bweeding , can i close this (or maybe move to discussions)?

[bweeding]

Yep, close away!