diff --git a/modulus/utils/zenith_angle.py b/modulus/utils/zenith_angle.py index 2bca1f2b5a..c03b27130d 100644 --- a/modulus/utils/zenith_angle.py +++ b/modulus/utils/zenith_angle.py @@ -32,12 +32,16 @@ import numpy as np +# can replace this import with zoneinfo from the standard library in python3.9+. +import pytz + # helper type dtype = np.float32 + T = TypeVar("T", np.ndarray, float) -TIMESTAMP_2000 = datetime.datetime(2000, 1, 1, 12, 0).timestamp() +TIMESTAMP_2000 = datetime.datetime(2000, 1, 1, 12, 0, tzinfo=pytz.utc).timestamp() def cos_zenith_angle( @@ -51,12 +55,16 @@ def cos_zenith_angle( be assumed to be in degrees, unless they have a units attribute that contains "rad"; in that case they will automatically be converted to having units of degrees. - Args: - time: time in UTC - lon: float or np.ndarray in degrees (E/W) - lat: float or np.ndarray in degrees (N/S) - Returns: - float, np.ndarray + + Parameters + ---------- + time: time in UTC + lon: float or np.ndarray in degrees (E/W) + lat: float or np.ndarray in degrees (N/S) + + Returns + -------- + float, np.ndarray Example: -------- @@ -67,7 +75,7 @@ def cos_zenith_angle( """ lon_rad = np.deg2rad(lon, dtype=dtype) lat_rad = np.deg2rad(lat, dtype=dtype) - julian_centuries = _days_from_2000(time) / 36525.0 + julian_centuries = _datetime_to_julian_century(time) return _star_cos_zenith(julian_centuries, lon_rad, lat_rad) @@ -101,6 +109,200 @@ def cos_zenith_angle_from_timestamp( return _star_cos_zenith(julian_centuries, lon_rad, lat_rad) +def irradiance( + t, + S0=1361, + e=0.0167, + perihelion_longitude=282.895, + mean_tropical_year=365.2422, + newton_iterations: int = 3, +): + """The flux of solar energy in W/m2 towards Earth + + The default orbital parameters are set to 2000 values. + Over the period of 1900-2100 this will result in an error of at most 0.02%, + so can be neglected for many applications. + + Parameters + ---------- + t: linux timestamp + S0: the solar constant in W/m2. This is the mean irradiance received by + earth over a year. + e: the eccentricity of earths elliptical orbit + perihelion_longitude: spatial angle from moving vernal equinox to perihelion with Sun as angle vertex. + Perihelion is moment when earth is closest to sun. vernal equinox is + the longitude when the Earth crosses the equator from South to North. + newton_iterations: number of iterations for newton solver for elliptic anomaly + + + Notes + ----- + + TISR can be computed from Berger's formulas: + + Berger, A. (1978). Long-Term Variations of Daily Insolation and + Quaternary Climatic Changes. Journal of the Atmospheric Sciences, + 35(12), 2362–2367. + https://doi.org/10.1175/1520-0469(1978)035<2362:LTVODI>2.0.CO;2 + + NASA Example computing the orbital parameters: https://data.giss.nasa.gov/modelE/ar5plots/srorbpar.html. From 1900-2100 these are the ranges:: + Orbital Parmameters + + Long. of + Year Eccentri Obliquity Perihel. + (A.D.) city (degrees) (degrees) + ------ -------- --------- -------- + 1900 0.016744 23.4528 281.183 + 2000 0.016704 23.4398 282.895 + 2100 0.016663 23.4268 284.609 + + """ + seconds_per_solar_day = 86400 + mean_tropical_year = mean_tropical_year * seconds_per_solar_day + + year_2000_equinox = datetime.datetime(2000, 3, 20, 7, 35, tzinfo=pytz.utc) + + # from appendix of Berger 1978 + M = (t - year_2000_equinox.timestamp()) % mean_tropical_year + M = M / mean_tropical_year * 2 * np.pi + M -= np.deg2rad(perihelion_longitude) + + # to get the elliptic anomaly E from the "mean anomaly" M + # use eq. 6.37 + # https://link.springer.com/book/10.1007/978-3-662-53045-0) + # r / a = (1 - e cos E ) + # E - e sin(E) = M + def f(E): + return E - e * np.sin(E) - M + + def fp(E): + return 1 - e * np.cos(E) + + # newton iterations + # initial guess + E = M + for _ in range(newton_iterations): + E = E - f(E) / fp(E) + + rho = 1 - e * np.cos(E) + return S0 / rho**2 + + +def toa_incident_solar_radiation_accumulated( + t, + lat, + lon, + interval=3600, + S0=1361, + e=0.0167, + perihelion_longitude=282.895, + mean_tropical_year=365.2422, +): + """Approximate ECMWF TISR with analytical formulas + + According to the ECWMF docs, the TISR variable is integrated over the + preceeding hour. Error is about 0.1% different from the ECMWF TISR + variable. + + Parameters + ---------- + t: linux timestamp + lat, lon: latitude and longitude in degrees + interval: the integral length in seconds over which the irradiance is integrated + S0: the solar constant in W/m2. This is the mean irradiance received by + earth over a year. + e: the eccentricity of earths elliptical orbit + perihelion_longitude: spatial angle from moving vernal equinox to perihelion with Sun as angle vertex. + Perihelion is moment when earth is closest to sun. vernal equinox is + the longitude when the Earth crosses the equator from South to North. + + + Returns + ------- + TOA incident solar radiation accumulated from [t-inteval, t] in J/m2 + + Notes + ----- + + We make some approximations: + + The default orbital parameters are set to 2000 values. + Over the period of 1900-2100 this will result in an error of at most 0.02%, + so can be neglected for many applications. + + The irradiance is constant over the ``interval``. + + From ECWMF [docs](https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation#ERA5:datadocumentation-Meanrates/fluxesandaccumulations) + + Such parameters, which are only available from forecasts, have undergone particular types of statistical processing (temporal mean or accumulation, respectively) over a period of time called the processing period. In addition, these parameters may, or may not, have been averaged in time, to produce monthly means. + + The accumulations (over the accumulation/processing period) in the short forecasts (from 06 and 18 UTC) of ERA5 are treated differently compared with those in ERA-Interim and operational data (where the accumulations are from the beginning of the forecast to the validity date/time). In the short forecasts of ERA5, the accumulations are since the previous post processing (archiving), so for: + + reanalysis: accumulations are over the hour (the accumulation/processing period) ending at the validity date/time + ensemble: accumulations are over the 3 hours (the accumulation/processing period) ending at the validity date/time + Monthly means (of daily means, stream=moda/edmo): accumulations have been scaled to have an "effective" processing period of one day, see section Monthly means + Mean rate/flux parameters in ERA5 (e.g. Table 4 for surface and single levels) provide similar information to accumulations (e.g. Table 3 for surface and single levels), except they are expressed as temporal means, over the same processing periods, and so have units of "per second". + + Mean rate/flux parameters are easier to deal with than accumulations because the units do not vary with the processing period. + The mean rate hydrological parameters (e.g. the "Mean total precipitation rate") have units of "kg m-2 s-1", which are equivalent to "mm s-1". They can be multiplied by 86400 seconds (24 hours) to convert to kg m-2 day-1 or mm day-1. + Note that: + + For the CDS time, or validity time, of 00 UTC, the mean rates/fluxes and accumulations are over the hour (3 hours for the EDA) ending at 00 UTC i.e. the mean or accumulation is during part of the previous day. + Mean rates/fluxes and accumulations are not available from the analyses. + Mean rates/fluxes and accumulations at step=0 have values of zero because the length of the processing period is zero. + + """ # noqa + lat = np.deg2rad(lat) + lon = np.deg2rad(lon) + + century = _timestamp_to_julian_century(t) + ra, dec = _right_ascension_declination(century) + interval_radians = interval / 86400 * 2 * np.pi + # 0 <= h1 < 2 pi + h1 = _local_hour_angle(century, lon, ra) + h0 = h1 - interval_radians + A = np.sin(lat) * np.sin(dec) + B = np.cos(lat) * np.cos(dec) + + # assume irradiance is constant over the interval + S = irradiance(t, S0, e, perihelion_longitude, mean_tropical_year) + sec_per_rad = 86400 / (2 * np.pi) + return S * _integrate_abs_cosz(A, B, h0, h1) * sec_per_rad + + +def _integrate_abs_cosz(A, B, h0, h1): + """Analytically integrate max(A + B cos(h), 0) from h=h0 to h1""" + + hc = np.arccos(-A / B) + + def integrate_cosz(left, right): + return A * (right - left) + B * (np.sin(right) - np.sin(left)) + + def integrate_abs_cosz_from_zero_to(a): + root1 = -hc + 2 * np.pi + T = np.pi * 2 + + # how many periods + n = a // T + + # if there is a root + a = a % T + C = integrate_cosz(0, np.where(a < hc, a, hc)) + D = np.where(root1 < a, integrate_cosz(root1, a), 0) + total = integrate_cosz(0, hc) + integrate_cosz(root1, 2 * np.pi) + return C + D + total * n + + return np.where( + np.isnan(hc), + np.maximum(integrate_cosz(h0, h1), 0), + integrate_abs_cosz_from_zero_to(h1) - integrate_abs_cosz_from_zero_to(h0), + ) + + +def _datetime_to_julian_century(time: datetime.datetime) -> float: + return _days_from_2000(time) / 36525.0 + + def _days_from_2000(model_time): """Get the days since year 2000. @@ -110,13 +312,16 @@ def _days_from_2000(model_time): >>> _days_from_2000(model_time) 731.0 """ + if isinstance(model_time, datetime.datetime): + model_time = model_time.replace(tzinfo=pytz.utc) + date_type = type(np.asarray(model_time).ravel()[0]) if date_type not in [datetime.datetime]: raise ValueError( f"model_time has an invalid date type. It must be " f"datetime.datetime. Got {date_type}." ) - return _total_days(model_time - date_type(2000, 1, 1, 12, 0)) + return _total_days(model_time - date_type(2000, 1, 1, 12, 0, 0, tzinfo=pytz.utc)) def _total_days(time_diff): diff --git a/pyproject.toml b/pyproject.toml index 2a7f56d93e..da20f408b0 100755 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,6 +21,7 @@ dependencies = [ "nvidia_dali_cuda110>=1.16.0", "setuptools>=67.6.0", "certifi>=2023.7.22", + "pytz>=2023.3", ] classifiers = [ "Programming Language :: Python :: 3", diff --git a/test/utils/test_zenith_angle.py b/test/utils/test_zenith_angle.py index a622883cbe..a5051d23e2 100644 --- a/test/utils/test_zenith_angle.py +++ b/test/utils/test_zenith_angle.py @@ -29,10 +29,14 @@ import numpy as np import pytest +from pytz import utc from modulus.utils.zenith_angle import ( + _datetime_to_julian_century, + _timestamp_to_julian_century, cos_zenith_angle, cos_zenith_angle_from_timestamp, + toa_incident_solar_radiation_accumulated, ) @@ -48,6 +52,7 @@ ), ) def test_zenith_angle(time, lon, lat, expected): + time = time.replace(tzinfo=utc) assert cos_zenith_angle(time, lon, lat) == pytest.approx(expected, abs=1e-10) timestamp = time.timestamp() assert cos_zenith_angle_from_timestamp(timestamp, lon, lat) == pytest.approx( @@ -61,3 +66,25 @@ def test_zenith_angle_array(): lon = np.array([0.0])[None, None, :] out = cos_zenith_angle_from_timestamp(timestamp, lon, lat) assert out.shape == (3, 2, 1) + + +@pytest.mark.parametrize( + "t", + [ + datetime(2020, 7, 6, 9, 0, 0), + datetime(2000, 1, 1, 12, 0, 0), + datetime(2000, 7, 1, 12, 0, 0), + datetime(2000, 7, 1, 12, 0, 0, tzinfo=utc), + ], +) +def test_timestamp_to_julian_centuries(t): + a = _datetime_to_julian_century(t) + b = _timestamp_to_julian_century(t.replace(tzinfo=utc).timestamp()) + assert a == b + + +def test_toa(): + t = datetime(2000, 7, 1, 12, 0, 0, tzinfo=utc).timestamp() + lat, lon = 0.0, 0.0 + ans = toa_incident_solar_radiation_accumulated(t, lat, lon) + assert ans >= 0