#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mar 21 15:57:50 2020

@author: JKH

This function helps with matching up EM27/SUN measurement times with GPS and 
meteorological measurement times. It replaces the core of egi_create_in_file.py 
and em27_fix_sunrun.py from EGI2014.
I'm making the intermediate file instead of just a direct *.in file for I2S so 
it can work for both interferograms, or spectra.


It gets called in by matchmet and does the following:
1) Reads in list of file sizes, flags which ones to further processes
2) Reads in *_sigtime.dat file and calculates arbitrary signal
3) Based on signal, flags which files are likely cloudy or too low of signal to 
    be much use anyways.
4) An output file of just ifgs with sufficient signal is output.
5) Of the "uncloudy" files, met & location are matched up
6) A prediction of the EM27 configuration is also made (dual, single, 
 single-extended, or mid). Except for single/single-extended the type can
 be inferred from an OPUS header.
7) An output file is written with the instrument type and the matched-up data
    plus flag values.

Inputs:
fname_pre: full file path to directory where header example, *_sigtime.dat, 
    and met and gps files are for the day, also includes the xxYYYYMMDD prefix
location_file: path to the file_loc.lst

"""


############ import some sub-functions
import numpy as np
import os
import pandas as pd
#import copy

curr_dir = os.getcwd()
os.chdir('../python_subs')
import read_hdr
import egi_subs as egs
import readusbgpsfiles as rgps
import read_timesig_fb as rtime
import findclosest_n as fclsn
import readzenofiles as rzeno
os.chdir(curr_dir)





def metmatch_sub(fname_pre, location_file):
    ### example inputs
#    fname_pre = '/home/osboxes/research/tccon/egi2020/results/daily/20140613/matchup/cn20140613'
#    fname_pre = '/home/osboxes/research/tccon/egi2020/results/xc/daily/20180716/matchup/xc20180716'
#    location_file = [34.14,-118.13,237]
    
    
    
    ############# read in header file to get the UTC offset #############
    hdr_raw = read_hdr.read_hdr(fname_pre + '_hdr1.dat')
    if 'TIM1' in hdr_raw:
      ifg_utc_off = float(hdr_raw['TIM1'].split('GMT')[1].split(')')[0]) #E.g., Harvard=-4
    else: #likely I2S files, NOTE: changed "GMT" to "UTC" as well
      ifg_utc_off = float(hdr_raw['TIM'].split('UTC')[1].split(')')[0]) #E.g., Harvard=-4

    
    ############## set approximate longitude for approximate timezone #############
    long_in = -118.1 #default to Caltech...
    if not(location_file == 0 or egs.is_number(location_file[0])): # first character not a number, filepath given
        gps_loc_raw = rgps.readusbgpsfiles(fname_pre + '_cmb_gps.gps')
        if 'CompSrlDate' in gps_loc_raw:
            gps_time = gps_loc_raw['CompSrlDate'] #serial date (matlab version) exists.  Same time as ifgs. 
        elif gps_loc_raw['utc_lcl'] == 2: #local (or rather computer) time flag
            gps_time_p = egs.dt_from6vec(gps_loc_raw['Cyear'],gps_loc_raw['Cmonth'],gps_loc_raw['Cday'],gps_loc_raw['Chour'],gps_loc_raw['Cmin'],gps_loc_raw['Csec']) #matches spectrometer/computer time
            gps_time = gps_time_p #+ 719529 #converts to matlab type
        elif gps_loc_raw['utc_lcl'] == 1: #utc time
            gps_time_p = egs.dt_from6vec(gps_loc_raw['Uyear'],gps_loc_raw['Umonth'],gps_loc_raw['Uday'],gps_loc_raw['Uhour'],gps_loc_raw['Umin'],gps_loc_raw['Usec']) + ifg_utc_off/24 #matches spectrometer/computer time
            gps_time = gps_time_p #+ 719529 #converts to matlab type
        
        long_in = np.nanmedian(gps_loc_raw['Long'])
                 
    
    elif len(location_file) == 3: #lat, long, alt_m given
        long_in = location_file[1]
    
    
    
    ############ read in time/sig file #############
    time_sig = rtime.read_timesig_fb(fname_pre + '_sigtime.dat', long_in, -1*ifg_utc_off)
    n_intfgs = len(time_sig['intfg_name'])
    filedates = time_sig['raw_time'] # + 719529 #converts to matlab type 
    
    
    ########### Read in file sizes #################
    # This is used as a diagnostic
    fsize_dat = pd.read_csv(fname_pre + '_fsizes.dat', names=['fname', 'bytes'], delimiter=',')
    fsize_flag = np.ones(len(fsize_dat['fname'])) #initially ALL are flagged
    fname_all = [None]*len(fsize_dat['fname']) #includes bad/small files
    for i1, iFile in enumerate(fsize_dat['fname']):
        fname_all[i1] = os.path.basename(iFile)
        if fname_all[i1] in time_sig['intfg_name']:
            fsize_flag[i1] = 0
    
    
    ############# replace location data with correct location ############
    # set Default as Caltech (Should this even be here anymore??? I guess users could put their own "home station" in...)
    gps_avg_corr_data = np.outer( np.ones([n_intfgs,1]), [34.13625, -118.12697, 237])
    if not(location_file == 0 or egs.is_number(location_file[0])): # first character not a number, filepath given
        gps_closedate_p = fclsn.findclosest_n(filedates, gps_time, 3) #up to 3 points
        gps_closedate = gps_closedate_p[0]
        gpsind_win_min = abs(gps_closedate['closest_ind_diff']) <= float(1)/(24*60) # %indexes where data was taken within one minute
        gps_avg_corr_data = -99.9*np.ones([len(filedates),3]) #pre-allocating to see what's been written
        
        
        if 'masl' in gps_loc_raw:
            alt_all = gps_loc_raw['masl'] 
        else:
            alt_all = gps_loc_raw['masl_inacc'] + gps_loc_raw['geoidal_WGS84_sep'] + 18; #the +18 seems like it would make the measurement more accurate (from an empirical observation at Caltech)
        
        i1 = 0
        for fdate_i1 in filedates:
            gps_data_indx_i1_p2 = gps_closedate['closest_ind'][i1, :]
            gps_data_indx_i1_p1 = gps_data_indx_i1_p2[gpsind_win_min[i1,:]]
            gps_data_indx_i1 = gps_data_indx_i1_p1.astype(int)
            
            ## Average up to 3 points within +/- 1 minute             
            gps_avg_corr_data[i1,0] = egs.nanmean( gps_loc_raw['Lat'][gps_data_indx_i1])
            gps_avg_corr_data[i1,1] = egs.nanmean(gps_loc_raw['Long'][gps_data_indx_i1])
            gps_avg_corr_data[i1,2] = egs.nanmean(alt_all[gps_data_indx_i1])
            
            i1 += 1 #iterates by 1
         
        gps_avg_corr_data[np.isnan(gps_avg_corr_data)] = -999 #converts NaNs to -999
        
        
    elif len(location_file) == 3: #lat, long, alt_m given
        # this is probably the most common option... 
        gps_avg_corr_data = np.outer( np.ones([n_intfgs,1]), location_file)
    
#    elif location_file == 0: #Default to Caltech, already formatted
#        gps_avg_corr_data = np.outer( np.ones([n_intfgs,1]), [34.13625, -118.12697, 237])
#        
#    else: #all others default to Caltech
#        gps_avg_corr_data = np.outer( np.ones([n_intfgs,1]), [34.13625, -118.12697, 237])
    
    
    
    ########### Match up the met data #############
    ## Zeno (or rather meteorological) data read in and average and coordinate with times interferograms were taken
    
    zenofraw = rzeno.readzenofiles(fname_pre + '_cmb_met.zno')
    
    if 'CompSrlDate' in zenofraw:
        zeno_date = zenofraw['CompSrlDate'] #serial date (matlab version) exists.  Same time as ifgs. 
    elif zenofraw['utc_lcl'] == 2: #local time
        zeno_date_p = egs.dt_from6vec(zenofraw['Cyear'],zenofraw['Cmonth'],zenofraw['Cday'],zenofraw['Chour'],zenofraw['Cmin'],zenofraw['Csec']) #matches spectrometer/computer time
        zeno_date = zeno_date_p #+ 719529 #converts to matlab type
    elif zenofraw['utc_lcl'] == 1: #utc time
        zeno_date_p = egs.dt_from6vec(zenofraw['Uyear'],zenofraw['Umonth'],zenofraw['Uday'],zenofraw['Uhour'],zenofraw['Umin'],zenofraw['Usec']) + ifg_utc_off/24 #matches spectrometer/computer time
        zeno_date = zeno_date_p #+ 719529 #converts to matlab type
        
    
        
    # add in -9 if no WSPD and WDIR data
    if not('WSPD' in zenofraw):
        zenofraw['WSPD'] = np.ones([len(zeno_date),1])*0
    
    if not('WDIR' in zenofraw):
        zenofraw['WDIR'] = np.ones([len(zeno_date),1])*0
    
    if not('RH' in zenofraw): #if no RH data in the met files
        zenofraw['RH'] = np.ones([len(zeno_date),1])*0
    
    zeno_closedate_p = fclsn.findclosest_n(filedates, zeno_date, 5) #up to 5 points
    zeno_closedate = zeno_closedate_p[0]
    zenoind_win_2min = abs(zeno_closedate['closest_ind_diff']) <= float(2)/(24*60) # %indexes where data was taken within two minutes. Needs float or Python 2 rounds.
    
    zeno_avg_corr_data = -99.9*np.ones([len(filedates),5]) #data from zeno that goes to input file includes temp, pres, RH, winddir, windspd - pre-allocating to see what's been written
    
    if 'Temp' in zenofraw:
        zeno_temp = zenofraw['Temp']
    elif 'Tout' in zenofraw:
        zeno_temp = zenofraw['Tout']
    else: #no external Temperature data
        zeno_temp = np.ones([len(zeno_date),1])*-99.9
    
    
    i1 = 0
    for fdate_i1 in filedates: #over each ifg/timepoint
        zeno_data_indx_i1_p2 = zeno_closedate['closest_ind'][i1, :]
        zeno_data_indx_i1_p1 = zeno_data_indx_i1_p2[zenoind_win_2min[i1,:]]
        zeno_data_indx_i1 = zeno_data_indx_i1_p1.astype(int)
         
        zeno_avg_corr_data[i1,0] = egs.nanmean(zeno_temp[zeno_data_indx_i1])
        zeno_avg_corr_data[i1,1] = egs.nanmean(zenofraw['RH'][zeno_data_indx_i1])
        zeno_avg_corr_data[i1,2] = egs.nanmean(zenofraw['Pout'][zeno_data_indx_i1])
        wdirspd_avg = egs.avg_wdir(zenofraw['WDIR'][zeno_data_indx_i1], zenofraw['WSPD'][zeno_data_indx_i1])
        zeno_avg_corr_data[i1,3] = wdirspd_avg[1] #
        zeno_avg_corr_data[i1,4] = wdirspd_avg[0]
        
        i1 += 1 #iterates by 1
    
    ### Round to a more reasonable number of digits
    zeno_avg_corr_data[:,0] = np.round(zeno_avg_corr_data[:,0],1) #Tout
    zeno_avg_corr_data[:,1] = np.round(zeno_avg_corr_data[:,1]) #RH
    zeno_avg_corr_data[:,2] = np.round(zeno_avg_corr_data[:,2],2) #Pins
    zeno_avg_corr_data[:,3] = np.round(zeno_avg_corr_data[:,3],1) #WSPD
    zeno_avg_corr_data[:,4] = np.round(zeno_avg_corr_data[:,4]) #WDIR
    
    
    n_a = np.isnan(zeno_avg_corr_data)
    zeno_avg_corr_data[n_a] = -999 #switch NaNs to -999
    
    ######### Get vectors of time #########
    time_utc_mat = egs.dt_to6vec(time_sig['UTC_time'])
    yyyymmdd_vec = time_utc_mat[:,0]*10000 + time_utc_mat[:,1]*100 + time_utc_mat[:,2]
    HHMMSSm_vec = time_utc_mat[:,3]*100000 + time_utc_mat[:,4]*1000 + np.round(time_utc_mat[:,5]*10)
    
    
    ######### Flag for if filtered because of clouds ############
    # flag of 1=bad file size, flag of 2=cloudy
    # Files that have already been skipped due to size are included at all
    init_flag = np.zeros([len(time_sig['intfg_name'])])
    init_flag[time_sig['arb_sia']<25] = 2 #25 for a cutoff is conservative, meant to pass most ifgs
    
    
    ####### Put fields of interest into data frame, prep for file writing #######
    part_df = pd.DataFrame.from_dict( {'Iflag': init_flag,
        'filename':  time_sig['intfg_name'], 
        'fANDb': np.sum([time_sig['GBW']>0, time_sig['GFW']>0], 0),
        'UYYYYMMDD': yyyymmdd_vec, 
        'UHHMMSSm': HHMMSSm_vec, 
        'lat': gps_avg_corr_data[:,0], 
        'lon': gps_avg_corr_data[:,1], 
        'alt': gps_avg_corr_data[:,2], 
        'sia': np.round(time_sig['arb_sia'], 1), 
        'Tout': zeno_avg_corr_data[:,0], 
        'RH':  zeno_avg_corr_data[:,1], 
        'Pout': zeno_avg_corr_data[:,2], 
        'WSPD': zeno_avg_corr_data[:,3], 
        'WDIR': zeno_avg_corr_data[:,4], 
        'Tins': time_sig['TSC']} )
    
    
    ### Add in extra (empty) rows due to small file sizes
    row_fmt = [np.int64(10), -9, np.int64(0), np.int64(-9999999), np.int64(-999999), -99.99, -99.99, -99.99, 
                  -99.99, -99.99, np.int64(-99), -99.99, -99.99, np.int64(-99), -99.99]
    out_df = pd.DataFrame(data=row_fmt*np.ones([len(fsize_flag), 1 ]), columns=part_df.columns )
    for iCol in part_df.columns:
        out_df.loc[fsize_flag==0, iCol] = part_df.loc[:, iCol]
    
    
    out_df.loc[fsize_flag==1, 'Iflag'] = 1
    out_df.loc[:, 'filename'] = fname_all
    
    ### Convert some fields to integers
    for iCol in ['Iflag', 'fANDb', 'RH', 'WDIR', 'UYYYYMMDD', 'UHHMMSSm']:
        out_df.loc[:, iCol] = np.int64(out_df.loc[:, iCol])
    
    
    
    ######## Make a prediction of the EM27 configuration (dual, single, 
    # mid) from the OPUS header. The single/single-extended distinction is 
    # not made here (can't be determined solely from header).
    expected_spec = np.sum(out_df.loc[out_df.loc[:, 'Iflag']==0, 'fANDb']) #count the number of expected spectra
    
    if hdr_raw['INS'] == 'EM27/SUN MIR': #seems like rarest
        inst_config = 'em27_mid'
    elif 'NPT2' in hdr_raw and np.int64(hdr_raw['NPT2'])>0: 
        # The presence of a second channel likely means it's in a dual configuration
        inst_config = 'em27_dual'
        expected_spec = expected_spec*2 #double the number of spectra ('a' and 'd')
    else: #everything else...
        inst_config = 'em27_single'
#    elif hdr_raw['INS'] == 'EM27SUN' or hdr_raw['INS'] == 'EM27/SUN':
#        inst_config = 'em27_single'
#    else:
#        inst_config = 'unknown' 
    
    
    #########  Get median of each of the location bins for day
    # (Most of the time this should just be one location, but on occasion the EM27
    # may have made measurements in different 1 deg grid cells)
    ll_unflagged = np.array([out_df.loc[out_df.loc[:, 'Iflag']==0, 'lat'], out_df.loc[out_df.loc[:, 'Iflag']==0, 'lon']])  #2 x n array
    _, ll_binidx = np.unique(np.round(ll_unflagged), return_inverse=True, axis=1) #unique rows 
    mod_lls = [] #lats and lon for the modtran input
    for i1 in range(np.max(ll_binidx)+1): #plus 1 because of python end is exclusive
      mod_lls.append(np.array([np.median(ll_unflagged[0, ll_binidx==i1]), 
               np.median(ll_unflagged[1, ll_binidx==i1])]))
    
    #### Add in days before / after if it appears measurements may cross over into
    # previous or next UTC day. The ephem toolbox would be more accurate, but is not
    # natively included in Anaconda Python
#    utc_noon = np.array(-1*out_df.loc[out_df.loc[:, 'Iflag']==0, 'lon']/15 + 12) #approximate UTC solar noon
#    # See: https://en.wikipedia.org/wiki/Position_of_the_Sun
#    day_in_year = (time_sig['UTC_time'] - 730486) % 365.24 #approximate date in year (730486 is for Jan 1, 2000)
#    solar_decl = -1*np.arcsin(0.39779 * np.cos(np.pi/180*(0.98565*(day_in_year+10) + 1.914*np.sin(0.98565*np.pi/180*(day_in_year-2))  ))  )
#    hour_angle_rad = np.arccos(-1*np.tan(np.pi/180*out_df.loc[out_df.loc[:, 'Iflag']==0, 'lat'])*np.squeeze(np.tan(solar_decl)) )
#    hour_angle_hr = np.array(hour_angle_rad/(2*np.pi)*24)
#    utc_srise = utc_noon - hour_angle_hr - 4 #padding by 4 for up to an hour error in calculation, and 3 hours between *.mod files
#    utc_sset = utc_noon + hour_angle_hr + 4
    
    utc_dt_day = egs.dt_from6vec(np.float64(fname_pre[-8:-4]), np.float64(fname_pre[-4:-2]), np.float64(fname_pre[-2:]), 0, 0, 0)
    utc_dt_tod = np.squeeze(time_sig['UTC_time'] - utc_dt_day)*24 #hour
    
    days_add = [0]
    if np.any(utc_dt_tod<3.5): #within 3.5 hours of beginning of day
      days_add.append(-1)
    if np.any(utc_dt_tod>20.5): #within 3.5 hours of end of day
      days_add.append(1)
    
    
    
    ######### Actual file writing
    with open(fname_pre + '_mmtch.dat', 'w') as fid:
      fid.write('Inst_config: ' + inst_config + '\n')
      fid.write('nRaw: ' + '%d' % np.sum(out_df.loc[:, 'Iflag']==0) + 
                  ' nSpc: ' + '%d' % expected_spec + '\n')
      for iAdd in days_add: #loop accounts for +/-1 UTC days
        vec_day = egs.dt_to6vec(utc_dt_day + iAdd)
        
        num_day = '%4d%02d%02d' % (vec_day[0,0], vec_day[0,1], vec_day[0,2]) #convert back to a string format YYYYMMDD
        for iLoc in mod_lls: #loop accounts for different grids in same day
                fid.write(fname_pre[-10:-8] + ' ' + '%s' % num_day + 
                      ' %7.3f' % iLoc[0] + 
                      ' %8.3f' % iLoc[1] + '\n')   
        
      fid.write('\n') #spacer row
      out_df.to_csv(fid, index=False, float_format='%.4f', sep=' ') #4 decimals on float esp. for lat/lon
    
    
    ######## Write a simplar *_opus.gnd with just the unfiltered files
    # This would be an input to create_sunrun
    with open(fname_pre + '_opus.gnd', 'w') as fid:
        out_df.loc[out_df.loc[:, 'Iflag']==0, 'filename'].to_csv(fid, index=False, header=False) #
    
    
    out = fname_pre + '_mmtch.dat'
    return out
    
    
    
    
    
    
