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