# -*- coding: utf-8 -*-
"""
Created on Mon Aug 17 10:16:43 2015

@author: Jacob Hedelius
"""

#!/bin/python
# Filename egi_create_in_file.py

#__author__ = 'Jacob Hedelius'




## import some sub-functions
import numpy as np
import pandas as pd
import os
curr_dir = os.getcwd()
os.chdir('../python_subs')
import read_top_params as rtop
#    from read_top_params import *
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 egi_create_in_file2020(fname_pre):
    """
    Creates *.in file for EGI 2020. I've modified things so that the met and 
    location matching is done before in metmatch_sub.py.
    I'm assuming filenames based on an input path name.
    All that's left to do is just format the output and concatenate with the *.top
    file.
    
    example input
    fname_pre = '/home/osboxes/research/tccon/egi2020/results/xd/daily/20180716/matchup/xd20180716'
    """
    
    top_fname = fname_pre + '-i2s.top'
    mm_fname = fname_pre + '_mmtch.dat'
    output_fname = fname_pre + '-i2s.in'
    
    #### Read in the met match file
#    fobj = open(mm_fname, 'r') # gets file object for reading
#    file_lines = fobj.readlines() #reads in entire file all at once, line by line
#    fobj.close()
    
    # Find the column header line by looking for line beginning with "Iflag"
    with open(mm_fname, 'r') as fobj: # gets file object for reading
      mm_raw = fobj.readlines() 
    n_mm_hdr = next(i1 for i1, val in enumerate(mm_raw) if val[0:5] == 'Iflag')-1
    
    mm_df = pd.read_csv(mm_fname, delim_whitespace=True, header=n_mm_hdr) 
    not_flag = mm_df['Iflag'] == 0  #only unflagged files get processed
    mm_noflag = mm_df.loc[not_flag,:] #get rid of flagged rows
    mm_noflag = mm_noflag.reset_index(drop=True)
    
    
    ### Format output strings 
    yr_list1 = np.floor(mm_noflag['UYYYYMMDD']/10000).astype(np.int32).tolist()
    yr_list = [str(l_val).zfill(4) for l_val in yr_list1] #no need to add leading zeros on year, but converts to stings
    
    mth_list1 = np.floor((mm_noflag['UYYYYMMDD'] % 10000)/100 ).astype(np.int32).tolist()
    mth_list = [str(l_val).zfill(2) for l_val in mth_list1] #adds a leading zero  if needed
    
    day_list1 = np.floor((mm_noflag['UYYYYMMDD'] % 100) ).astype(np.int32).tolist()
    day_list = [str(l_val).zfill(2) for l_val in day_list1] #adds a leading zero if needed
    
    run_num = np.arange(1, np.sum(not_flag)*2, 2) #run numbers as odd
    rn_list1 = run_num.tolist() #converting the array to a list
    rn_list = [str(l_val).zfill(4) for l_val in rn_list1] #adds leading zeros as necessary
    
    
    lat_list =  ["%8.4f" % x for x in  mm_noflag['lat'] ]
    lon_list = ["%9.4f" % x for x in mm_noflag['lon'] ]
    alt_list = ["%4.f" % x for x in mm_noflag['alt'] ]
    tins_list = ["%4.1f" % x for x in mm_noflag['Tins'] ] 
    pins_list = ["%7.2f" % x for x in mm_noflag['Pout'] ]
    hins_list = ["%5.1f" % x for x in mm_noflag['RH'] ]
    tout_list = ["%4.1f" % x for x in mm_noflag['Tout'] ]
    pout_list = pins_list
    hout_list = hins_list
    sia_list = ["%6.1f" % x for x in mm_noflag['sia'] ]
    
    ## Put in spacers for FVSI
    n9999_null_p = np.ones([np.sum(not_flag),1])*-0.9999
    fvsi_list = ["%6.4f" % x for x in n9999_null_p]
    
    wspd_list = ["%5.1f" % x for x in mm_noflag['WSPD'] ]
    wdir_list = ["%5.1f" % x for x in mm_noflag['WDIR'] ]
    
    
    outtext = [mm_noflag['filename'], yr_list, mth_list, day_list, rn_list, lat_list, 
                lon_list, alt_list, tins_list, pins_list, hins_list, tout_list, 
                pout_list, hout_list, sia_list, fvsi_list, wspd_list, wdir_list]
    
        
    ## Write to output file
    with open(output_fname, 'w') as fileID:
        ## Write in the *.top part
        with open(top_fname) as infile:
            fileID.write(infile.read())
        
        ## Write in the bottom part
        for iLine in range(len(yr_list)): 
            line_txt = outtext[0][iLine]
            for iCol in range(1, len(outtext)):
                line_txt = line_txt + ' ' + outtext[iCol][iLine]
                
            fileID.write(line_txt + '\n')
    
    
    return output_fname








def egi_create_in_file(filenametop, intfg_timesigfile, zenofilename, location_file):
    """
    This script is meant to create i2s files.  It combines the correct file names
    # and adds in weather and location data where necessary.
    #
    # Pre-operations needed include: (1) get a Zeno weather file
    # (2) Extract time headers using the read_timesig_fb script.  
    # (3) a formatted input .top file
    # (4) (optional) get a file for location (GPS) note that this is not
    # necessary if at a single location.  Enter in a "0" for Caltech, or enter
    # in longitude and latitude and altitude (m) (ex. [-118.1, 34.2, 237.] for 118.1 west and 34.2
    # north)
    # These pre-operations call all be performed by running the shell script
    # i2sprep
    
    # NOTE: This function is based on a similar function to fix sunruns.
    # Variable names for the most part have been changed...
    
    ## Test values and code
    #    filenametop = 'C:\\Users\\Jacob Hedelius\\Documents\\Research\\python\\example_data\\cn20140613opus-ipp.top'
    #    intfg_timesigfile = 'C:\\Users\\Jacob Hedelius\\Documents\\Research\\python\\example_data\\cnintf_sigtime20140602.dat' 
    #    zenofilename = 'C:\\Users\\Jacob Hedelius\\Documents\\Research\\python\\example_data\\cn_cmb_zeno.zno'
    #    location_file = [-118.1, 34.2, 237]
    # version 2017-Nov-11
    """
    
    
    ## setup file locations
    files_dir_ind_l = max(intfg_timesigfile.rfind('/'), intfg_timesigfile.rfind('\\')) + 1 #last occurance
    files_dir = intfg_timesigfile[0:files_dir_ind_l] #directory for where .in files will go (where time_sig files came from)
    
    files_dir_ind_l2 = max(intfg_timesigfile.rfind('/', 0, files_dir_ind_l), intfg_timesigfile.rfind('\\', 0, files_dir_ind_l)) #second to last
    files_dir_ind_l3 = max(intfg_timesigfile.rfind('/', 0, files_dir_ind_l2), intfg_timesigfile.rfind('\\', 0, files_dir_ind_l2))+1 #3rd to last
    
    day_folder = intfg_timesigfile[files_dir_ind_l3:files_dir_ind_l2] #typically in egi the day in yyyymmdd (formerly named "month_folder" for some reason...)
    siteid = intfg_timesigfile[(files_dir_ind_l):(files_dir_ind_l+2)]; #2 letter site identifier
    
    filenameout_botm_list = files_dir + siteid + day_folder + 'i2s.botm'
    infilenameout = files_dir + siteid + day_folder + 'opusi2s.in' #name of .in file
    SIAfilenameout = files_dir + siteid + day_folder + '_sias.dat' #list of (arbitrary) SIA values
    
    
    ## read in .top file (for time correction, already written in by egiprepdaily)
    top_params = rtop.read_top_params(filenametop)
    time_fix2utc = top_params['add_get_utc_19']
    
    
    ## set longitude
    if not(location_file == 0 or egs.is_number(location_file[0])): # first character not a number, filepath given
        gps_loc_raw = rgps.readusbgpsfiles(location_file)
        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 time
            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']) - time_fix2utc/24 #matches spectrometer/computer time
            gps_time = gps_time_p #+ 719529 #converts to matlab type
        
        
        long_in = np.nanmedian(gps_loc_raw['Long'])
#        long_in = gps_loc_raw.Long[0] #manually sets first longitude
        
    elif location_file == 0: #Default to Caltech, already formatted
        long_in = -118.1 
    
    elif len(location_file) == 3: #lat, long, alt_m given
        long_in = location_file[1]
    
    else: #all others default to Caltech
        long_in = -118.1 
    
    ## read in time/sig file
    time_sig = rtime.read_timesig_fb(intfg_timesigfile, long_in, time_fix2utc)
    n_intfgs = len(time_sig['intfg_name'])
    
    ## replace location data with correct location
    if not(location_file == 0 or egs.is_number(location_file[0])): # first character not a number, filepath given
        pass        
        #gps_avg_corr_data[0,1] = gps_loc_raw.Long[0] # may be unused
        
    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])
    
    elif len(location_file) == 3: #lat, long, alt_m given
        gps_avg_corr_data = np.outer( np.ones([n_intfgs,1]), location_file)
    
    else: #all others default to Caltech
        gps_avg_corr_data = np.outer( np.ones([n_intfgs,1]), [34.13625, -118.12697, 237])
    
    ## finish GPS cooridnation if needed
    filedates = time_sig['raw_time'] # + 719529 #converts to matlab type 
    
    if not(location_file == 0 or egs.is_number(location_file[0])): 
        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
        
        # gps_avg_corr_data_orig = gps_avg_corr_data; %original values, without nans - UNUSED
        
        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)
            # alt_i1 = gps_loc_raw.data(gps_data_indx_i1, 8) + gps_loc_raw.data(gps_data_indx_i1, 9);
             
            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
    
    
    ## read in arbitrary solar signals
    
    outtext = time_sig['outtext']
    arb_sia1 = np.round(time_sig['arb_sia'],1) #gives 1 decimal
    arb_sia2 = ["%6.1f" % x for x in arb_sia1] #converts to list of 6.1 floating format of strings
#    intfg_names = time_sig['intfg_name']
    
    
    
    ## replace GPS
    lat_out = ["%8.4f" % x for x in gps_avg_corr_data[:,0] ]
    long_out = ["%9.4f" % x for x in gps_avg_corr_data[:,1] ]
    alt_out = ["%4.f" % x for x in gps_avg_corr_data[:,2] ]
    
    outtext = outtext + [lat_out, long_out, alt_out]
    
    
    ## Zeno data read in and average and coordinate with times interferograms were taken
    
    zenofraw = rzeno.readzenofiles(zenofilename)
    
    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']) - time_fix2utc/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 Temperature data
        zeno_temp = np.ones([len(zeno_date),1])*0
    
    
    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
    
    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],1) #Pout
    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] = -99.9
    
    
    ## TODO: This is where the internal EM27 data on pressure, humidity, etc. could be read in
    #NONE IS AVALIABLE
    
    ## Replace false weather data with Zeno data
    tins_list1 = time_sig['TSC'].tolist()
    tins_list = ["%4.1f" % x for x in tins_list1] #format list properly
    
    Tout_v = zeno_avg_corr_data[:,0] #tout original 
    tout_list1 = Tout_v.tolist()
    tout_list = ["%4.1f" % x for x in tout_list1]
    
    RH_v = zeno_avg_corr_data[:,1]
    rh_list1 = RH_v.tolist()
    rh_list = ["%5.1f" % x for x in rh_list1]
    
    Pout_v = zeno_avg_corr_data[:,2]
    pout_list1 = Pout_v.tolist()
    pout_list = ["%6.1f" % x for x in pout_list1]
    
    WSPD_v = zeno_avg_corr_data[:,3]
    wspd_list1 = WSPD_v.tolist()
    wspd_list = ["%5.1f" % x for x in wspd_list1]
    
    wdir_v = zeno_avg_corr_data[:,4]
    wdir_list1 = wdir_v.tolist()
    wdir_list = ["%5.1f" % x for x in wdir_list1]
    
    outtext = outtext + [tins_list, pout_list, rh_list, tout_list, pout_list, rh_list]
    
    ## Put in spacers for SIA and FVSI
    n9999_null_p = np.ones([n_intfgs,1])*-0.9999
    n9999_null = ["%6.4f" % x for x in n9999_null_p]
    outtext = outtext + [n9999_null, n9999_null]
    # outtext = outtext + [arb_sia2, n9999_null] #use this line instead of above line for GGG201x
    # note: changed SIA null from -1.0 to -.9999 and %4.1f to %6.4f
    #when both these values are null (-.9999) then the variation in DC is put in as FVSI.  Later on (if using egi) and arbitrary SIA is put in.
    
    
    ## Save an SIA file seperately
    fileID = open(SIAfilenameout, 'w')
    i1 = 0
    for sia_i1 in arb_sia2: 
        fileID.write(sia_i1 + '\n')
        
        if (time_sig['GFW'][i1] >= 1 and time_sig['GBW'][i1] >= 1):
            fileID.write(sia_i1 + '\n') #print row twice if forward and reverse pair exists
        
        i1 += 1 #iterates by 1
    
    fileID.close()
    
    ## Windspd and Winddir
    outtext = outtext + [wspd_list, wdir_list]
    
    
    ## Actual file writing
    fileID = open(filenameout_botm_list, 'w') #open file for writing
    i1 = 0
    out_cols = len(outtext)
    for ifg_i1 in outtext[0]: #over each ifg name
        i2 = 1
        i1_line = ifg_i1
        while i2<out_cols:
            i1_line = i1_line + ' ' + outtext[i2][i1]
            i2 += 1
        
        fileID.write(i1_line + '\n')
        
        i1 += 1 #iterates by 1
    
    fileID.close()
    
    
    ## catenate files
    cat_filenames = [filenametop, filenameout_botm_list]
    with open(infilenameout, 'w') as outfile:
        for fname in cat_filenames:
            with open(fname) as infile:
                outfile.write(infile.read())
    
    out = {'filename': filenameout_botm_list}
    return out



