iridl / python-maprooms

Dash maprooms and tools
0 stars 4 forks source link

eq invalid at lat higher than 55 in Winter months #378

Closed remicousin closed 1 year ago

remicousin commented 1 year ago

equations look all the same in Eunjin python code and in FAO doc (and in my code) however, FAO doc indicates that:

"For the winter months in latitudes greater than 55° (N or S), the equations for Ra have limited validity."

I assume that is where the arccos error came from (I actually could go as far as 70 or 75 before getting it).

aaron-kaplan commented 1 year ago

I think I figured it out: north of the arctic circle and south of the antarctic circle, there are days of the year when the sun doesn't rise or set, so the sunset hour angle is undefined.

I confirmed numerically that at a solar declination of 23.45 (the maximum, which happens on the solstice) and 66.6 degrees latitude (the arctic circle), tan(lat)*tan(solar_declination) comes out to 1.0. At higher latitudes it's greater than 1.

I think the 55 degree thing is something else, because "limited validity" is not the same as "undefined." My guess is that above 55 degrees the thickness of the atmosphere becomes too big a factor to ignore.

Let's leave some notes for the next person who wonders about this. At the least, add the FAO's note about 55 degrees to the docstring. Would it be appropriate also to raise an exception if lat is > 55, or if lat is > 66.6?

remicousin commented 1 year ago

Since this is meant to be applied vectorially on doy and lat, I don't want to raise an exception for the values that should then return nan. But I still get a warning of wrong inputs into arccos... how is that possible?

aaron-kaplan commented 1 year ago

where(cond, x, y) works by first evaluating x and y, and then picking values from them according to cond. So line 500 still evaluates arccos over the entire vector, and the warning is still generated.

I think you want something like

arccos_input = (-1 * np.tan(lat) * np.tan(solar_declination)).where(lambda x: abs(x) <= 1, NaN)
sunset_hour_angle = np.arccos(arccos_input)

(untested)

remicousin commented 1 year ago

@aaron-kaplan please have another look (it crossed my mind that where would evaluation everywhere and then I thought not so)

remicousin commented 1 year ago

@aaron-kaplan this is not urgent but should be trivial to move on.

aaron-kaplan commented 1 year ago

I'll take your word for it. I still haven't worked through all the geometry, and I don't think it's worth my time to do it given that it's of no practical importance.

We could emit a warning (instead of an Exception) if latitudes > 55 or < -55 are used. Also not of great practical importance.