I calculate surface visibility from WRF output using a calculation that I adapted from DTC's Unified Post Processor, specifically from their Fortran routine found in UPPV2.2/src/unipost/CALVIS.f. The calculation is based on hydrometeor mixing ratios, and air temperature and pressure, all from the lowest model layer. The documentation in the original code reads:
This routine computes horizontal visibility at the
surface or lowest model layer, from qc, qr, qi, and qs.
qv--water vapor mixing ratio (kg/kg)
qc--cloud water mixing ratio (kg/kg)
qr--rain water mixing ratio (kg/kg)
qi--cloud ice mixing ratio (kg/kg)
qs--snow mixing ratio (kg/kg)
tt--temperature (k)
pp--pressure (Pa)
If iice=0:
qprc=qr qrain=qr and qclw=qc if T>0C
qcld=qc =0 =0 if T<0C
qsnow=qs and qclice=qc if T<0C
=0 =0 if T>0C
If iice=1:
qprc=qr+qs qrain=qr and qclw=qc
qcld=qc+qi qsnow=qs and qclice=qc
Independent of the above definitions, the scheme can use different
assumptions of the state of hydrometeors:
meth='d': qprc is all frozen if T<0, liquid if T>0
meth='b': Bocchieri scheme used to determine whether qprc
is rain or snow. A temperature assumption is used to
determine whether qcld is liquid or frozen.
meth='r': Uses the four mixing ratios qrain, qsnow, qclw,
and qclice
The routine uses the following
expressions for extinction coefficient, beta (in km**-1),
with C being the mass concentration (in g/m**3):
cloud water: beta = 144.7 * C ** (0.8800)
rain water: beta = 2.24 * C ** (0.7500)
cloud ice: beta = 327.8 * C ** (1.0000)
snow: beta = 10.36 * C ** (0.7776)
These expressions were obtained from the following sources:
for cloud water: from Kunkel (1984)
for rainwater: from M-P dist'n, with No=8e6 m**-4 and
rho_w=1000 kg/m**3
for cloud ice: assume randomly oriented plates which follow
mass-diameter relationship from Rutledge and Hobbs (1983)
for snow: from Stallabrass (1985), assuming beta = -ln(.02)/vis
The extinction coefficient for each water species present is
calculated, and then all applicable betas are summed to yield
a single beta. Then the following relationship is used to
determine visibility (in km), where epsilon is the threshhold
of contrast, usually taken to be .02:
vis = -ln(epsilon)/beta [found in Kunkel (1984)]
I have adapted the code from this routine to a Python function below, which you can use for your purposes:
def calculate_visibility(qv,qc,qr,qi,qs,T,p):
"""
Calculates visibility based on the UPP algorithm.
See documentation in UPPV2.2/src/unipost/CALVIS.f for the description of
input arguments and references.
"""
Rd = 287.
COEFLC = 144.7
COEFLP = 2.24
COEFFC = 327.8
COEFFP = 10.36
EXPLC = 0.88
EXPLP = 0.75
EXPFC = 1.
EXPFP = 0.7776
Tv = T * (1+0.61*qv) # Virtual temperature
rhoa = p/(Rd*Tv) # Air density [kg m^-3]
rhow = 1e3 # Water density [kg m^-3]
rhoi = 0.917e3 # Ice density [kg m^-3]
vovmd = (1+qv)/rhoa + (qc+qr)/rhow + (qi+qs)/rhoi
conc_lc = 1e3*qc/vovmd
conc_lp = 1e3*qr/vovmd
conc_fc = 1e3*qi/vovmd
conc_fp = 1e3*qs/vovmd
# Make sure all concentrations are positive
conc_lc[conc_lc < 0] = 0
conc_lp[conc_lp < 0] = 0
conc_fc[conc_fc < 0] = 0
conc_fp[conc_fp < 0] = 0
betav = COEFFC*conc_fc**EXPFC\
+ COEFFP*conc_fp**EXPFP\
+ COEFLC*conc_lc**EXPLC\
+ COEFLP*conc_lp**EXPLP+1E-10
vis = -np.log(0.02)/betav # Visibility [km]
vis[vis > 24.135] = 24.135
return vis