Add routines to compute TISR (#194)

* Fix zenith angle time zone handling

python datetime objects do not handle the timezone by default, so
computing day differences, etc can give unreliable answers. This was
leading to the cos zenith angle tests failing on my Linux desktop.

Fixed by explicitly adding the UTC timezone info.

* Add tisr function

* Use more accurate formula for elliptic anomaly

* Skip the tisr test

* Change initial guess to E=M for simplicity

* Delete tisr test

* convert to numpy doc

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):

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",

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