"""
The module contains various physics calculations needed for estimating
the thermal radition and associated values.
"""
import numpy as np
import os
on_rtd = os.environ.get('READTHEDOCS') == 'True'
if on_rtd:
IPW = '.' # placehold while building the docs
elif 'IPW' not in os.environ:
IPW = '/usr/local/bin'
else:
IPW = os.environ['IPW'] # IPW executables
# define some constants
MAXV = 1.0 # vis albedo when gsize = 0
MAXIR = 0.85447 # IR albedo when gsize = 0
IRFAC = -0.02123 # IR decay factor
VFAC = 500.0 # visible decay factor
VZRG = 1.375e-3 # vis zenith increase range factor
IRZRG = 2.0e-3 # ir zenith increase range factor
IRZ0 = 0.1 # ir zenith increase range, gsize=0
STEF_BOLTZ = 5.6697e-8 # stephman boltzman constant
EMISS_TERRAIN = 0.98 # emissivity of the terrain
EMISS_VEG = 0.96 # emissivity of the vegitation
FREEZE = 273.16 # freezing temp K
BOIL = 373.15 # boiling temperature K
STD_LAPSE_M = -0.0065 # lapse rate (K/m)
STD_LAPSE = -6.5 # lapse rate (K/km)
SEA_LEVEL = 1.013246e5 # sea level pressure
RGAS = 8.31432e3 # gas constant (J / kmole / deg)
GRAVITY = 9.80665 # gravity (m/s^2)
MOL_AIR = 28.9644 # molecular weight of air (kg / kmole)
[docs]def thermal_correct_terrain(th, ta, viewf):
"""
Correct the thermal radiation for terrain assuming that
the terrain is at the air temperature and the pixel and
a sky view
Args:
th: thermal radiation
ta: air temperature [C]
viewf: sky view factor from view_f
Returns:
corrected thermal radiation
20150611 Scott Havens
"""
# thermal emitted from the terrain
terrain = STEF_BOLTZ * EMISS_TERRAIN * np.power(ta + 273.15, 4)
# correct the incoming thermal
return viewf * th + (1 - viewf) * terrain
[docs]def thermal_correct_canopy(th, ta, tau, veg_height, height_thresh=2):
"""
Correct thermal radiation for vegitation. It will only correct
for pixels where the veg height is above a threshold. This ensures
that the open areas don't get this applied. Vegitation temp
is assumed to be at air temperature
Args:
th: thermal radiation
ta: air temperature [C]
tau: transmissivity of the canopy
veg_height: vegitation height for each pixel
height_thresh: threshold hold for height to say that there is veg in
the pixel
Returns:
corrected thermal radiation
Equations from Link and Marks 1999 :cite:`Link&Marks:1999`
20150611 Scott Havens
"""
# thermal emitted from the canopy
veg = STEF_BOLTZ * EMISS_VEG * np.power(ta + FREEZE, 4)
# pixels with canopy above the threshold
ind = veg_height > height_thresh
# correct incoming thermal
th[ind] = tau[ind] * th[ind] + (1 - tau[ind]) * veg[ind]
return th
[docs]def precipitable_water(ta, ea):
"""
Estimate the precipitable water from Prata (1996) :cite:`Prata:1996`
"""
return 4650*ea/ta
[docs]def Dilly1998(ta, ea):
"""
Estimate clear-sky downwelling long wave radiation from Dilley & O'Brian
(1998) :cite:`Dilley&OBrian:1998` using the equation:
.. math::
L_{clear} = 59.38 + 113.7 * \\left( \\frac{T_a}{273.16} \\right)^6 + 96.96 \\sqrt{w/25}
Where :math:`T_a` is the air temperature and :math:`w` is the amount of
precipitable water. The preipitable water is estimated as
:math:`4650 e_o/T_o` from Prata (1996) :cite:`Prata:1996`.
Args:
ta: distributed air temperature [degree C]
ea: distrubted vapor pressure [kPa]
Returns:
clear sky long wave radiation [W/m2]
20170509 Scott Havens
"""
ta = ta + FREEZE # convert to K
w = precipitable_water(ta, ea) # precipitable water
th = 59.38 + 113.7 * (ta/273.16)**6 + 96.96*np.sqrt(w/25)
return th
[docs]def calc_long_wave(e, ta):
"""
Apply the Stephan-Boltzman equation for longwave
"""
return e * STEF_BOLTZ * np.power(ta, 4)
[docs]def Prata1996(ta, ea):
"""
Estimate clear-sky downwelling long wave radiation from Prata (1996)
:cite:`Prata:1996` using the equation:
.. math::
\epsilon_{clear} = 1 - (1 + w) * exp(-1.2 + 3w)^{1/2}
Where :math:`w` is the amount of precipitable water. The preipitable
water is estimated as :math:`4650 e_o/T_o` from Prata (1996)
:cite:`Prata:1996`.
Args:
ta: distributed air temperature [degree C]
ea: distrubted vapor pressure [kPa]
Returns:
clear sky long wave radiation [W/m2]
20170509 Scott Havens
"""
ta = ta + FREEZE # convert to K
w = precipitable_water(ta, ea) # precipitable water
# clear sky emmissivity
ec = 1 - (1 + w) * np.exp(-np.sqrt(1.2 + 3*w))
return calc_long_wave(ec, ta)
[docs]def Angstrom1918(ta, ea):
"""
Estimate clear-sky downwelling long wave radiation from Angstrom (1918)
:cite:`Angstrom:1918` as cited by Niemela et al (2001)
:cite:`Niemela&al:2001` using the equation:
.. math::
\\epsilon_{clear} = 0.83 - 0.18 * 10^{-0.067 e_a}
Where :math:`e_a` is the vapor pressure.
Args:
ta: distributed air temperature [degree C]
ea: distrubted vapor pressure [kPa]
Returns:
clear sky long wave radiation [W/m2]
20170509 Scott Havens
"""
ta = ta + FREEZE # convert to K
e = 0.83 - 0.18 * np.power(10, -0.067*ea)
return calc_long_wave(e, ta)
[docs]def hysat(pb, tb, L, h, g, m):
"""
integral of hydrostatic equation over layer with linear temperature
variation
Args:
pb: base level pressure
tb: base level temp [K]
L: lapse rate [deg/km]
h: layer thickness [km]
g: grav accel [m/s^2]
m: molec wt [kg/kmole]
Returns:
hydrostatic results
20151027 Scott Havens
"""
# the factors 1.e-3 and 1.e3 are for units conversion
if L == 0:
return pb * np.exp(-g * m * h * 1.e3/(RGAS * tb))
else:
return pb * np.power(tb/(tb + L * h), g * m/(RGAS * L * 1.e-3))
[docs]def satw(tk):
"""
Saturation vapor pressure of water. from IPW satw
Args:
tk: temperature in Kelvin
Returns:
saturated vapor pressure over water
20151027 Scott Havens
"""
# remove bad values
tk[tk < 0] = np.nan
l10 = np.log(10.0)
btk = BOIL/tk
x = -7.90298*(btk - 1.0) + 5.02808*np.log(btk)/l10 - \
1.3816e-7*(np.power(10.0, 1.1344e1*(1.0 - tk/BOIL))-1.) + \
8.1328e-3*(np.power(10.0, -3.49149*(btk - 1.0)) - 1.0) + \
np.log(SEA_LEVEL)/l10
x = np.power(10.0, x)
return x
[docs]def sati(tk):
"""
saturation vapor pressure over ice. From IPW sati
Args:
tk: temperature in Kelvin
Returns:
saturated vapor pressure over ice
20151027 Scott Havens
"""
# remove bad values
tk[tk < 0] = np.nan
# preallocate
x = np.empty(tk.shape)
# vapor above freezing
ind = tk > FREEZE
x[ind] = satw(tk[ind])
# vapor below freezing
l10 = np.log(10.0)
x[~ind] = 100.0 * np.power(10.0, -9.09718*((FREEZE/tk[~ind]) - 1.0) -
3.56654*np.log(FREEZE/tk[~ind])/l10 +
8.76793e-1*(1.0 - (tk[~ind]/FREEZE)) +
np.log(6.1071)/l10)
return x
[docs]def brutsaert(ta, l, ea, z, pa):
"""
Calculate atmosphere emissivity from Brutsaert (1975):cite:`Brutsaert:1975`
Args:
ta: air temp (K)
l: temperature lapse rate (deg/m)
ea: vapor pressure (Pa)
z: elevation (z)
pa: air pressure (Pa)
Returns:
atmosphericy emissivity
20151027 Scott Havens
"""
t_prime = ta - (l * z)
rh = ea / sati(ta)
rh[rh > 1] = 1
e_prime = (rh * sati(t_prime))/100.0
air_emiss = (1.24*np.power(e_prime/t_prime, 1./7.0))*pa/SEA_LEVEL
air_emiss[air_emiss > 1.0] = 1.0
return air_emiss
[docs]def topotherm(ta, tw, z, skvfac):
"""
Calculate the clear sky thermal radiation. topotherm calculates thermal
radiation from the atmosphere corrected for topographic effects, from near
surface air temperature Ta, dew point temperature DPT, and elevation. Based
on a model by Marks and Dozier (1979) :citeL`Marks&Dozier:1979`.
Args:
ta: air temperature [C]
tw: dew point temperature [C]
z: elevation [m]
skvfac: sky view factor
Returns:
Long wave (thermal) radiation corrected for terrain
20151027 Scott Havens
"""
# convert ta and tw from C to K
ta = ta + FREEZE
tw = tw + FREEZE
# if below zero set to nan
tw[tw < 0] = np.nan
ta[ta < 0] = np.nan
# calculate theoretical sea level
# atmospheric emissivity
# from reference level ta, tw, and z
ind = tw > ta
tw[ind] = ta[ind]
ea = sati(tw)
emiss = brutsaert(ta, STD_LAPSE_M, ea, z, SEA_LEVEL)
# calculate sea level air temp
T0 = ta - (z * STD_LAPSE_M)
# adjust emiss for elev, terrain
# veg, and cloud shading
press = hysat(SEA_LEVEL, T0, STD_LAPSE, z/1000.0, GRAVITY, MOL_AIR)
# elevation correction
emiss *= press/SEA_LEVEL
# terrain factor correction
emiss = (emiss * skvfac) + (1.0 - skvfac)
# check for emissivity > 1.0
emiss[emiss > 1.0] = 1.0
# calculate incoming lw rad
return emiss * STEF_BOLTZ * np.power(ta, 4)
[docs]def Garen2005(th, cloud_factor):
"""
Cloud correction is based on the relationship in Garen and Marks (2005)
:cite:`Garen&Marks:2005` between the cloud factor and measured long
wave radiation using measurement stations in the Boise River Basin.
.. math::
L_{cloud} = L_{clear} * (1.485 - 0.488 * cloud\_factor)
Args:
th: clear sky thermal radiation [W/m2]
cloud_factor: fraction of sky that are not clouds, 1 equals no clouds,
0 equals all clouds
Returns:
cloud corrected clear sky thermal radiation
20170515 Scott Havens
"""
return th * (1.485 - 0.488 * cloud_factor)
[docs]def Unsworth1975(th, ta, cloud_factor):
"""
Cloud correction is based on Unsworth and Monteith (1975)
:cite:`Unsworth&Monteith:1975`
.. math::
\epsilon_a = (1 - 0.84) \epsilon_{clear} + 0.84c
where :math:`c = 1 - cloud\_factor`
Args:
th: clear sky thermal radiation [W/m2]
ta: temperature in Celcius that the clear sky thermal radiation was
calcualted from [C] cloud_factor: fraction of sky that are not
clouds, 1 equals no clouds, 0 equals all clouds
Returns:
cloud corrected clear sky thermal radiation
20170515 Scott Havens
"""
c = 1 - cloud_factor
ta = ta + FREEZE
# calculate the clear sky emissivity
ec = th/(STEF_BOLTZ * np.power(ta, 4))
# adjsut the atmospheric emissivity for clouds
ea = (1 - 0.84*c) * ec + 0.84*c
return calc_long_wave(ea, ta)
[docs]def Kimball1982(th, ta, ea, cloud_factor):
"""
Cloud correction is based on Kimball et al. (1982) :cite:`Kimball&al:1982`
.. math::
L_d &= L_{clear} + \\tau_8 c f_8 \sigma T^{4}_{c}
\\tau_8 &= 1 - \epsilon_{8z} (1.4 - 0.4 \epsilon_{8z})
\epsilon_{8z} &= 0.24 + 2.98 \\times 10^{-6} e^2_o exp(3000/T_o)
f_8 &= -0.6732 + 0.6240 \\times 10^{-2} T_c - 0.9140 \\times 10^{-5} T^2_c
where the original Kimball et al. (1982) :cite:`Kimball&al:1982` was for
multiple cloud layers, which was simplified to one layer. :math:`T_c` is
the cloud temperature and is assumed to be 11 K cooler than :math:`T_a`.
Args:
th: clear sky thermal radiation [W/m2]
ta: temperature in Celcius that the clear sky thermal radiation was
calcualted from [C]
ea: distrubted vapor pressure [kPa]
cloud_factor: fraction of sky that are not clouds, 1 equals no clouds,
0 equals all clouds
Returns:
cloud corrected clear sky thermal radiation
20170515 Scott Havens
"""
ta = ta + FREEZE
Tc = ta - 11
f8 = -0.6732 + 0.6240 * 10**(-2) * Tc - 0.9140 * 10**(-5) * Tc**2
e8z = 0.24 + 2.98 * 10**(-6) * ea**2 * np.exp(3000/ta)
t8 = 1 - e8z * (1.4 - 0.4 * e8z)
return th + calc_long_wave(f8 * e8z * t8, ta)
[docs]def Crawford1999(th, ta, cloud_factor):
"""
Cloud correction is based on Crawford and Duchon (1999)
:cite:`Crawford&Duchon:1999`
.. math::
\epsilon_a = (1 - cloud\_factor) + cloud\_factor * \epsilon_{clear}
where :math:`cloud\_factor` is the ratio of measured solar radiation to the
clear sky irradiance.
Args:
th: clear sky thermal radiation [W/m2]
ta: temperature in Celcius that the clear sky thermal radiation was
calcualted from [C]
cloud_factor: fraction of sky that are not clouds, 1 equals no clouds,
0 equals all clouds
Returns:
cloud corrected clear sky thermal radiation
20170515 Scott Havens
"""
return calc_long_wave(1 - cloud_factor, ta + FREEZE) + cloud_factor * th