Laptop changes

setup.cfg

@@ -69,10 +69,11 @@ rst-roles=
     jstor
 docstring-convention=numpy
 doctests=True
-exclude=.git,__pycache__,*_flymake.py,.tox,.eggs,*.egg,build,dist
+exclude=.git,__pycache__,*_flymake.py,.tox,.eggs,*.egg,build,dist,examples
 max-complexity=15
 
 [tool:pytest]
 addopts=
     --doctest-modules --ignore-glob=*_flymake.py
-    --cov --cov-append --cov-report=
+    --cov --cov-append --cov-report=
+    --ignore=examples --ignore=test_fisher

src/olsen_randerson/__init__.py

@@ -3,11 +3,11 @@
 
 from .__version__ import VERSION as __version__  # noqa: F401
 
-NEP_TO_GPP_FACTOR = 2
+NPP_TO_GPP_FACTOR = 2
 """Conversion factor to estimate GPP from NEE
 
 The downscaling needs :abbr:`GPP (Gross Primary Productivity)`, but
-often only :abbr:`NEE (Net Ecosystem Exchange)` is available.  This
+often only :abbr:`NPP (Net Primary Productivity)` is available.  This
 describes how to turn the latter into an estimate of the former.
 """
 Q10 = 1.5
@@ -23,13 +23,22 @@
 """
 
 
-def olsen_randerson_once(flux_nep, temperature, par):
-    """Perform the Olson Randerson downscaling.
+def olsen_randerson_once(
+        flux_npp,
+        flux_rh,
+        temperature,
+        par
+):
+    """Perform the Olson Randerson downscaling for a single month.
 
     Parameters
     ----------
-    flux_nep : np.ndarray[...]
-        Biogenic :abbr:`NEP (Net Ecosystem Productivity)`, usually at
+    flux_npp : np.ndarray[...]
+        Biogenic :abbr:`NPP (Net Primary Productivity)`, usually at
+        monthly timescale.  Units must include time in the
+        denominator.
+    flux_rh : np.ndarray[...]
+        Biogenic :abbr:`Rh (heterotrophic respiration)`, usually at
         monthly timescale.  Units must include time in the
         denominator.
     temperature : np.ndarray[N, ...]
@@ -45,6 +54,7 @@ def olsen_randerson_once(flux_nep, temperature, par):
     -------
     flux_nee : np.ndarray[N, ...]
         The downscaled :abbr:`NEP (Net Ecosystem Productivity)`.
+        (positive is uptake by plants).
 
     References
     ----------
@@ -60,20 +70,29 @@ def olsen_randerson_once(flux_nep, temperature, par):
     >>> par[par < 0] = 0
     >>> temperature = 10 - 10 * np.cos(2 * np.pi * time)
     >>> # First item in row alternates between midnight and noon
-    >>> olsen_randerson_once(np.array(5), temperature, par)
+    >>> olsen_randerson_once(np.array(5), np.array(0), temperature, par)
     array([-3.20043607, -3.4581163 , -4.2353102 ,  4.10768819, 18.33559721,
            23.70071827, 18.33559721,  4.10768819, -4.2353102 , -3.4581163 ,
            -3.20043607, -3.4581163 , -4.2353102 ,  4.10768819, 18.33559721,
            23.70071827, 18.33559721,  4.10768819, -4.2353102 , -3.4581163 ,
            -3.20043607, -3.4581163 , -4.2353102 ,  4.10768819, 18.33559721,
            23.70071827, 18.33559721,  4.10768819, -4.2353102 , -3.4581163 ])
+
     """
-    estimated_gpp = NEP_TO_GPP_FACTOR * flux_nep
+    assert par.shape == temperature.shape
+    assert flux_npp.shape == flux_rh.shape
+    # It is plausible that NPP < 0 in some seasons; I should probably
+    # use different assumptions then.  GPP = - NPP, Rauto = -2 NPP?
+    estimated_gpp = NPP_TO_GPP_FACTOR * flux_npp
     flux_gpp = olsen_randerson_gpp_once(
         estimated_gpp, par
     )
     flux_resp = olsen_randerson_resp_once(
-        estimated_gpp - flux_nep, temperature
+        # Rauto
+        estimated_gpp - flux_npp +
+        # Rh
+        flux_rh,
+        temperature
     )
     return flux_gpp - flux_resp
 

src/olsen_randerson/fisher.py

@@ -9,7 +9,7 @@
 because that is easy to get pandas to do.
 """
 
-from . import NEP_TO_GPP_FACTOR, Q10, T0
+from . import NPP_TO_GPP_FACTOR, Q10, T0
 
 INPUT_FREQUENCY = "1M"
 """The frequency at which the input data are given.
@@ -19,7 +19,7 @@
 """
 
 
-def downscale_timeseries(flux_nee, temperature, par):
+def downscale_timeseries(flux_npp, flux_rh, temperature, par):
     """Downscale the columns of flux_nee.
 
     The parts of the downscaled flux corresponding to the first and
@@ -29,8 +29,12 @@ def downscale_timeseries(flux_nee, temperature, par):
 
     Parameters
     ----------
-    flux_nee : pd.DataFrame[N_large, M]
-        :abbr:`NEE (Net Ecosystem Exchange)`, at the large timesteps.
+    flux_npp : pd.DataFrame[N_large, M]
+        :abbr:`NPP (Net Primary Production)`, at the large timesteps.
+        Must have datetime index.  Positive indicates carbon is
+        leaving the atmosphere.  Units must have time in denominator.
+    flux_rh : pd.DataFrame[N_large, M]
+        :abbr:`Rh (Heterotrophic Respiration)` at the large timesteps.
         Must have datetime index.  Positive indicates carbon is
         entering the atmosphere.  Units must have time in denominator.
     temperature : pd.DataFrame[N, M]
@@ -47,23 +51,37 @@ def downscale_timeseries(flux_nee, temperature, par):
     flux_nee : pd.DataFrame[N, M]
         The downscaled :abbr:`NEE (Net Ecosystem Exchange)`.
 
+    Notes
+    -----
+    NEE = GPP - Reco = GPP - Ra - Rh = NPP - Rh
+
     References
     ----------
     Fisher, J. B., Sikka, M., Huntzinger, D. N., Schwalm, C., and Liu,
     J., 2016: Technical note: 3-hourly temporal downscaling of monthly
     global terrestrial biosphere model net ecosystem exchange,
     *Biogeosciences*, vol. 13, no. 14, 4271--4277,
     :doi:`10.5194/bg-13-4271-2016`.
+
     """
-    estimated_gpp = -NEP_TO_GPP_FACTOR * flux_nee
+    estimated_gpp = NPP_TO_GPP_FACTOR * flux_npp
     flux_gpp = downscale_gpp_timeseries(
         estimated_gpp, par
     )
     flux_resp = downscale_resp_timeseries(
-        estimated_gpp + flux_nee, temperature
+        estimated_gpp - flux_npp + flux_rh, temperature
     )
     downscaled_nee = flux_resp - flux_gpp
-    return downscaled_nee
+    original_nee = (flux_npp - flux_rh).resample("1MS").first()
+    difference = (
+        downscaled_nee.rolling(
+            "30D", min_periods=1
+        ).mean() -
+        original_nee.resample(par.index.freq).ffill().rolling(
+            "30D", min_periods=1
+        ).mean().loc[par.index[0]:par.index[-1], :]
+    )
+    return downscaled_nee + difference
 
 
 def downscale_gpp_timeseries(flux_gpp, par):

tests/test_fisher.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 import numpy.testing as np_tst
+from numpy.testing import assert_allclose
 import pandas as pd
 
 from hypothesis import given, assume
@@ -113,7 +114,15 @@ def test_downscale_resp_timeseries(flux_resp, temperature):
 @given(
     arrays(
         float, (len(MONTH_CENTER_INDEX), len(COLUMNS)),
-        elements=floats(min_value=-1e30, max_value=1e30)
+        elements=floats(min_value=0, max_value=1e30)
+    ).map(
+        functools.partial(pd.DataFrame,
+                          index=MONTH_CENTER_INDEX,
+                          columns=COLUMNS)
+    ),
+    arrays(
+        float, (len(MONTH_CENTER_INDEX), len(COLUMNS)),
+        elements=floats(min_value=0, max_value=1e30)
     ).map(
         functools.partial(pd.DataFrame,
                           index=MONTH_CENTER_INDEX,
@@ -136,13 +145,17 @@ def test_downscale_resp_timeseries(flux_resp, temperature):
                           columns=COLUMNS)
     )
 )
-def test_downscale_nee_timeseries(flux_nee, temperature, par):
+def test_downscale_timeseries(flux_npp, flux_rh, temperature, par):
     """Test downscaling of NEE."""
     flux_nee_downscaled = olsen_randerson.fisher.downscale_timeseries(
-        flux_nee, temperature, par
+        flux_npp, flux_rh, temperature, par
     )
     assert flux_nee_downscaled.shape == temperature.shape
     flux_nee_downscaled_upscaled = flux_nee_downscaled.resample(
         olsen_randerson.fisher.INPUT_FREQUENCY
     ).sum()
-    assert flux_nee_downscaled_upscaled.shape == flux_nee.shape
+    assert flux_nee_downscaled_upscaled.shape == flux_npp.shape
+    assert_allclose(
+        (flux_npp - flux_rh).values[1:-1, :],
+        flux_nee_downscaled_upscaled
+    )

tests/test_olsen_randerson.py

@@ -15,6 +15,13 @@
 
 Needed to provide bounds for fluxes
 """
+RTOL_FOR_ATOL = 1e-10
+"""rtol to use with NPP/Rh to find atol for NEP/NEE.
+
+Combining NPP and Rh or GPP and Reco to get NEP/NEE subtracts two
+large numbers to get a small number.  Use the original numbers to find
+reasonable precision.
+"""
 
 
 # Not entirely sure what units Photosynthetically Active Radiation is
@@ -78,6 +85,11 @@ def test_olsen_randerson_resp_once(flux_resp, temperature):
         elements=floats(min_value=-UNREASONABLY_LARGE_FLUX_MAGNITUDE,
                         max_value=+UNREASONABLY_LARGE_FLUX_MAGNITUDE)
     ),
+    arrays(
+        np.float, (3, 5),
+        elements=floats(min_value=0,
+                        max_value=+UNREASONABLY_LARGE_FLUX_MAGNITUDE)
+    ),
     arrays(
         np.float, (TEST_LENGTH, 3, 5),
         elements=floats(min_value=-100, max_value=100)
@@ -89,15 +101,19 @@ def test_olsen_randerson_resp_once(flux_resp, temperature):
         lambda par: np.all(np.any(par != 0, axis=0))
     )
 )
-def test_olsen_randerson_once(flux_nee, temperature, par):
+def test_olsen_randerson_once(flux_npp, flux_rh, temperature, par):
     """Test single downscaling of NEE."""
     assume(np.all(np.any(par != 0, axis=0)))
     flux_nee_downscaled = olsen_randerson.olsen_randerson_once(
-        flux_nee, temperature, par
+        flux_npp, flux_rh, temperature, par
     )
     assert flux_nee_downscaled.shape == temperature.shape
     flux_nee_downscaled_upscaled = flux_nee_downscaled.sum(axis=0)
-    assert flux_nee_downscaled_upscaled.shape == flux_nee.shape
+    assert flux_nee_downscaled_upscaled.shape == flux_npp.shape
+    atol = RTOL_FOR_ATOL * max(
+        flux_npp.max(),
+        flux_rh.max(),
+    )
     np_tst.assert_allclose(flux_nee_downscaled_upscaled,
-                           flux_nee * TEST_LENGTH,
-                           atol=1e-100)
+                           (flux_npp - flux_rh) * TEST_LENGTH,
+                           atol=atol)