import numpy as np
from scipy.special import gamma
from pycaps.util.util import abstract
[docs]class PSD(object):
"""
Represents one or more particle size distributions for any hydrometeor.
Args:
q_hydro (np.array): Hydrometeor mixing ratio (kg/kg)
density_hydro (np.array): Hydrometeor density (kg/m^3)
alpha_hydro (np.array): Shape parameter for the PSD
n0_hydro (np.array): Intercept parameter for the PSD
nt_hydro (np.array): Total number concentration for the PSD
"""
_dielectric_water = 0.93
_dielectric_ice = 0.176
def __init__(self, q_hydro, density_hydro, alpha_hydro=0., n0_hydro=None, nt_hydro=None):
"""
Constructor for the PSD class
"""
self._q_hydro = q_hydro
self._density_hydro = density_hydro
self._n0_hydro = n0_hydro
self._nt_hydro = nt_hydro
if nt_hydro is not None and n0_hydro is not None:
print "PSD object needs either an intercept parameter (n0_hydro) or total number concentration (nt_hydro)"
import sys; sys.exit()
self._alpha_hydro = alpha_hydro
return
@property
def density(self):
"""
Hydrometeor density.
"""
return self._density_hydro
@property
def mixingrat(self):
"""
Hydrometeor mixing ratio.
"""
return self._q_hydro
[docs] def interceptParam(self, density_air):
"""
Get the intercept parameter for the distribution.
Args:
density_air (np.array): Density of the air (kg/m^3).
Returns:
The intercept parameter for the PSD.
"""
if self._n0_hydro is None:
self._n0_hydro = self._computeN0(density_air)
return self._n0_hydro
def _computeN0(self, density_air):
"""
Compute the intercept parameter given air density and total number concentration.
Args:
density_air (np.array): Density of the air (kg/m^3)
Returns:
The intercept parameter for the PSD. Units are 1/m^4.
"""
gamma1 = gamma(1. + self._alpha_hydro)
gamma4 = gamma(4. + self._alpha_hydro)
lamda = np.where((self._q_hydro > 0.),
((gamma4 / gamma1) * (np.pi / 6. * self._density_hydro) * self._nt_hydro / (density_air * self._q_hydro)) ** (1. / 3.),
0
)
N0 = self._nt_hydro / gamma1 * lamda ** (1. + self._alpha_hydro)
# print "Max intercept parameter:", np.nanmax(N0)
return N0
@abstract
[docs] def reflectivityZhang(self):
"""
Compute this PSD's contribution to linear reflectivity using the Zhang method.
"""
return
@abstract
[docs] def reflectivityJZX(self):
"""
Compute this PSD's contribution to linear reflectivity using the method of Jung, Zhang, and Xue (2008).
"""
return
@staticmethod
def _waterFractionJZX(water, frozen, max_mixed_frac):
"""
Compute the water fraction in mixed phase regions for the Jung, Zhang, and Xue (2008) method.
Args:
water (PSD): The liquid PSD to use for the water fraction calculation.
frozen (PSD): The frozen PSD to use for the water fraction calculation.
max_mixed_frac (float): The maximum fraction of mass of either PSD that can be of mixed phase.
Returns:
A tuple containing the water fraction, mass of water tied up in mixed phase, mass of ice tied up in mixed
phase, total mass of mixed phase, and density of that mixture.
"""
np_err_settings = np.geterr()
np.seterr(all='ignore')
mixed_frac = np.where((water.mixingrat > 0) & (frozen.mixingrat > 0),
max_mixed_frac * np.minimum(frozen.mixingrat / water.mixingrat, water.mixingrat / frozen.mixingrat) ** 0.3,
0
)
water_melt_mass = mixed_frac * water.mixingrat
ice_melt_mass = mixed_frac * frozen.mixingrat
total_melt_mass = water_melt_mass + ice_melt_mass
water_frac = np.where(total_melt_mass > 0,
water_melt_mass / total_melt_mass,
0
)
density_mix = water.density * water_frac ** 2. + (1. - water_frac ** 2.) * frozen.density
np.seterr(**np_err_settings)
return water_frac, water_melt_mass, ice_melt_mass, total_melt_mass, density_mix
@staticmethod
def _partialHReflectivityFrozenJZX(T_mat_coeff, T_mat_alpha, density_air, mixingrat_spec, density_spec, N0_spec):
"""
Compute the horizontal reflectivity contribution for a frozen water species using the method of Jung, Zhang, and Xue (2008).
Args:
T_mat_coeff (tuple): T-Matrix coefficients
T_mat_alpha (tuple): T-Matrix alphas
density_air (np.array): Density of the air (kg/m^3).
mixingrat_spec (np.array): Mixing ratio for the hydrometeor species (kg/kg)
density_spec (np.array): Density of the hydrometeor species (kg/m^3)
N0_spec (np.array): Intercept parameter for the hydrometeor species
Returns:
Linear reflectivity contribution for an ice species
"""
coef = 5.6212976e26 # Note: this may assume a 10-cm wavelength radar
A, B, C = T_mat_coeff
alpha_a, alpha_b = T_mat_alpha
Z_spec = np.where(N0_spec > 0,
coef * (A * alpha_a ** 2 + B * alpha_b ** 2 + 2 * C * alpha_a * alpha_b) * (density_air * mixingrat_spec / density_spec) ** 1.75 / N0_spec ** 0.75,
0
)
# print "Max linear reflectivity for this species:", np.nanmax(Z_spec)
return Z_spec
[docs]class RainPSD(PSD):
"""
Represents one or more particle size distributions for rain.
Args:
qrain (np.array): Rain mixing ratio (kg/kg)
density (np.array): Hydrometeor density (kg/m^3). Optional, default is 1000 kg/m^3.
shape_param (np.array): Shape parameter for the PSD. Optional, default is 0.
intercept_param (np.array): Intercept parameter for the PSD in 1/m^4.
num_concentration (np.array): Total number concentration for the PSD in #/m^3.
"""
def __init__(self, qrain, density=1000., shape_param=0., intercept_param=None, num_concentration=None):
"""
Constructor for the RainPSD class
"""
super(RainPSD, self).__init__(qrain, density, shape_param, intercept_param, num_concentration)
return
[docs] def reflectivityZhang(self, density_air, is_wet):
"""
Compute the rain PSD's contribution to linear reflectivity using the Zhang method.
Args:
density_air (np.array): Density of the air (kg/m^3)
is_wet (np.array): Array of integers (either 0 or 1) telling whether or not this region is above 0 C.
Returns:
Rain contribution to linear reflectivity
"""
N0 = self.interceptParam(density_air)
Zr_factor = 720 * 1e18 / ((np.pi * self.density) ** 1.75 * N0 ** 0.75)
Zr = Zr_factor * (density_air * self.mixingrat) ** 1.75
return Zr
[docs] def reflectivityJZX(self, density_air, *frozen_psds):
"""
Compute the rain PSD's contribution to linear reflectivity using the method of Jung, Zhang, and Xue (2008).
Args:
density_air (np.array): Density of the air (kg/m^3)
frozen_psds (list): List of tuples. Each tuple contains a frozen PSD to compute water fractions with and
the maximum mixed fraction. May be called by reflectivityJZX(... *list_of_psds_and_mixed_fractions).
Returns:
Rain contribution to linear reflectivity
"""
# print
# print "Computing reflectivity contribution from rain"
water_mix_masses = []
for fpsd, mmf in frozen_psds:
water_frac, water_mix_mass, frozen_mix_mass, total_mix_mass, density_mix = PSD._waterFractionJZX(self, fpsd, mmf)
water_mix_masses.append(water_mix_mass)
qr_pure = np.maximum(self.mixingrat - np.sum(water_mix_masses, axis=0), 0)
mm3todBZ = 1e9
wavelen = 107.
alphaa = 4.28e-4
gamma7_08 = 836.7818
N0r = self.interceptParam(density_air)
Zr_const = np.where(N0r > 0,
mm3todBZ * (4 * wavelen ** 4 * alphaa ** 2 * gamma7_08) / (np.pi ** 5.77 * PSD._dielectric_water * self.density ** 1.77 * (N0r * 1e-12) ** 0.77),
0
)
Zr = Zr_const * (qr_pure * density_air) ** 1.77
# print "Max linear reflectivity for rain:", np.nanmax(Zr)
return Zr
[docs]class SnowPSD(PSD):
"""
Represents one or more particle size distributions for snow.
Args:
qsnow (np.array): Snow mixing ratio (kg/kg)
density (np.array): Hydrometeor density (kg/m^3). Optional, default is 100 kg/m^3.
shape_param (np.array): Shape parameter for the PSD. Optional, default is 0.
intercept_param (np.array): Intercept parameter for the PSD in 1/m^4.
num_concentration (np.array): Total number concentration for the PSD in #/m^3.
"""
def __init__(self, qsnow, density=100., shape_param=0., intercept_param=None, num_concentration=None):
"""
Constructor for the SnowPSD class
"""
super(SnowPSD, self).__init__(qsnow, density, shape_param, intercept_param, num_concentration)
return
[docs] def reflectivityZhang(self, density_air, is_wet):
"""
Compute the snow PSD's contribution to linear reflectivity using the Zhang method.
Args:
density_air (np.array): Density of the air (kg/m^3)
is_wet (np.array): Array of integers (either 0 or 1) telling whether or not this region is above 0 C.
Returns:
Snow contribution to linear reflectivity
"""
density_ice = 917.
N0 = self.interceptParam(density_air)
Zs_wet_factor = 720 * 1e18 / ((np.pi * self.density) ** 1.75 * N0 ** 0.75)
Zs_dry_factor = 720 * 1e18 * PSD._dielectric_ice * self.density ** 0.25 / (np.pi ** 1.75 * PSD._dielectric_water * N0 ** 0.75 * density_ice ** 2)
Zs = Zs_wet_factor * is_wet * (density_air * self.mixingrat) ** 1.75 + Zs_dry_factor * (1 - is_wet) * (density_air * self.mixingrat) ** 1.75
return Zs
[docs] def reflectivityJZX(self, density_air, rain_psd, max_mixed_frac):
"""
Compute the snow PSD's contribution to linear reflectivity using the method of Jung, Zhang, and Xue (2008).
Args:
density_air (np.array): Density of the air (kg/m^3)
rain_psd (PSD): Rain PSD to compute water fractions with.
max_mixed_frac (float): Maximum mixed fraction for the snow PSD.
Returns:
Snow contribution to linear reflectivity
"""
# print
# print "Computing reflectivity contribution from snow"
water_frac, water_mix_mass, snow_mix_mass, total_mix_mass, density_mix = PSD._waterFractionJZX(rain_psd, self, max_mixed_frac)
qs_wet = snow_mix_mass
qs_dry = np.maximum(self.mixingrat - qs_wet, 0)
A = 0.8140
B = 0.0303
C = 0.0778
D = 0.4221
Ck = 0.7837
alpha_a = (0.194 + 7.094 * water_frac + 2.135 * water_frac ** 2 - 5.225 * water_frac ** 3) * 1e-4
alpha_b = (0.191 + 6.916 * water_frac - 2.841 * water_frac ** 2 - 1.160 * water_frac ** 3) * 1e-4
Zs_dry = PSD._partialHReflectivityFrozenJZX((A, B, C), (alpha_a, alpha_b), density_air, qs_dry, self.density, self.interceptParam(density_air))
Zs_wet = PSD._partialHReflectivityFrozenJZX((A, B, C), (alpha_a, alpha_b), density_air, qs_wet, density_mix, self.interceptParam(density_air))
return Zs_dry + Zs_wet
[docs]class HailPSD(PSD):
"""
Represents one or more particle size distributions for hail.
Args:
qhail (np.array): Hail mixing ratio (kg/kg)
density (np.array): Hydrometeor density (kg/m^3). Optional, default is 913 kg/m^3.
shape_param (np.array): Shape parameter for the PSD. Optional, default is 0.
intercept_param (np.array): Intercept parameter for the PSD in 1/m^4.
num_concentration (np.array): Total number concentration for the PSD in #/m^3.
"""
def __init__(self, qhail, density=913., shape_param=0., intercept_param=None, num_concentration=None):
"""
Constructor for the HailPSD class
"""
super(HailPSD, self).__init__(qhail, density, shape_param, intercept_param, num_concentration)
return
[docs] def reflectivityZhang(self, density_air, is_wet):
"""
Compute the hail PSD's contribution to linear reflectivity using the Zhang method.
Args:
density_air (np.array): Density of the air (kg/m^3)
is_wet (np.array): Array of integers (either 0 or 1) telling whether or not this region is above 0 C.
Returns:
Hail contribution to linear reflectivity
"""
N0 = self.interceptParam(density_air)
Zh_factor = 720 * 1e18 / ((np.pi * self.density) ** 1.75 * N0 ** 0.75)
Zh = (Zh_factor * (density_air * self.mixingrat) ** 1.75) ** 0.95
return Zh
[docs] def reflectivityJZX(self, density_air, rain_psd, max_mixed_frac):
"""
Compute the hail PSD's contribution to linear reflectivity using the method of Jung, Zhang, and Xue (2008).
Args:
density_air (np.array): Density of the air (kg/m^3)
rain_psd (PSD): Rain PSD to compute water fractions with.
max_mixed_frac (float): Maximum mixed fraction for the hail PSD.
Returns:
Hail contribution to linear reflectivity
"""
# print
# print "Computing reflectivity contribution from hail"
water_frac, water_mix_mass, hail_mix_mass, total_mix_mass, density_mix = PSD._waterFractionJZX(rain_psd, self, max_mixed_frac)
qh_wet = hail_mix_mass
qh_dry = np.maximum(self.mixingrat - qh_wet, 0)
sigma = np.pi / 3 * (1 - 0.8 * water_frac)
A = .125 * (3 + 4 * np.exp(-2 * sigma ** 2) + np.exp(-8 * sigma ** 2))
B = .125 * (3 - 4 * np.exp(-2 * sigma ** 2) + np.exp(-8 * sigma ** 2))
C = .125 * (1 - np.exp(-8 * sigma ** 2))
D = .125 * (3 + np.exp(-8 * sigma ** 2))
Ck = np.exp(-2 * sigma ** 2)
alpha_a = (0.191 + 2.39 * water_frac - 12.57 * water_frac ** 2 + 38.71 * water_frac ** 3 - 65.53 * water_frac ** 4 + 56.16 * water_frac ** 5 - 18.98 * water_frac ** 6) * 1e-3
alpha_b = (0.165 + 1.72 * water_frac - 9.92 * water_frac ** 2 + 32.15 * water_frac ** 3 - 56.0 * water_frac ** 4 + 48.83 * water_frac ** 5 - 16.69 * water_frac ** 6) * 1e-3
Zh_dry = PSD._partialHReflectivityFrozenJZX((A, B, C), (alpha_a, alpha_b), density_air, qh_dry, self.density, self.interceptParam(density_air))
Zh_wet = PSD._partialHReflectivityFrozenJZX((A, B, C), (alpha_a, alpha_b), density_air, qh_wet, density_mix, self.interceptParam(density_air))
return Zh_dry + Zh_wet
[docs]class GraupelPSD(PSD):
"""
Represents one or more particle size distributions for graupel.
Args:
qgraupel (np.array): Graupel mixing ratio (kg/kg)
density (np.array): Hydrometeor density (kg/m^3). Optional, default is 400 kg/m^3.
shape_param (np.array): Shape parameter for the PSD. Optional, default is 0.
intercept_param (np.array): Intercept parameter for the PSD in 1/m^4.
num_concentration (np.array): Total number concentration for the PSD in #/m^3.
"""
def __init__(self, qgraupel, density=400., shape_param=0., intercept_param=None, num_concentration=None):
"""
Constructor for the GraupelPSD class
"""
super(GraupelPSD, self).__init__(qgraupel, density, shape_param, intercept_param, num_concentration)
return
[docs] def reflectivityZhang(self, density_air, iswet):
"""
Compute the graupel PSD's contribution to linear reflectivity using the Zhang method.
Args:
density_air (np.array): Density of the air (kg/m^3)
is_wet (np.array): Array of integers (either 0 or 1) telling whether or not this region is above 0 C.
Returns:
Graupel contribution to linear reflectivity. Will be 0 for the Zhang method.
"""
return 0.
[docs] def reflectivityJZX(self, density_air, rain_psd, max_mixed_frac):
"""
Compute the graupel PSD's contribution to linear reflectivity using the method of Jung, Zhang, and Xue (2008).
Args:
density_air (np.array): Density of the air (kg/m^3)
rain_psd (PSD): Rain PSD to compute water fractions with.
max_mixed_frac (float): Maximum mixed fraction for the graupel PSD.
Returns:
Graupel contribution to linear reflectivity
"""
# print
# print "Computing reflectivity contribution from graupel"
water_frac, water_mix_mass, graupel_mix_mass, total_mix_mass, density_mix = PSD._waterFractionJZX(rain_psd, self, max_mixed_frac)
qg_wet = graupel_mix_mass
qg_dry = np.maximum(self.mixingrat - qg_wet, 0)
sigma = np.pi / 3 * (1 - 0.8 * water_frac)
A = .125 * (3 + 4 * np.exp(-2 * sigma ** 2) + np.exp(-8 * sigma ** 2))
B = .125 * (3 - 4 * np.exp(-2 * sigma ** 2) + np.exp(-8 * sigma ** 2))
C = .125 * (1 - np.exp(-8 * sigma ** 2))
D = .125 * (3 + np.exp(-8 * sigma ** 2))
Ck = np.exp(-2 * sigma ** 2)
alpha_a = (0.081 + 2.04 * water_frac - 7.39 * water_frac ** 2 + 18.14 * water_frac ** 3 - 26.02 * water_frac ** 4 + 19.37 * water_frac ** 5 - 5.75 * water_frac ** 6) * 1e-3
alpha_b = (0.076 + 1.74 * water_frac - 7.52 * water_frac ** 2 + 20.22 * water_frac ** 3 - 30.42 * water_frac ** 4 + 23.31 * water_frac ** 5 - 7.06 * water_frac ** 6) * 1e-3
Zg_dry = PSD._partialHReflectivityFrozenJZX((A, B, C), (alpha_a, alpha_b), density_air, qg_dry, self.density, self.interceptParam(density_air))
Zg_wet = PSD._partialHReflectivityFrozenJZX((A, B, C), (alpha_a, alpha_b), density_air, qg_wet, density_mix, self.interceptParam(density_air))
return Zg_dry + Zg_wet
[docs]class PSDCollection(object):
"""
Uses several PSDs to compute total statistics for particles in a volume.
Args:
rain (PSD): The rain PSD in the collection. Optional, default is no rain PSD.
snow (PSD): The rain PSD in the collection. Optional, default is no snow PSD.
hail (PSD): The rain PSD in the collection. Optional, default is no hail PSD.
graupel (PSD): The rain PSD in the collection. Optional, default is no graupel PSD.
"""
def __init__(self, rain=None, snow=None, hail=None, graupel=None):
self._psds = dict( (k, v) for k, v in [ ('rain', rain), ('snow', snow), ('hail', hail), ('graupel', graupel)] if v is not None )
return
[docs] def reflectivity(self, *args, **kwargs):
"""
Computes the total logarithmic reflectivity for this PSD collection using different methods. The methods are
described below.
Args:
*args: Differ based on which method you're using. When method='zhang', air temperature (K) and air density
(kg/m^3) are required. For method='jzx', only air density (kg/m^3) is required.
method (string): The method to use when computing reflectivity. 'zhang' uses the Zhang reflectivity
calculations, intended for single-moment microphysics with no graupel category. 'jzx' uses the method
of Jung, Zhang, and Xue (2008).
Returns:
Logarithmic reflectivity (dBZ) computed for all PSDs in this collection.
"""
err = np.geterr()
np.seterr(all='ignore')
if 'method' not in kwargs or kwargs['method'].lower() == "zhang":
refl = self._reflectivityZhang(*args)
elif kwargs['method'].lower() == "jzx":
refl = self._reflectivityJZX(*args)
np.seterr(**err)
return refl
def _reflectivityZhang(self, temperature, density_air):
"""
Computes the total logarithmic reflectivity for this PSD collection using the Zhang method.
Args:
temperature (np.array): Air temperature (K)
density_air (np.array): Density of the air (kg/m^3)
Returns:
Logarithmic reflectivity (dBZ) computed for all PSDs in this collection.
"""
is_wet = np.where(temperature > 273.15, 1, 0)
reflec_factors = [ p.reflectivityZhang(density_air, is_wet) for p in self._psds.itervalues() ]
return 10 * np.log10(np.sum(reflec_factors, axis=0))
def _reflectivityJZX(self, density_air):
"""
Computes the total logarithmic reflectivity for this PSD collection using the method of Jung, Zhang, and Xue
(2008).
Args:
density_air (np.array): Density of the air (kg/m^3)
Returns:
Logarithmic reflectivity (dBZ) computed for all PSDs in this collection.
"""
_fox = {
# Hail, Graupel
(True, True): {'snow':0.2, 'hail':0.25, 'graupel':0.2},
(True, False): {'snow':0.5, 'hail':0.3, 'graupel':0.0},
(False, True): {'snow':0.5, 'hail':0.0, 'graupel':0.3},
(False, False): {'snow':0.5, 'hail':0.0, 'graupel':0.0},
}
max_melting_fracs = _fox[tuple( p in self._psds for p in [ 'hail', 'graupel' ] )]
frozen_psds = [ (v, max_melting_fracs[k]) for k, v in self._psds.iteritems() if k != 'rain' ]
rain_psd = [ v for k, v in self._psds.iteritems() if k == 'rain' ][0]
frozen_reflec_factors = [ p.reflectivityJZX(density_air, rain_psd, mmf) for p, mmf in frozen_psds ]
liquid_reflec_factor = rain_psd.reflectivityJZX(density_air, *frozen_psds)
return 10 * np.log10(np.sum(frozen_reflec_factors, axis=0) + liquid_reflec_factor)
if __name__ == "__main__":
from computeQuantities import computeReflectivity, theta2Temperature
import Nio as nio
from datetime import datetime, timedelta
hdf = nio.open_file("/caps2/tsupinie/05June2009/1kmf-z04vr=30dBZ/ena001.hdf014400", mode='r', format='hdf')
p = hdf.variables['p'][:]
pt = hdf.variables['pt'][:]
temperature = theta2Temperature(p=p, pt=pt)
density_air = p / (temperature * 287.)
start = datetime.now()
old_refl_zhang = computeReflectivity(**dict( (k, hdf.variables[k][:]) for k in ['pt', 'p', 'qr', 'qh', 'qs'] ))
old_zhang_time = datetime.now() - start
psdc = PSDCollection(
rain=RainPSD(hdf.variables['qr'][:], intercept_param=4e5),
snow=SnowPSD(hdf.variables['qs'][:], intercept_param=3e6),
hail=HailPSD(hdf.variables['qh'][:], intercept_param=4e4)
)
start = datetime.now()
new_refl_zhang = psdc.reflectivity(temperature, density_air, method='zhang')
new_zhang_time = datetime.now() - start
start = datetime.now()
refl_jzx_sm = psdc.reflectivity(density_air, method='jzx')
jzx_sm_time = datetime.now() - start
hdf.close()
hdfmm = nio.open_file("/caps2/tsupinie/05June2009/1kmf-my2/ena001.hdf014400", mode='r', format='hdf')
p = hdfmm.variables['p'][:]
pt = hdfmm.variables['pt'][:]
temperature = theta2Temperature(p=p, pt=pt)
density_air = p / (temperature * 287.)
psdcmm = PSDCollection(
rain=RainPSD(hdfmm.variables['qr'][:], num_concentration=hdfmm.variables['nr'][:]),
snow=SnowPSD(hdfmm.variables['qs'][:], num_concentration=hdfmm.variables['ns'][:]),
hail=HailPSD(hdfmm.variables['qh'][:], num_concentration=hdfmm.variables['nh'][:]),
graupel=GraupelPSD(hdfmm.variables['qg'][:], num_concentration=hdfmm.variables['ng'][:])
)
start = datetime.now()
refl_jzx_dm = psdcmm.reflectivity(density_air, method='jzx')
jzx_dm_time = datetime.now() - start
print
print "New Zhang calculation gives the same answer as the old one:", np.all(old_refl_zhang == new_refl_zhang)
print "Max reflectivity from Zhang:", new_refl_zhang.max()
print "Max reflectivity from JZX SM:", refl_jzx_sm.max()
print "Max reflectivity from JZX MM:", np.nanmax(refl_jzx_dm)
print "Time to compute Zhang reflectivity the old way:", old_zhang_time
print "Time to compute Zhang reflectivity the new way:", new_zhang_time
print "Time to compute JZX reflectivity for SM:", jzx_sm_time
print "Time to compute JZX reflectivity for MM:", jzx_dm_time