OLD_NOTES.txt

Fuzzy-Logic Hydrometeor Identification

#### BELOW ARE COPIED COMMENTS AND INFO FROM ORIGINAL IDL/WIENS CODE #####
#### MANY OF THESE ARE NO LONGER APPLICABLE; PROVIDED ONLY FOR POSTERITY #####
#### UPDATED COMMENTS ARE IN DOC STRINGS FOR PYTHON FUNCTIONS #####
#+
# NAME:
# cdf_fhc
#
# PURPOSE:
#  Perform the fuzzy logic hydrometeor classification (FHC)
#  using input polarimetric variables and temperature.  Return the FHC
#  types in an array of the same size as the input polarimetric variables.
#
#  This method is based on the method of Liu and Chandrasekar (2000),
#  J. Atmos. Oceanic Tech., _17_, 140-164.  This assumes summer convective
#  storms.
#  However, my method uses a weighted sum of each truth value instead
#  of the straight multiplication method of L&C (2000).
#
# Truth value of each hydrotype for each input variable is a beta
# function.  The total truth value for each hydro type is thus a
# weighted sum of the truth value of each input variable, with two
# exceptions:  The horizontal reflectivity and temperature (if used) truth
# values are multiplied by the weighted sum of the others, thus effectively
# weighting the final output more strongly on the horizontal reflectivity
# and temperature (if present).
#
# Uses parameters defined in 'cdf_beta_functions.pro' to define the
# m,a,b parameters for each hydro type's beta function.  The independent
# variable 'x' is one of the radar measurements.
#
# Basic form of a Membership Beta Function is:
#
#                         1
#               ---------------------
# beta(x) =         |              |^b
#                1+ |  ((x-m)/a)^2 |
#                   |              |
#
#  Where x = input data (Zh, Zdr, etc)
#        m = center
#        a = width
#        b = slope
#
# Anticipate having up to 6 input data fields and 11 fuzzy sets (HID types) for
# summer storms.
#
# Input data fields:           HID types:           Fuzzy set:
#-------------------           ----------           ----------
#     Z_h                      Drizzle                  1
#     Z_dr                     Rain                     2
#     K_dp                     Dry snow                 3
#     LDR                      Wet snow                 4
#     rho_hv[0]                Vertical Ice             5
#     Temperature              Dry graupel              6
#                              Wet graupel              7
#                              Small hail               8
#                              Large hail               9
#                              Sm. Hail & Rain         10
#                              Lg. Hail & Rain         11
#
#  So there should be 6 X 11 = 66 beta functions...one for each pair
#  of inputs and fuzzy sets
#
# The final output for each grid point is the fuzzy set which has the
# highest cumulative truth value (i.e., the most likely hydrometeor type).
#
# If there is an error, the function will return the integer zero instead
# of a gridded array.
#
# USAGE:
# output = cdf_fhc(radar_data,mbf_sets,fhc_vars[,T,weights=weights,
#                  use_temp=use_temp,compare=compare])
#
# Terms in brackets [] are optional.
#
# INPUT PARAMETERS:
#        REQUIRED INPUTS
# radar_data -- A structure containing the radar data
#               (Should not matter if in radial
#                or gridded format).  The structure must
#               have one of the following tag names: DZ,DR,KD,LD,RH.
#               These correspond to horiz. reflect, diff. reflect.,
#               spec. diff. phase, linear depol. ratio, and corr. coeff.,
#               respectively.  For example, to use all these variables in
#               the classification, define this structure as:
#
#               radar_data = {DZ:dz_data,$
#                             DR:DR_data,$
#                             KD:KD_data,$
#                             LD:LD_data,$
#                             RH:RH_data}
#
#               Then call the program:
#               result = cdf_fhc(radar_data,...)
#
#               Or if you have only horiz. reflect. and
#               diff. reflect. define the structure like this:
#
#               radar_data = {DZ:dz_data,$
#                             DR:DR_data}
#
#               Only the tag names are important, e.g., the two-character
#               string before each colon must be DZ, DR, etc.  The
#               actual data (e.g., dz_data) should be 1,2,or3-d floating
#               point array of the appropriate data.  Each data field is
#               expected to have the same dimensions.
# mbf_sets -- Beta function parameters for each fuzzy set and input
#             variable (defined in 'cdf_beta_functions.pro')
# fhc_vars -- A Boolean structure of which input variables to use in the
#             FHC.  Should be a 6-element structure of either 0 (Don't
#             Include) or 1 (Do Include)
#             fhc_vars = {Zh: 0 or 1,
#                         Zdr:0 or 1,
#                         Kdp:0 or 1,
#                         LDR:0 or 1,
#                         rho_hv: 0 or 1,
#                         T: 0 or 1}
#             The use for this array is to do sensitivity tests where you'd
#             explicity state which radar fields to use in the
#             classification.  If the corresponding data is not present
#             in the radar_data input structure, then this program will
#             override what you passed it as 'fhc_vars' and set
#             its boolean value to 0.  For example, if there is no LD field
#             in the radar_data structure, then this program will do the
#             following:
#
#             fhc_vars.LDR = 0
#
#             And hence will not try to use LDR in the classification.
#
#       OPTIONAL INPUTS
# T -- Temperature, either single value or a 1-D vertical profile with the
#      temperature defined at each vertical grid point of the radar data.
#        --BD 2-2012: Changed to require this be the same size as the input data.
# weights -- Weights for the different variables. Added this as a keyword such
#           that you don't need a speciifc cdf_fhc for different radars.
#           If weights is not passed, then default values are used.
#           This is a structure like fhc_vars:
#        Weights= { W_Zh: 1.5 ,
#                   W_Zdr: 0.8,
#                   W_Kdp: 1.0,
#                   W_rho: 0.8,
#                   W_ldr: 0.5,
#                   W_T: 0.4}
#
#
# KEYWORDS:
# use_temp -- If set, use temperature in the FHC.  If you set this keyword,
#             make sure you also supply a temperature (T) input.
# compare -- If set, a structure will be returned intstead of just the FHT
#            variable.  It will contain the first
#             and second best scoring categories, as well as their scores.
#             fdat = { FHT, FHT2, mu_best, mu_sec}
#
#BD 2-2012: CHANGED SO NONE OF THESE ARE AVAILABLE ANYMORE. By making temperature
#            The same size as the other variables, no longer need /rhi or /ppi.
#            Sum_it never quite worked right for me anyway, so stick with the
#           Hybrid method for now.
# rhi -- IF set, assume temperature is an array of temperatures at
#        each z level of the radar field arrays.  Program assumes that the
#        second dimension of the input radar data corresponds to height.
# ppi -- If set, assume temperature is a single number for a given
#        height level.
# sum_it -- If set, use a weighted sum method for all variables.  If not
#           set, use a weighted sum for Zdr, Kdp, LDR, and rho_hv; then
#           multiply this sum by the scores for Zh and Temperature.
#
# You need to set either the 'rhi' or 'ppi' keyword only if you are passing
# a 2-D cross-section of data to this program and you are including
# temperature.  Otherwise, if you are passing a 3-D array and including
# temperature, the program will assume the 3rd dimension of the input radar
# data corresponds to height.
#
#IF radar_data is just single values, then the scores for each category given
#   those inputs will be printed out.
#
#IF the compare keyword is set, then the output will be a structure that contains
#   the first and second best scoring
#    HID types, as well as the socres associated with those two categories.
#
# AUTHOR:
# Kyle C. Wiens
# kwiens@atmos.colostate.edu
#
# DATE:
# 26 April 2002.
#
# MODIFICATIONS:
#
# 28 April 2002:  Changed the aggregation operations.  Instead of
# taking a sum of all the MBFs, I take a sum of the polarimetric MBFs
# and then multiply this sum by the MBFs for temperature and
# horizontal reflectivity.
#
# 4 May 2002:  Added weighting factors.
#              Also added ability to use any combination of the input
#              variables in the FHC by multiplying each MBF by its
#              associated fhc_vars value.  This is for doing sensitivity
#              tests or if you don't have all 5 of the radar measurements.
#
# 5 May 2002:  Made weighting factors for Zdr, Kdp, LDR, and rho_hv
# dependent on the hydro type.
#
# 15 May 2002: Added division by the sum of the weighting factors so
# that the scores are normalized with respect to these factors.  I
# should have done this a long time ago!!!
#
# 24 May 2002: Changed the weighting factors back to being the same for each
# hydrometeor type.
#
# 23 July 2002: Greatly improved the efficiency of this program by
# performing all the beta function calculations on the entire 2 or 3
# dimensional data sets, all at once.  I used to have 4 embedded FOR
# loops to do this.  It took some brain power to figure out how to
# keep the maximum scores from each iteration, but apparently, my
# brain was up to the task.  This runs MUCH faster now, with identical
# results.
#
# 24 July 2002: Decomposed rain/hail mix into large and small hail
#               components.
#
# 7 April 2003: Added a few lines to allow the user to do the
#               classification using any combination of input variables.
#
# 23 April 2003:  Changed the calling sequence so that the input radar data
# is a structure instead of each individual field.  This should make the
# program more modular by allowing it to work with any combination of radar
# variables.   I also added more comments.
#
# 21 January 2004:  Added option to use a weighted sum for all variables by
#                   passing the /Sum_It keyword.
# 22 February 2012: Modified to be more modular, and no longer require so many
#                   keywords.
#                    Also set up to pass in a single set of radar observations to
#                   calculate the scores.
#                    B.D. 2-2012
#
#--------
# This is a sub-function of the main program.
#
# FUNCTION hid_beta
# Compute the value of the beta function from inputs x,a,b,m
# Where,
#        x = independent variable (DZ,KDP,etc.)
#        a = width.
#        b = slope.
#        m = midpoint
#--------

########## END OF ORIGINAL COMMENTS AND INFORMATION ##########

Blended Rainfall Calculation

# See original comments below
#***************************************************************************
#***************************************************************************
# Based on Code from A. Rowe and T. Lang from L. Nelson and L. Carey.
# Brenda Dolan
# February 27, 2012
# bdolan@atmos.colostate.edu
#***************************************************************************
#***************************************************************************

#############################
ORIGINAL FUNCTION NOTES BELOW
KEPT FOR POSTERITY ONLY
MANY POINTS OUT OF DATE
#############################
# Brody Fuchs, Oct 2014, converting the IDL code to python here
# This program essentially dublicates the fortran program
# hybridrain_mod5.f.  It determine best polarimetric
# rainfall estimator, based on thresholds of ZH, ZDR, and KDP,
# as well as Ice fraction.
# The input files are 2-D cdf (gridded to a specific ht).
#
# The program also does hydro ID, using the appropriate subroutines.
# Because the input cdfs are 2D (1 ht), temp info at the height of
# interest should be input for the event of interest

# fi=ice fraction
# dz=radar reflectivity
# dr=differential reflectivity
# kd=specific differntial phase
# rkd=R(kdp)
# rkddr=R(kdp,zdr)
# rzhdr=R(zh,Zdr)

# create zk "blended" rainfall (zk) array and an array to identify which
# method is used
# 1 = R(kdp,zdr)
# 2 = R(kdp)
# 3 = R(zh,zdr)
# 4 = R-Z NXRAD (total reflectivity)
# 5 = R-Z NEXRAD (rain only reflectivity)

#cdf_array=intarr(3)

### IF CALLING FROM BRORADAR, ALREADY HAVE DATA, DONT NEED TO READ IN
# open the cdf files
#    dir = '/data3/brfuchs'
#cdffile = findfile(dir+'*.cdf',count=numfile)
#    cdffile = '%s/CHL20120606_220558_cedunf.cdf'

#    vars = ['DZ', 'DR', 'KD', 'LH', 'RH']

#    radar_data = {}
#    nc = Dataset(filename)

# READ IN NETCDF FILE AND GET DATA
#    for var in vars:
#       radar_data[var] = nc.variables[var][:]
# create x,y arrays that will be used in the output cdf files
#    lon=xfirst+xdelta*findgen(cdf_array[0])
#    lat=yfirst+ydelta*findgen(cdf_array[1])

# BF have x and y arrays already

# Run Kyle W.'s fuzzy logic HID program

# first call the membership routine that does the fuzzification...
# use the 2004 updated version of the MBFs with the "nomix" option
#   cdf_beta_functions_kwiens,mbf_sets,/use_temp,/NoMix

# now run the HID program
# first have to define a couple of structures
# 1st structure is the radar parameter tags and array
# that the tag corresponds to - this is described in the
# HID subroutine cdf_fhc.pro
#
# for chill rain maps, input cdf files have only DZ,DR, and KD
#   radar_data = {DZ:dz,DR:dr,KD:kd}

# define a structure of which input variables to actually use in the HID
# don't really need this since we want to use all the variables in the HID
# process
#   fhc_vars = {Zh:1, Zdr:1, Kdp:1, LDR:0, rho_hv:0,T:1}

# CHILL rainmaps are 2D at 1 km height (CAPPI or ~ 1km height for
# SPRINT  gridppi option).  So,input a T at 1km AGL for event of
# interest (based on sounding file)
#   T=dz
#   T[*,*]=17.8

# now run the HID program over the grid and stuff the results into a
# new array (hid)
#   hid= cdf_fhc(radar_data,mbf_sets,fhc_vars,T,/use_temp,/ppi)

# here's the number of hids - useful for defining more arrays...
#   nhid=n_elements(mbf_sets.zh_set.a)
    nhid = 10

# set and initialize arrays for rr estimation

#############################
END ORIGINAL FUNCTION NOTES
#############################

#############################
ORIGINAL FUNCTION NOTES BELOW
KEPT FOR POSTERITY ONLY
MANY POINTS OUT OF DATE
#############################
###### HID rain rate calculations ########

#This is a blended polarimetric rainfall algorithm that selects
#rainfall estimators based on a hydrometeor identification scheme.
#
#Modified by Brenda Dolan
#bdolan@atmos.colostate.edu
#June 2010
#Colorado State University

#Ported over to python by Brody Fuchs, Oct 2014

# first step is to fit the consensus of hid from all 6 inputs.  There
# are 3 categories to consider: rain, mix, ice. the
# category with the biggest number wins.
# 1=rain(c_hid_consensus=1)
# 2=mix(c_hid_consensus=2)
# 3=ice(c_hid_consensus=3)
# 4=rain-mix tie (average of rain rr and mix rr is used)
# (c_hid_consensus=4)
# 5 rain-ice tie (rain value is used in rr calculation)
# (c_hid_consensus=5)
# 6 mix-ice tie (mix value is used in rr calculation)
# (c_hid_consensus=6)

# this is a variable that keeps track of which hid is actually
# used for the site in question, based on the relative totals
# of the 6 values.
#   c_hid_consensus=dz
#   c_hid_consensus[*]=0.

# here is variable to store the actual rr estimate method for hidzk
# the values go from 1-5 just like for zk; however, there is also a method
# 6 which represents the result when mix and rain HID end in a tie:
# in this case, an average between the value calculated using the mix (r(kd))
# and rain methods (could be r(kd), r(kd,dr), r(zh) or r(zh,dr))
#############################
END ORIGINAL FUNCTION NOTES
#############################

Liquid/Ice Mass Calculations

################################
OLD FUNCTION NOTES BELOW
MAINTAINED FOR POSTERITY'S SAKE, MAY BE OUT OF DATE
################################
#In this case, dbz and zdr are 3-D arrays. Z is an array of the vertical grid.
# CHECK TO SEE IF THIS IS IN GRAMS OR KILOGRAMS????
# for comparison: Brett is getting order 10^9 kg for ice mass
#    delta_thresh = fltarr(n_elements(z))
#This uses methodology similar to Larry Carey's that makes it more and more
#restrictive to identify ice below the melting level...based upon personal
#communication with Larry.  Ice and Water mass are calculated using the Zdp
#methodology outlined in Bringi and Chandra (2001)
# Define empirical thresholds at each vertical level to restrict ice below
#   melting level
# These delta thresholds are set up for a 0.5 km grid. These are to make it more
#   difficult to get ice in the lowest levels of the grid.
#                delta_thresh[*] = 2.0
#                delta_thresh[0] = 11.3333  ; 0.5 km
#                delta_thresh[1] = 10.0000  ; 1.0 km
#                delta_thresh[2] = 8.66667  ; 1.5 km
#                delta_thresh[3] = 7.33333  ; 2.0 km
#                delta_thresh[4] = 6.00000  ; 2.5 km
#                delta_thresh[5] = 4.66667  ; 3.0 km
#                delta_thresh[6] = 3.33333  ; 3.5 km
#
#    #These delta thresholds are set up for a 1.0 km grid starting at 0 km.
#    delta_thresh = np.zeros_like(z)
#    delta_thresh[:] = 1.1
#    delta_thresh[0] = 12.6667 #  0 km
#    delta_thresh[1] = 10.0000 #  1 km
#    delta_thresh[2] = 7.33333 #  2 km
#    delta_thresh[3] = 4.66667 #  3 km
#    delta_thresh[4] = 2.0  #     4 km
################################
#END OLD NOTES
################################