"""
The Sap flow toolbox
====================
Sap velocity and sap flow is a very interesting means to monitor xylem water
dynamics in higher plants (i.e. trees in our case). Since a mere flow conversion
is likely to result in erroneous assumptions (Čermák et al., 2004), we have
compiled a couple of sapwood-related functions to estimate active sapwood area
for sap velocity conversion to sap flow.
.. note::
This toolbox is by no means complete nor exhaustive. Please regard it as
helper functions which require throughout testing and deserve substantial
extension to further application cases.
.. note::
To get started and for direct application (tested for beech trees) use rootwater.sapflow.sap_calc
and provide a pandas.DataFrame with measured sap velocity from East30 sensors
(in cm/h).
References
----------
Čermák, J., J. Kučera, and N. Nadezhdina (2004), Sap flow measurements with some
thermodynamic methods, flow integration within trees and scaling up from sample
trees to entire forest stands, Trees, 18(5), 529–546, doi:10.1007/s00468-004-0339-6.
"""
import json
import sys
import os
import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar
sys.path.append(os.path.dirname(__file__))
from gebauer_params import gp
#gp : dictionary for all valid tree names. Each name
# has to be key to a nested dict that defines the
# four Weibull parameters a,b,c,d
[docs]def roessler(r, tree='beech'):
r"""Estimate bark thickness
Returns the estimated bark thickness as published by Rössler (2008)
Parameters
----------
r : float
tree radius at breast height (in cm)
tree : str
Tree name, for which to calculate bark thickness.
Can be one of ['beech', 'oak']
Returns
-------
db : float
bark thickness (in mm)
Raises
------
NotImplementedError : if tree is not in ('oak', 'beech')
References
----------
Rössler, G.: Rindenabzug richtig bemessen, Forstzeitung, 4, p. 21, 2008.
"""
if tree=='beech':
db = 2.61029 + 0.28522 * 2 * r
elif tree=='oak':
db = 9.88855 + 0.56734 * 2 * r
else:
raise NotImplementedError("Only 'beech' and 'oak' supported.")
return db / 10
[docs]def gebauer(r, tree='beech'):
r"""Sap-wood thickness
Calculates sap-wood thickness as published by Gebauer et al. (2008)
Parameters
----------
r : float
tree radius at breast height (in cm)
tree : str
Tree name, for which to calculate bark and sapwood thickness.
Can be one of ['beech', 'oak']
Returns
-------
th : float
sap-wood thickness (in mm)
Raises
------
NotImplementedError : if tree is not in ('oak', 'beech')
References
----------
Gebauer, T., Horna, V., and Leuschner, C.: Variability in radial sap flux
density patterns and sapwood area among seven co-occurring temperate
broad-leaved tree species, Tree Physiol., 28, 1821–1830, 2008.
"""
r = r - roessler(r, tree=tree) / 2.
if tree=='beech':
As = 0.778 * (2*r)**1.917
elif tree=='oak':
As = 0.065 * (2*r)**2.264
else:
raise NotImplementedError("Only 'beech' and 'oak' supported.")
return -1.*(np.sqrt((np.pi*r**2 - As)/np.pi)-r)
[docs]def gebauer_weibull(x,a,b,c,d):
r"""4-parameter Weibull function after Gebauer
Calls Weibull distribution function as published by
Gebauer et al. (2008).
Parameters
----------
x : float or numpy.ndarray
realtive sampling points for distribution function
a, b, c, d : float
Weibull function parameters, which are tree-specific in this case
References
----------
Gebauer, T., Horna, V., and Leuschner, C.: Variability in radial sap flux
density patterns and sapwood area among seven co-occurring temperate
broad-leaved tree species, Tree Physiol., 28, 1821–1830, 2008.
"""
#Weibull function with 4 tree-specific parameters abcd
return (c-1)/c + (a*((c-1)/c)**((c-1)/c))*np.exp(-1.*((x-d)/b + (c-1/c)**(1/c))**c) * ((x-d)/b + (c-1/c)**(1/c))**(c-1)
[docs]def get_default_gp():
r"""read default gp
Loads default Weibull distribution parameters as published by
Gebauer et al. (2008). Can be used by
:func:`gebauer_weibull <rootwater.sf.gebauer_weibull>` to
calculate sap velocity distribution in sapwood
References
----------
Gebauer, T., Horna, V., and Leuschner, C.: Variability in radial sap flux
density patterns and sapwood area among seven co-occurring temperate
broad-leaved tree species, Tree Physiol., 28, 1821–1830, 2008.
"""
with open(GP_PATH, 'r') as fs:
params = json.load(fs)
# if needed, you could validate the params here
return params
[docs]def gebauer_rel(r, tree='beech', n_points=50):
r"""relative flux density
Calculates relative flux density as a function
of depth on sapwood for n_points.
Parameters
----------
r : float
tree radius at breast height (in cm)
tree : str
Tree name, for which to calculate Weibull function.
Tree name has to be in gp.keys()
n_points : int
Number of points for solving Weibull.
This is the resolution over depth.
Returns
-------
sv : numpy.ndarray
relative flux density at n_points
References
----------
Gebauer, T., Horna, V., and Leuschner, C.: Variability in radial sap flux
density patterns and sapwood area among seven co-occurring temperate
broad-leaved tree species, Tree Physiol., 28, 1821–1830, 2008.
"""
x=np.arange(n_points)/ n_points *gebauer(r)
p = gp.get(tree)
if p is None:
raise ValueError('Tree %s is unknown' % tree)
return gebauer_weibull(x, p['a'], p['b'], p['c'], p['d'])
[docs]def recko(r,hydra=False):
r"""sap-wood thickness after Račko et al. (2018)
Calculates sap-wood thickness as a function of the tree radius
at breast height.
Parameters
----------
r : float
tree radius at breast height (in cm)
hydra : bool
selects if only the hydrated area is returned (when True)
Returns
-------
sapwood : float
sapwood thickness (in cm)
References
----------
Račko, V., O. Mišíková, P. Hlaváč, and V. Deáková (2018),
Can bark stripping cause red heartwood formation in beech stems?
iForest - Biogeosciences and Forestry, 11(2), 251–258, doi:10.3832/ifor2147-011.
"""
if hydra:
#only hydrated area
return 0.34*r - 2.714
else:
return 0.3748*r
[docs]def galvac(r):
r"""sap-wood thickness after Galvac et al. (1990)
Calculates sap-wood thickness as a function of the tree radius
at breast height.
Parameters
----------
r : float
tree radius at breast height (in cm)
Returns
-------
sapwood : float
sapwood thickness (in cm)
References
----------
Glavac, V., Koenies, H. & Ebben, U. Holz als Roh- und Werkstoff (1990) 48: 437. https://doi.org/10.1007/BF02627628
"""
#returns sap-wood thickness
dbh = r*2
As = 0.6546*dbh**2 + 0.5736*dbh - 40.069
return -1.*(np.sqrt((np.pi*r**2 - As)/np.pi)-r)
[docs]def gebauer_act(r,perc=0.95,tree='beech'):
r"""Active sapwood area based on percentile of Weibull distribution
Calculates the "zero" sap velocity limit as given percentile of relative
flux velocity distribution as a Weibull function after Gebauer et al. (2008)
Parameters
----------
r : float or numpy.ndarray
tree radius at breast height (in cm)
perc : float
percentile to define the "zero" sap velocity limit
tree : str
Tree name, for which to calculate Weibull function.
Tree name has to be in gp.keys()
Returns
-------
act_sap : float or numpy.ndarray
depth of "zero" sap velocity limit as active sapwood in tree (in cm)
References
----------
Gebauer, T., Horna, V., and Leuschner, C.: Variability in radial sap flux
density patterns and sapwood area among seven co-occurring temperate
broad-leaved tree species, Tree Physiol., 28, 1821–1830, 2008.
"""
if type(r)==float:
return np.where(np.cumsum(gebauer_rel(r,tree))/sum(gebauer_rel(r,tree))>perc)[0][0]/50.*gebauer(r,tree)
else:
act_sap = r*np.nan
for i in np.arange(len(r))[1:]:
act_sap[i] = np.where(np.cumsum(gebauer_rel(r[i],tree))/sum(gebauer_rel(r[i],tree))>perc)[0][0]/50.*gebauer(r[i],tree)
return act_sap
[docs]def sap_volume(r,s1,s2,vout=False,perc=0.95,tree='beech'):
r"""Estimate sap flow from sap velocity in inner sapwood measured with East30 sensors
Calculates the sap flow after Gebauer et al. (2008) based on sap velocity measurements
by fitting of Gebauer-Weibull function to measured sap velocity at mid and inner point
through a scaling factor (but not changing the empirical, tree-specific parameters).
Parameters
----------
r : float
tree radius at breast height (in cm)
s1 : float or pandas.Series (datetime index is preferable)
sap velocity measurement at mid point of East30 sensor (cm/time)
s2 : float or pandas.Series (datetime index is preferable)
sap velocity measurement at inner point of East30 sensor (cm/time)
vout : bool
True returns aggregated volume flux (cm3/time), False returns velocity distribution (cm/time)
perc : float
percentile to define the "zero" sap velocity limit
tree : str
Tree name, for which to calculate bark thickness and Weibull function.
Tree name has to be in gp.keys()
Returns
-------
return : float or pandas.Series
aggregated volume flux (cm3/time) (if vout is True), or
returns velocity distribution (cm/time) (if vout is False)
References
----------
Gebauer, T., Horna, V., and Leuschner, C.: Variability in radial sap flux
density patterns and sapwood area among seven co-occurring temperate
broad-leaved tree species, Tree Physiol., 28, 1821–1830, 2008.
"""
xi = np.arange(50)/50.*gebauer(r,tree)
def aply_geb(s1x):
dummy = s1x*gebauer_rel(r,tree)
er1 = (dummy[np.arange(50)/50.*gebauer(r,tree) >= 1.8][0]-s1)**2
er2 = (dummy[np.arange(50)/50.*gebauer(r,tree) >= 3.][0]-s2)**2
return np.sqrt(0.2*er1+er2)
res = minimize_scalar(aply_geb)
dummy = res.x*gebauer_rel(r,tree)
v3 = dummy[(np.arange(50)/50.*gebauer(r,tree) > 2.4) & (np.arange(50)/50.*gebauer(r,tree) <= gebauer_act(r,perc,tree))]
rx = (np.arange(50)/50.*gebauer(r,tree))[(np.arange(50)/50.*gebauer(r,tree) > 2.4) & (np.arange(50)/50.*gebauer(r,tree) <= gebauer_act(r,perc,tree))]
Ax = rx*np.nan
for i in np.arange(len(Ax)):
Ax = A_circ(r,[rx[i]-0.01*gebauer(r,tree),rx[i]+0.01*gebauer(r,tree)])
if vout:
return dummy
else:
return np.sum(Ax * v3)
[docs]def A_circ(r,sens=[0.,1.1],tree='beech'):
r"""Calculate area of circular ring
Simple geometrical calculation of a circular ring area as reference cross-section
for sap flow calculation.
Parameters
----------
r : float
tree radius at breast height (in cm)
sens : list of floats
outer and inner point of ring (in cm)
tree : str
Tree name, for which to calculate and subtract bark thickness.
Can be one of ['beech', 'oak'].
Returns
-------
return : float
area of ring as
"""
#East30 thermocouplers location at needles 5, 18 and 30 mm
r = r - roessler(r, tree)/2. #correct for bark
return np.pi*((r-sens[0])**2) - np.pi*((r-sens[1])**2)
[docs]def sap_calc(SV,r,perc=0.95,tree='beech'):
r"""Wrapper for sap flow calculation with rootwater.sapflow.sap_volume
Calculates the sap flow after Gebauer et al. (2008) based on measured sap velocity
with East30 sensors.
Parameters
----------
SV : pandas.DataFrame
sap velocity (in cm/h) in three columns ordered inner, mid, outer point
r : float
tree radius at breast height (in cm)
perc : float
percentile to define the "zero" sap velocity limit
tree : str
Tree name, for which to calculate bark thickness and Weibull function.
Tree name has to be in gp.keys()
Returns
-------
return : pandas.DataFrame
sap volume flux (cm3/h)
"""
Sap = SV.copy()*np.nan
colx = SV.columns[:3]
for i in Sap.index:
Sap.loc[i,colx[0]] = sap_volume(r,SV.loc[i,colx[1]],SV.loc[i,colx[0]],False,perc,tree)
Sap.loc[i,colx[1]] = SV.loc[i,colx[1]]*A_circ(r,[1.1,2.4],tree)
Sap.loc[i,colx[2]] = SV.loc[i,colx[2]]*A_circ(r,[0.,1.1],tree)
return Sap
def stackplot(A):
import matplotlib.pyplot as plt
r"""plot stacked time series (of first three columns of the provided dataframe)
Parameters
----------
A : pandas.DataFrame (preferrably with datetime index)
DataFrame with three columns to be stacked
Returns
-------
plot
"""
#plot stacked time series of first three columns of the dataframe A
# These are the "Tableau 20" colors as RGB.
tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),
(44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),
(148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),
(227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),
(188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]
# Scale the RGB values to the [0, 1] range, which is the format matplotlib accepts.
for i in range(len(tableau20)):
r, g, b = tableau20[i]
tableau20[i] = (r / 255., g / 255., b / 255.)
tableau10=tableau20[0::2]
plt.fill_between(A.index,A.iloc[:,2],facecolor=tableau10[0],alpha=0.7,color='b',lw=0,label=A.columns[2])
plt.fill_between(A.index,A.iloc[:,[1,2]].sum(axis=1),A.iloc[:,2],facecolor=tableau10[2],alpha=0.7,color='g',lw=0,label=A.columns[1])
plt.fill_between(A.index,A.iloc[:,:3].sum(axis=1),A.iloc[:,[1,2]].sum(axis=1),facecolor=tableau10[3],alpha=0.7,color='y',lw=0,label=A.columns[0])