# -*- coding: utf-8 -*-
"""
Created on Thu Jan 28 14:52:44 2021

@author: ppantina
"""

##Import the required libraries
import glob
import sys
import os
import h5py
import datetime
import numpy          as     np
from matplotlib       import pyplot as plt
from scipy            import ndimage  #for convolutions
from silx.io.dictdump import dicttoh5 #for hdf writing       
from datetime         import timezone
import L1B_sub_add_horiz_wind_hrrr_fullflight    #to calculate winds
import L1B_sub_along_cross_winds2     #to calculate along/cross winds
import L1B_sub_plotter
import sub_cmap        #for colorbar
import sub_params      #for parameters
import sub_flightLines #for flightLines
import L1B_sub_EXRAD_SCAN
##Change these variables as needed
radName       = input('Which radar? CRS, EXRAD, HIWRAP_KU, HIWRAP_KA: ')
flightDate    = input('Which date? yyyymmdd (append flight number if needed) ')
versionLetter = input('Which version?: A, B, C, ... ')
makePlots     = input('Run the plots?: yes, no, ... ')


##Next, read the parameter file to fine radar-relevant variables.
fileNameParams = 'radar_params_irma.xlsx'
config = sub_params.sub_params(fileNameParams, flightDate, radName)

##Define the colorbar
cmap_n0q = sub_cmap.make_cmap('n0q', bit = True)

######Define paths and filenames######
##First locate L1A paths/files
L1A_dir  = config['Directory_L1A']
if radName== 'EXRAD_SCAN':
    L1A_file = np.sort(glob.glob(L1A_dir + '*' + radName + '*' + flightDate + '*.h5'))
    L1A      = h5py.File(L1A_file[0], 'r')

else:
    L1A_file = glob.glob(L1A_dir + '*' + radName + '*' + flightDate + '*.h5')[0]
    print(L1A_file)
    L1A      = h5py.File(L1A_file, 'r')

#

##Next locate L1B paths/files. Create the output dir if necessary.
L1B_dir      = config['Directory_L1B']
L1B_file     = config['experimentName'] + '_' + radName + '_L1B_Rev'+ versionLetter + '_' + flightDate + '.h5'
L1B_filepath = L1B_dir + L1B_file
L1B_plotpath = L1B_dir + 'plots/' + radName + '/' + flightDate + '/'
if not(os.path.isdir(os.path.dirname(L1B_filepath))): os.makedirs(os.path.dirname(L1B_filepath)) ##make the output dir if needed.
if not(os.path.isdir(os.path.dirname(L1B_plotpath))): os.makedirs(os.path.dirname(L1B_plotpath)) ##make the output dir if needed.

##Finally locate the HRRR filepath
m_filepath   = config['Directory_hrrr'] + '/nc3d/*.nc'

##Define some vars to index the ranges.
rangeMin    = config['rangeMin'   ]
rangeMax    = config['rangeMax'   ]
deg_offset  = config['degOffset'  ] ##roll_offset for wind correction.
##calCoeffs defined separately. For KA, there is a second offset.
##this is a roundabout way to add tuples.
extCalCoeff = np.array(eval(config['extCalCoeff'])) + np.array(eval(config['extCalCoeff2']))

##Set up range indices and calibration coefficient arrays.
##For HIWRAP, we need a dictionary of ranges, one for each channel
if radName == 'HIWRAP_KU':
    #% The -0.10 (Ku) -0.33 (Ka) is for attenuation in the external calibration
    #% The 0.1, +0.2 (Ka) is for calibration matching.
    #%L1B RevA calibration Modification - Lowering Ka-band by 1.25 dB.  
    rangeIndex = {}
    for ch in np.arange(0,4): rangeIndex[ch] = np.where((L1A['Products/Channel' + str(ch) + '/Information/Range'][:]>rangeMin) & (L1A['Products/Channel' + str(ch) + '/Information/Range'][:]<rangeMax))[0]
if radName == 'HIWRAP_KA':
    rangeIndex = {}
    for ch in np.arange(0,4): rangeIndex[ch] = np.where((L1A['Products/Channel' + str(ch) + '/Information/Range'][:]>rangeMin) & (L1A['Products/Channel' + str(ch) + '/Information/Range'][:]<rangeMax))[0]

##Read in the appropriate flight lines.
flightLines = sub_flightLines.sub_flightLines (flightDate)


#######################Add /Information to h5#######################
##To write to dict, all strings must be placed in arrays, as far as I can tell.
print ('Adding Information')
L1B                = {}
L1B['Information'] = {}
L1B['Information/ExperimentName'   ] = L1A['ExperimentName' ]
L1B['Information/L1A_ProcessDate'  ] = L1A['L1A_ProcessDate']
L1B['Information/FlightDate'       ] = [flightDate          ]
L1B['Information/RadarName'        ] = [radName             ]
L1B['Information/L1B_ProcessDate'  ] = [datetime.datetime.now().strftime("%Y%m%d")]
L1B['Information/L1B_Revision'     ] = ['Rev' + versionLetter]
L1B['Information/DataContact'      ] = [config['dataContact' ]]
L1B['Information/MissionPI'        ] = [config['missionPI'   ]]
L1B['Information/Aircraft'         ] = [config['aircraft'    ]]
L1B['Information/InstrumentPI'     ] = [config['instrumentPI']]
L1B['Information/L1B_Revision_Note'] = [config['revisionNote']]

##For HIWRAP, save data every 0.5 seconds instead of every 0.25. It is too large otherwise. For other radars, save all times.
if 'HIWRAP' in radName: times2save = np.where((L1A['Time/Data/TimeUTC'][:]%1 == 0) | (L1A['Time/Data/TimeUTC'][:]%1 == 0.5))[0]
else:                   times2save = np.where( L1A['Time/Data/TimeUTC'][:] > 0)[0]

if ((radName == 'HIWRAP_KU')&(flightDate == '20220203')):
    times2save = np.where(((L1A['Time/Data/TimeUTC'][:]%1 == 0) | (L1A['Time/Data/TimeUTC'][:]%1 == 0.5)) & (L1A['Time/Data/TimeUTC'][:]>1643902120))[0]
    flightLines  = flightLines[-2:] ## Only the last two legs are good.
    print ('Special flightLine case')   

#######################Add Time/Data and Time/Information#######################
print ('Adding Time')
flightYear = flightDate[0:4] #extract the year from flightDate. Use it for printing below.
L1B['Time']             = {}
L1B['Time/Data']        = {}
L1B['Time/Information'] = {}
L1B['Time/Data/TimeUTC'                          ] = L1A['Time/Data/TimeUTC'][:][times2save] ##Select only the relevant times
L1B['Time/Information/TimeUTC_units'             ] = ['seconds']
L1B['Time/Information/TimeUTC_description'       ] = L1A['Time/Information/Description']
L1B['Time/Information/TimeUTC_01Jan' + flightYear] = [datetime.datetime(int(flightYear),1,1).replace(tzinfo=timezone.utc).timestamp()]
L1B['Time/Information/TimeUTC_01Jan' + flightYear + '_description'] = ['Time of 0 UTC, Jan 01, ' + flightYear + ', for reference if the user does not have an easy linux time converter']

##Define the channel names and frequencies we are using for HIWRAP.
##Could use np.arange(3) for the channel numbers, but explicitly writing them in case they change.
##For KA, channels at L0 are 4, 5, 6. They have been reindexed for L1A/L1B.
if radName == 'HIWRAP_KU':
    chNum  = [ 0,            1,             2     ]
    chFrq  = ['Ku',         'Ku',          'Ku'   ]
    chName = ['LowResPulse','HighResPulse','Chirp']
if radName == 'HIWRAP_KA':
    chNum  = [  0,            1,             2     ]
    chFrq  = [ 'Ka',         'Ka',          'Ka'   ]
    chName = [ 'LowResPulse','HighResPulse','Chirp']

#######################Add Products/Information#######################
print ('Adding Product Information')
L1B['Products'            ] = {}
#L1B['Products/Data'       ] = {}
L1B['Products/Information'] = {}
##HIWRAP requires Information for both frequencies
if radName == 'HIWRAP':
    #L1B['Products/Data'       ] = {} ##We dont actually need this for HIWRAP
    for k, ch in enumerate (chNum): 
        ##Create empty dicts for Chirp/Low/High
        L1B['Products/' + chName[k]]                 = {}    
        L1B['Products/' + chName[k]+ '/Information'] = {}
        L1B['Products/' + chName[k]+ '/Data'       ] = {}

##Redefine the ranges (channels are different between radars), then calculate GateSpacing
if 'HIWRAP' in radName       :L1B['Products/Information/Range'] = L1A['Products/Channel0/Information/Range'][:][rangeIndex[0]] #choose first range index
ranges = np.copy(L1B['Products/Information/Range'])            ##Use this as a shortcut later
L1B['Products/Information/GateSpacing'] = np.diff(ranges[0:2]) ##Just the difference between first two gates

##Populate aircraft motion. Use the smoothed motion where available.
if 'HIWRAP' in radName:
    L1B['Products/Information/AircraftMotion'            ] = L1A['Products/Information/AircraftMotion_smooth'][:][times2save]
    L1B['Products/Information/AircraftMotion_units'      ] = ['m/s']
    L1B['Products/Information/AircraftMotion_description'] = ['Estimated aircraft motion normal to the beam subtracted from Doppler estimate. Smoothed to a 1.5-second average motion.']

##Added turn flag Dec 2022
turnFlag = np.zeros_like(L1B['Time/Data/TimeUTC'][:], dtype = np.int8)
for ff in flightLines:
    idx = np.where((L1B['Time/Data/TimeUTC'][:]>ff[0])&(L1B['Time/Data/TimeUTC'][:]<ff[1]))
    turnFlag[idx] = 1
L1B['Products/Information/AircraftTurnFlag'            ] = turnFlag
L1B['Products/Information/AircraftTurnFlag_units'      ] = ['0 or 1']
L1B['Products/Information/AircraftTurnFlag_description'] = ['Flag is 1 when plane is flying level without turns. Flag is 0 elsewhere.']

##These are common to all radars
L1B['Products/Information/PRI'               ] = [str(int(config['L1B_HighPRF']*1e6)) + ' us / '+str(int(config['L1B_LowPRF']*1e6))+' us staggered']
L1B['Products/Information/AveragedPulses'    ] = np.int32  ([config['L1B_averagedPulses'       ]])
L1B['Products/Information/AntennaBeamwidth'  ] = np.float32([config['L1B_beamWidth'            ]])
L1B['Products/Information/AntennaSize'       ] = np.float32([config['AntennaSize'              ]])
L1B['Products/Information/Frequency'         ] = L1A['Products/Channel0/Information/Frequency'][:]
L1B['Products/Information/Wavelength'        ] = 3e8/L1B['Products/Information/Frequency'     ][:]

##Populate pulses, frequency, wavelength units/descriptions
L1B['Products/Information/AntennaBeamwidth_units'      ] = ['degrees']
L1B['Products/Information/AntennaBeamwidth_description'] = ['Antenna one-way 3dB beamwidth']
L1B['Products/Information/AntennaSize_units'           ] = ['meters']
L1B['Products/Information/AntennaSize_description'     ] = ['none'  ]
L1B['Products/Information/AveragedPulses_units'        ] = ['unitless']
L1B['Products/Information/AveragedPulses_description'   ] = ['Number of pulses averaged per profile, centered around profile time']
L1B['Products/Information/Frequency_units'             ] = ['Hz']
L1B['Products/Information/Frequency_description'       ] = ['Radar transmit frequency']
L1B['Products/Information/GateSpacing_units'           ] = ['meters']
L1B['Products/Information/GateSpacing_description'     ] = ['Range gate spacing']
L1B['Products/Information/PRI_units'                   ] = ['microseconds']
L1B['Products/Information/PRI_description'             ] = ['Pulse repetition interval.']
L1B['Products/Information/Range_units'                 ] = ['meters']
L1B['Products/Information/Range_description'           ] = ['Range in meters from the aircraft']
L1B['Products/Information/Wavelength_description'      ] = ['Radar transmit wavelength']
L1B['Products/Information/Wavelength_units'            ] = ['meters']
L1B['Products/Information/NominalAntennaPointing'      ] = ['Nadir' ]

#######################Populate the Navigation/Data #######################
print ('Adding Navigation')
L1B['Navigation'            ] = {}
L1B['Navigation/Data'       ] = {}
L1B['Navigation/Information'] = {}

##Add navigation data
L1B['Navigation/Data/Drift'        ] = L1A['Navigation/Data/Drift'        ][:][times2save]
L1B['Navigation/Data/EastVelocity' ] = L1A['Navigation/Data/EastVelocity' ][:][times2save]
L1B['Navigation/Data/Heading'      ] = L1A['Navigation/Data/Heading'      ][:][times2save]
L1B['Navigation/Data/Height'       ] = L1A['Navigation/Data/Height'       ][:][times2save]
L1B['Navigation/Data/Longitude'    ] = L1A['Navigation/Data/Longitude'    ][:][times2save]
L1B['Navigation/Data/NorthVelocity'] = L1A['Navigation/Data/NorthVelocity'][:][times2save]
L1B['Navigation/Data/Latitude'     ] = L1A['Navigation/Data/Latitude'     ][:][times2save]
L1B['Navigation/Data/Pitch'        ] = L1A['Navigation/Data/Pitch'        ][:][times2save]
L1B['Navigation/Data/Roll'         ] = L1A['Navigation/Data/Roll'         ][:][times2save]
L1B['Navigation/Data/Track'        ] = L1A['Navigation/Data/Track'        ][:][times2save]
L1B['Navigation/Data/UpVelocity'   ] = L1A['Navigation/Data/UpVelocity'   ][:][times2save]
L1B['Navigation/Data/dxdr'         ] = L1A['Navigation/Data/dx'           ][:][times2save]
L1B['Navigation/Data/dydr'         ] = L1A['Navigation/Data/dy'           ][:][times2save]
L1B['Navigation/Data/dzdr'         ] = L1A['Navigation/Data/dz'           ][:][times2save]

##Calculate speed and nominal distance
spd = np.abs(L1B['Navigation/Data/EastVelocity'][:] + 1j*L1B['Navigation/Data/NorthVelocity'][:])
spd[np.isnan(spd)] = np.nanmedian(spd) #replace missing values with median speed
L1B['Navigation/Data/NominalDistance'] = np.float32(np.diff(L1B['Time/Data/TimeUTC'][0:2])*np.cumsum(spd)) #distance = rate * time

##Add units and information
L1B['Navigation/Information/PitchCorrection'              ] = L1A['Navigation/Information/PitchCorrection'                ]
L1B['Navigation/Information/RollCorrection'               ] = L1A['Navigation/Information/RollCorrection'                 ] 
L1B['Navigation/Information/HeadingCorrection'            ] = L1A['Navigation/Information/HeadingCorrection'              ]
L1B['Navigation/Information/PitchCorrection_units'        ] = ['degrees'                                                  ]
L1B['Navigation/Information/PitchCorrection_description'  ] = ['Pitch offset applied relative to navigation data source'  ]
L1B['Navigation/Information/RollCorrection_units'         ] = ['degrees'                                                  ]
L1B['Navigation/Information/RollCorrection_description'   ] = ['Roll offset applied relative to navigation data source'   ]
L1B['Navigation/Information/HeadingCorrection_units'      ] = ['degrees'                                                  ]
L1B['Navigation/Information/HeadingCorrection_description'] = ['Heading offset applied relative to navigation data source']
L1B['Navigation/Information/NavigationSource'             ] = [config['L1A_navSource']+' Navigation data corrected to nominal radar antenna pointing' ]
L1B['Navigation/Information/Drift_units'                  ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Drift_description'            ] = ['Difference between track and heading'                                                 ]
L1B['Navigation/Information/EastVelocity_units'           ] = ['m/s'                                                                                  ]
L1B['Navigation/Information/EastVelocity_description'     ] = ['Eastward portion of velocity'                                                         ]
L1B['Navigation/Information/Heading_units'                ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Heading_description'          ] = ['Aircraft heading in degrees from north. 90 degrees is eastward pointing'              ]
L1B['Navigation/Information/Height_units'                 ] = ['meters'                                                                               ]
L1B['Navigation/Information/Height_description'           ] = ['Aircraft height above sea level'                                                      ]
L1B['Navigation/Information/NorthVelocity_units'          ] = ['m/s'                                                                                  ]
L1B['Navigation/Information/NorthVelocity_description'    ] = ['Northward portion of velocity'                                                        ]
L1B['Navigation/Information/Latitude_units'               ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Latitude_description'         ] = ['Latitude'                                                                             ]
L1B['Navigation/Information/Longitude_units'              ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Longitude_description'        ] = ['Longitude'                                                                            ]
L1B['Navigation/Information/Pitch_units'                  ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Pitch_description'            ] = ['Pitch'                                                                                ]
L1B['Navigation/Information/Roll_units'                   ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Roll_description'             ] = ['Roll'                                                                                 ]
L1B['Navigation/Information/Track_units'                  ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Track_description'            ] = ['Direction from motion in degrees from north. 90 degrees is eastward motion'           ]
L1B['Navigation/Information/UpVelocity_units'             ] = ['m/s'                                                                                  ]
L1B['Navigation/Information/UpVelocity_description'       ] = ['Upward Velocity'                                                                      ]
L1B['Navigation/Information/dxdr_units'                   ] = ['m/m'                                                                                  ]
L1B['Navigation/Information/dxdr_description'             ] = ['Data cross-track distance from aircraft per radar range. Positive is in the starboard direction.']
L1B['Navigation/Information/dydr_units'                   ] = ['m/m'                                                                                             ]
L1B['Navigation/Information/dydr_description'             ] = ['Data along-track distance from aircraft per radar range. Positive is in the forward direction.'  ]
L1B['Navigation/Information/dzdr_units'                   ] = ['m/m'                                                                                             ]
L1B['Navigation/Information/dzdr_description'             ] = ['Data vertical distance from aircraft per radar range. Positive is in the upward direction.'      ]    
L1B['Navigation/Information/NominalDistance_units'        ] = ['meters'                                                                                          ]
L1B['Navigation/Information/NominalDistance_description'   ] = ['Nominal aircraft travel distance estimated by cumulative summing instantaneous velocities.'      ]

##Define subroutine to calculate horizontal resolution based on beamwidth and range
def res (ranges, config):
    winHann = lambda N: 0.5 * (1-np.cos(2*np.pi*np.arange(1,N+1)/(N+1)))
    import numpy as np
    res = np.zeros(len(ranges), dtype = np.float32)  
    for r in range (len(ranges)): ##for all ranges
        #%x = conv(nominalmotion,antennabeam^2)
        #%res = x>0.25
        x    = np.convolve(np.ones(int(np.floor(200*(config['L1B_LowPRF']+config['L1B_HighPRF'])/2.*config['L1B_averagedPulses']))), winHann(np.floor(2*ranges[r]*np.sin(np.deg2rad(config['L1B_beamWidth']))))**2)
        x    = x/np.max(x) ##normalize to the max
        res[r] = np.sum(x>0.25)
        if r > 1:
            if  res[r] < res[r-1]: res[r] = res[r-1]
    return(res)

##Populate the horiz/vertical resolutions
horiz_res = res(ranges, config)
L1B['Products/Information/ResolutionHorizontal6dB'               ] = horiz_res
L1B['Products/Information/ResolutionHorizontal6dB_units'         ] = ['meters']
L1B['Products/Information/ResolutionHorizontal6dB_description'   ] = ['Approximate horizontal resolution defined as width of spatial weighting after averaging as a function of radar range']

##Subroutine to create a single channel mask variable from three individual L1A variables
def mask_sigma(L1A, channelNumber, rangeIndex):
    import numpy as np    
    array = np.zeros(L1A['Products/Channel' + channelNumber + '/Information/Mask1Sigma'].shape, dtype = 'int8')
    array[L1A['Products/Channel' + channelNumber + '/Information/Mask1Sigma'][:]==1] = 1
    array[L1A['Products/Channel' + channelNumber + '/Information/Mask2Sigma'][:]==1] = 2
    array[L1A['Products/Channel' + channelNumber + '/Information/Mask3Sigma'][:]==1] = 3
    array = array[:,rangeIndex] ##cut by rangeIndex
    return(array)
    
##Populate CoPol/CrPol Masks
if 'HIWRAP' in radName: ##CoPol only, different channels
    for k, ch in enumerate(chNum): 
        L1B['Products/' + chName[k]+ '/Information/MaskCoPol'            ] = mask_sigma(L1A, str(ch), rangeIndex[ch])[times2save]
        L1B['Products/' + chName[k]+ '/Information/MaskCoPol_description'] = ['Co-polarization SNR Mask: (Mask >= #) corresponds with (SNR > # sigma noise)']

#Function to calculate dBZe from Zuncal
def dbze  (zuncal, extcalCoeff, wavelength, k_factor): 
    import numpy as np
    db = (10*np.log10(zuncal) + extcalCoeff + 10*np.log10(wavelength**4*10**18/np.pi**5/k_factor)).astype(np.float32)
    return (db)

#Subroutine to process sigma0 from NRCSuncal
def sigma (NRCSuncal, extCalCoeff):
    import numpy as np
    NRCSuncal[NRCSuncal <= 0] = np.nan               #Filter negative numbers
    NRCSuncal = 10*np.log10(NRCSuncal) + extCalCoeff #convert to dB, add external cal coefficient
    return (NRCSuncal)

#######################ADD Products/Data#######################
##For CRS/EXRAD/HIWRAP(uncombined), add Product/Data, using above subroutines when needed
##For EXRAD and HIWRAP Velocities, remove the L1A AircraftMotion and add the smoothed L1B AircraftMotion
print ('Adding Product Data')
if 'HIWRAP' in radName:
    for k, ch in enumerate(chNum):     
        L1B['Products/' + chName[k]+ '/Data/dBZe'                 ] = (dbze (L1A['Products/Channel' + str(ch) + '/Data/Zuncal'] [:][times2save][:,rangeIndex[ch]], extCalCoeff[k], L1B['Products/Information/Wavelength'], config['L1B_Kfactor'])).astype('float32')
        L1B['Products/' + chName[k]+ '/Data/SpectrumWidth'        ] =       (np.real(L1A['Products/Channel' + str(ch) + '/Data/SpecWid'][:][times2save][:,rangeIndex[ch]])).astype('float32')
        L1B['Products/' + chName[k]+ '/Data/sigma0'               ] = sigma(L1A['Products/Channel' + str(ch) + '/Data/NRCSuncal']      [:][times2save], extCalCoeff[0])
        L1B['Products/' + chName[k]+ '/Data/Velocity_uncorrected' ] =       (L1A['Products/Channel' + str(ch) + '/Data/V']              [:][times2save][:,rangeIndex[ch]]).astype('float32')
        L1B['Products/' + chName[k]+ '/Information/SNR'           ] =       (L1A['Products/Channel' + str(ch) + '/Information/SNR'      ][:][times2save][:,rangeIndex[ch]]).astype('float32')
        L1B['Products/' + chName[k]+ '/Information/noiseFloor'    ] =       L1A['Products/Channel' + str(ch) + '/Information/EstimatedNoise'][:][times2save]

##Filter by 1 sigma for dBZe, Velocity, SpectrumWidth
if 'HIWRAP' in radName:
    for k, ch in enumerate(chNum):
        L1B['Products/' + chName[k]+ '/Data/dBZe'                ][L1B['Products/' +chName[k]+ '/Information/MaskCoPol']<1] = np.nan
        L1B['Products/' + chName[k]+ '/Data/Velocity_uncorrected'][L1B['Products/' +chName[k]+ '/Information/MaskCoPol']<1] = np.nan
        L1B['Products/' + chName[k]+ '/Data/SpectrumWidth'       ][L1B['Products/' +chName[k]+ '/Information/MaskCoPol']<1] = np.nan

##Add units and descriptions
if 'HIWRAP' in radName: 
    for k, ch in enumerate(chNum): 
        L1B['Products/' + chName[k]+ '/Information/dBZe_units'                      ] = ['10*log10(mm^6/m^3)']
        L1B['Products/' + chName[k]+ '/Information/dBZe_description'                ] = ['Equivalent reflectivity factor in dB with 1-sigma noise threshold applied. |K|^2 = ' + str(config['L1B_Kfactor'])]
        L1B['Products/' + chName[k]+ '/Information/Velocity_uncorrected_units'      ] = ['m/s']
        L1B['Products/' + chName[k]+ '/Information/Velocity_uncorrected_description'] = ['Doppler velocity with aircraft motion correction and 1-sigma noise threshold applied. Positive velocity is upward. Ka pitch offset has been removed.']
        L1B['Products/' + chName[k]+ '/Information/SpectrumWidth_units'             ] = ['m/s']
        L1B['Products/' + chName[k]+ '/Information/SpectrumWidth_description'       ] = ['Doppler velocity spectrum width estimate including aircraft motion. 1-sigma noise threshold applied.']
        L1B['Products/' + chName[k]+ '/Information/sigma0_units'                    ] = ['10*log10(m^2/m^2)' ]
        L1B['Products/' + chName[k]+ '/Information/sigma0_description'              ] = ['Ocean Normalized Radar Cross Section, only valid over ocean']
        L1B['Products/' + chName[k]+ '/Information/SNR_units'                       ] = ['Watt / Watt']
        L1B['Products/' + chName[k]+ '/Information/SNR_description'                 ] = ['Estimated Signal to Noise Ratio'    ]
        L1B['Products/' + chName[k]+ '/Information/noiseFloor_units'                ] = ['Relative power']
        L1B['Products/' + chName[k]+ '/Information/noiseFloor_description'          ] = ['Uncalibrated noise floor estimate'    ]
##For HIWRAP, add Products/Data for Combined channel
if 'HIWRAP' in radName:
    ChirpSidelobes = np.zeros(118*2+1)
    ChirpSidelobes[0:2    ] = 0
    ChirpSidelobes[2:19   ] = np.linspace(-90,-75,len(range(2,  19)))
    ChirpSidelobes[18:37  ] = np.linspace(-75,-75,len(range(18, 37)))
    ChirpSidelobes[36:41  ] = np.linspace(-75,-67,len(range(36, 41)))
    ChirpSidelobes[40:102 ] = np.linspace(-67,-58,len(range(40, 102)))
    ChirpSidelobes[101:105] = np.linspace(-58,-47,len(range(101,105)))
    ChirpSidelobes[104:108] = np.linspace(-47,-47,len(range(104,108)))
    ChirpSidelobes[-108:]   = ChirpSidelobes[0:108][::-1]
    
    if ((int(flightDate[0:4]) >= 2022) & (radName == 'HIWRAP_KU')): ##Added this to correct blending error Jun 2022. Updated on 20230614 to include 2023+ campaigns
      ChirpSidelobes = np.zeros(118*2+1)
      #ChirpSidelobes3[0:2    ] = 0
      ChirpSidelobes[0:19   ] = np.linspace(-80,-70,len(range(0,  19)))
      ChirpSidelobes[18:41  ] = np.linspace(-70,-67,len(range(18, 41)))
      ChirpSidelobes[40:102 ] = np.linspace(-67,-58,len(range(40, 102)))
      ChirpSidelobes[101:105] = np.linspace(-58,-47,len(range(101,105)))
      ChirpSidelobes[104:108] = np.linspace(-47,-47,len(range(104,108)))
      ChirpSidelobes[-108:]   = ChirpSidelobes[0:108][::-1]

      
    ChirpSidelobes = 10**(ChirpSidelobes/10)
    ChirpSidelobes[ChirpSidelobes==1] = 0
    ChirpSidelobes = ChirpSidelobes*10**(-5.71/10)#%ratio of peak to mainlobe power

    ChirpSNR = L1A['Products/Channel2/Information/SNR'][:][times2save][:,rangeIndex[2]]
  
    ChirpEstSL = (ndimage.filters.convolve1d(ChirpSNR,ChirpSidelobes,axis=1,mode='constant'))
    ChirpEstSL[ChirpEstSL<=0] = np.nan
    ChirpSNR  [ChirpSNR  <=0] = np.nan
    Chirpcombmask = np.zeros(ChirpSNR.shape, dtype = np.int32)#;%No signal
    Chirpcombmask[L1B['Products/LowResPulse/Information/MaskCoPol' ][:]>2] = 1
    Chirpcombmask[L1B['Products/HighResPulse/Information/MaskCoPol'][:]>2] = 2
    SSL = 10 #Signal-to-Sidelobe Level
    Chirpcombmask[((10*np.log10(ChirpSNR)-10*np.log10(ChirpEstSL)) > SSL) & (L1B['Products/Chirp/Information/MaskCoPol'][:]>2)] = 3 #chirp
    Chirpcombmask[L1B['Products/Chirp/Information/MaskCoPol'][:]<=2] = 0
    
    print ('Adding Combined Data')
    ##Declare empty dicts
    L1B['Products/Combined'            ] = {}
    L1B['Products/Combined/Data'       ] = {}
    L1B['Products/Combined/Information'] = {} 
    
    ##Add initialize Combined/Data arrays
    L1B['Products/Combined/Data/dBZe'                 ] = L1B['Products/Chirp/Data/dBZe'                ]*np.nan
    L1B['Products/Combined/Data/Velocity_uncorrected' ] = L1B['Products/Chirp/Data/Velocity_uncorrected']*np.nan
    L1B['Products/Combined/Data/SpectrumWidth'        ] = L1B['Products/Chirp/Data/SpectrumWidth'       ]*np.nan
    L1B['Products/Combined/Information/SNR'           ] = L1B['Products/Chirp/Information/SNR'          ]*np.nan
        
    ##Add Products, by mask
    L1B['Products/Combined/Data/dBZe'                 ][Chirpcombmask==1] = L1B['Products/LowResPulse/Data/dBZe'                  ][Chirpcombmask==1]
    L1B['Products/Combined/Data/dBZe'                 ][Chirpcombmask==2] = L1B['Products/HighResPulse/Data/dBZe'                 ][Chirpcombmask==2]
    L1B['Products/Combined/Data/dBZe'                 ][Chirpcombmask==3] = L1B['Products/Chirp/Data/dBZe'                        ][Chirpcombmask==3]
    L1B['Products/Combined/Data/Velocity_uncorrected' ][Chirpcombmask==1] = L1B['Products/LowResPulse/Data/Velocity_uncorrected'  ][Chirpcombmask==1]
    L1B['Products/Combined/Data/Velocity_uncorrected' ][Chirpcombmask==2] = L1B['Products/HighResPulse/Data/Velocity_uncorrected' ][Chirpcombmask==2]
    L1B['Products/Combined/Data/Velocity_uncorrected' ][Chirpcombmask==3] = L1B['Products/Chirp/Data/Velocity_uncorrected'        ][Chirpcombmask==3]
    L1B['Products/Combined/Data/SpectrumWidth'        ][Chirpcombmask==1] = L1B['Products/LowResPulse/Data/SpectrumWidth'         ][Chirpcombmask==1]
    L1B['Products/Combined/Data/SpectrumWidth'        ][Chirpcombmask==2] = L1B['Products/HighResPulse/Data/SpectrumWidth'        ][Chirpcombmask==2]
    L1B['Products/Combined/Data/SpectrumWidth'        ][Chirpcombmask==3] = L1B['Products/Chirp/Data/SpectrumWidth'               ][Chirpcombmask==3]
    L1B['Products/Combined/Information/SNR'           ][Chirpcombmask==1] = L1B['Products/LowResPulse/Information/SNR'            ][Chirpcombmask==1]
    L1B['Products/Combined/Information/SNR'           ][Chirpcombmask==2] = L1B['Products/HighResPulse/Information/SNR'           ][Chirpcombmask==2]
    L1B['Products/Combined/Information/SNR'           ][Chirpcombmask==3] = L1B['Products/Chirp/Information/SNR'                  ][Chirpcombmask==3]    
    
    ##Add channel masks, plus units/descriptions
    L1B['Products/Combined/Information/ChannelMask'                     ] = Chirpcombmask
    L1B['Products/Combined/Information/ChannelMask_description'         ] = ['Composite Image Channel Mask. 0: No signal, 1: Low-Resolution Pulse, 2: High-Resolution Pulse, 3: Chirp']
    L1B['Products/Combined/Information/dBZe_units'                      ] = ['10*log10(mm^6/m^3)']
    L1B['Products/Combined/Information/dBZe_description'                ] = ['Equivalent reflectivity factor in dB with 1-sigma noise threshold applied. |K|^2 = ' + str(config['L1B_Kfactor'])]
    L1B['Products/Combined/Information/Velocity_uncorrected_units'      ] = ['m/s']
    L1B['Products/Combined/Information/Velocity_uncorrected_description'] = ['Doppler velocity with aircraft motion correction and 1-sigma noise threshold applied. Positive velocity is upward. Ka pitch offset has been removed.']
    L1B['Products/Combined/Information/SpectrumWidth_units'             ] = ['m/s']
    L1B['Products/Combined/Information/SpectrumWidth_description'       ] = ['Doppler velocity spectrum width estimate including aircraft motion. 1-sigma noise threshold applied.']    
    L1B['Products/Combined/Information/LDR_units'                       ] = ['dB'] ##we will add LDR data later
    L1B['Products/Combined/Information/LDR_description'                 ] = ['Linear Depolarization Ratio thresholded at 3-sigma, using high-res pulse for co-polarization and chirp for cross-polarization. They are pretty well matched, but dont get tripped up right next to the surface.']
    L1B['Products/Combined/Information/SNR_units'                       ] = ['Watt / Watt']
    L1B['Products/Combined/Information/SNR_description'                 ] = ['Estimated Signal to Noise Ratio'    ]

    #%Okay, here's the deal with kacross. The Ka-band cross-pol channel didn't have a loopback path. Therefore we first need to undo the calibration. (multiply by calenergy7). 
    #Then we need to normalize by gain and tx power. Assuming the radiometric signal is SIMILAR Z7(*calEnergy7) ./ CalEnergy6 * (Noise6./Noise7)
    #CalEnergy 6 ~ Pt G6...... Noise6/Noise7 ~ G6/G7  
    if radName == 'HIWRAP_KU':
        L1B['Products/Combined/Data/LDR'] = (L1A['Products/Channel3/Data/Zuncal'][:][times2save][:,rangeIndex[3]]/L1A['Products/Channel1/Data/Zuncal'][:][times2save][:,rangeIndex[1]]).astype('float32')
        L1B['Products/Combined/Data/LDR'] = 10*np.log10(L1B['Products/Combined/Data/LDR'])      

    if radName == 'HIWRAP_KA':      
        kacross = ((L1A['Products/Channel3/Data/Zuncal'][:].T*L1A['Products/Channel3/Information/CalEnergy'][:]/L1A['Products/Channel2/Information/CalEnergy'][:]*L1A['Products/Channel2/Information/EstimatedNoise'][:]/L1A['Products/Channel3/Information/EstimatedNoise'][:]).T).astype('float32')
        L1B['Products/Combined/Data/LDR'] = kacross[times2save][:,rangeIndex[3]]/L1A['Products/Channel1/Data/Zuncal'][:][times2save][:,rangeIndex[1]]
        L1B['Products/Combined/Data/LDR'] = 10*np.log10(L1B['Products/Combined/Data/LDR'])    

    L1B['Products/Combined/Data/LDR'][L1A['Products/Channel3/Information/SNR'][:][times2save][:,rangeIndex[3]]<3/np.sqrt(config['L1B_averagedPulses'])] = np.nan
    L1B['Products/Combined/Data/LDR'][L1A['Products/Channel1/Information/SNR'][:][times2save][:,rangeIndex[1]]<3/np.sqrt(config['L1B_averagedPulses'])] = np.nan

##Subroutine to perform NUBF correction
def nubf (dist, spd, ranges, beam, tempdbze, gradkern):
    from scipy import ndimage
    import numpy  as np
    
    ##Constant defined in NUBF paper by Loftus and McLinden
    alpha = lambda spd, rng, beamwidth: spd*0.0207*rng*np.deg2rad(beamwidth)**2
    
    ##Make a 2d array (tile) of speed so that it is a function of range.
    ##Use this for the alpha calculation.
    spd2d = np.tile(spd, (len(ranges), 1)).T
    
    ##Calculate the constant based on alpha equation and the spd2d
    constant = alpha(spd2d, ranges/1000., beam)
    
    ##Calculate the dBZ gradient in km
    ##Start with DY (distance)
    DY = np.convolve(dist, -gradkern, 'same')/1000 #m -> km
    
    ##Plane flies at 200m/s for 2 secs. Expect around 0.4km distance. Filter out outliers.
    DY[DY<.3] = np.nan
    DY[DY>.6] = np.nan
    
    ##Now calculate dBZe/DY. ndimage allows convolution on a 2d array.   
    grad = (ndimage.filters.convolve1d(tempdbze.T, -gradkern, axis=1, mode='constant')/DY).T

    ##The offset is the product of the alpha constant * gradient
    offset = grad * constant
    return (offset)
    
##Find radar-specific NUBF variables, then do the calculation
if radName == 'HIWRAP_KU' : 
    gradkern = np.array([-1, 0, 0, 0, 1]) #4 sample, 2 second kernel
    beam     = L1B['Products/Information/AntennaBeamwidth']
    tempdbze = np.copy(L1B['Products/Combined/Data/dBZe'  ]) ##KU is already filtered
    offset = nubf(L1B['Navigation/Data/NominalDistance'      ], spd, ranges, beam, tempdbze, gradkern)
    
    L1B['Products/Combined/Information/Velocity_nubf_offset'               ] = offset.astype(np.float32)
    L1B['Products/Combined/Information/Velocity_nubf_offset_units'         ] = ['m/s']
    L1B['Products/Combined/Information/Velocity_nubf_offset_description'   ] = ['The nonuniform beam filling (NUBF) offset used to correct Doppler velocity [Vel_corr = Vel_uncorr - nubf_offset - horizwind_offset]']

##Save the L1B dict to HDF file
print ('Writing dict')
dicttoh5(L1B, L1B_filepath)
h = h5py.File(L1B_filepath, 'a')

##Calculate the HRRR Winds along the flight path
print ('Processing winds')
L1B_uwind, L1B_vwind = L1B_sub_add_horiz_wind_hrrr_fullflight.L1B_sub_add_horiz_wind ('HRRR', radName, L1B_filepath, m_filepath)

##Calculate the along/crosstrack winds and return the winds.
##The L1B[Vel] variable is for plotting only.
print ('Reprocessing doppler')
if 'HIWRAP' in radName:
    cross, along, vel_corr_cross, vel_corr_along = \
        L1B_sub_along_cross_winds2.L1B_sub_along_cross_winds (radName, flightDate, L1B, L1B_uwind, L1B_vwind, L1B['Products/Combined/Data/Velocity_uncorrected'][:], cmap_n0q, deg_offset)
else:
    cross, along, vel_corr_cross, vel_corr_along = \
        L1B_sub_along_cross_winds2.L1B_sub_along_cross_winds (radName, flightDate, L1B, L1B_uwind, L1B_vwind, L1B['Products/Data/Velocity_uncorrected'][:], cmap_n0q, deg_offset)
##Reopen the newly created file and add these new variables

print ('Adding cross/along wind information')
h.create_dataset('/Products/Information/HRRR_CrossWind',             data=cross)
h.create_dataset('/Products/Information/HRRR_AlongWind',             data=along)
h.create_dataset('/Products/Information/HRRR_CrossWind_units',       data=np.string_(['m/s']))
h.create_dataset('/Products/Information/HRRR_AlongWind_units',       data=np.string_(['m/s']))
h.create_dataset('/Products/Information/HRRR_CrossWind_description', data=np.string_(['Cross wind from HRRR reanalysis, interpolated to this grid']))
h.create_dataset('/Products/Information/HRRR_AlongWind_description', data=np.string_(['Along wind from HRRR reanalysis, interpolated to this grid']))

if (radName == 'HIWRAP_KU'): ##Do this for both KU (with NUBF) and KA (no NUBF)
    h.create_dataset('/Products/Combined/Data/Velocity_corrected',                             data=(L1B['Products/Combined/Data/Velocity_uncorrected'] -L1B['Products/Combined/Information/Velocity_nubf_offset'] - vel_corr_cross -vel_corr_along).astype('float32')) #NUBF and winds
    h.create_dataset('/Products/Combined/Information/Velocity_corrected_units',                data=np.string_(['m/s']))
    h.create_dataset('/Products/Combined/Information/Velocity_corrected_description',          data=np.string_(['Doppler velocity corrected to account for intrusion of horizontal reanalysis winds and nonuniform beam filling (NUBF)']))
    h.create_dataset('/Products/Combined/Information/Velocity_horizwind_offset',               data=(vel_corr_cross +vel_corr_along).astype('float32'))
    h.create_dataset('/Products/Combined/Information/Velocity_horizwind_offset_units',         data=np.string_(['m/s']))
    h.create_dataset('/Products/Combined/Information/Velocity_horizwind_offset_description',   data=np.string_(['The horizontal wind offset used to correct Doppler velocity [Vel_corr = Vel_uncorr - nubf_offset - horizwind_offset]']))
if (radName == 'HIWRAP_KA'): ##Do this for both KU (with NUBF) and KA (no NUBF)
    h.create_dataset('/Products/Combined/Data/Velocity_corrected',                             data=(L1B['Products/Combined/Data/Velocity_uncorrected'] - vel_corr_cross -vel_corr_along).astype('float32'))
    h.create_dataset('/Products/Combined/Information/Velocity_corrected_units',                data=np.string_(['m/s']))
    h.create_dataset('/Products/Combined/Information/Velocity_corrected_description',          data=np.string_(['Doppler velocity corrected to account for intrusion of horizontal reanalysis winds']))
    h.create_dataset('/Products/Combined/Information/Velocity_horizwind_offset',               data=(vel_corr_cross +vel_corr_along).astype('float32'))
    h.create_dataset('/Products/Combined/Information/Velocity_horizwind_offset_units',         data=np.string_(['m/s']))
    h.create_dataset('/Products/Combined/Information/Velocity_horizwind_offset_description',   data=np.string_(['The horizontal wind offset used to correct Doppler velocity [Vel_corr = Vel_uncorr - horizwind_offset]']))



if ((makePlots == 'yes') | (makePlots == 'y')):
    print ('Plotting')
    if 'HIWRAP' in radName:
        L1B_sub_plotter.L1B_sub_plotter (radName,       h, flightLines, L1B_file, L1B_plotpath, cmap_n0q, flightDate, versionLetter)
h.close()