#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on  Apr 13 2020

@author: JKH

Contains some plotting tools to work with GGG outputs in EGI standard directory
format and naming conventions.

"""
import glob
import os
import numpy as np
import re
import copy
import pandas as pd
import datetime as dt
import requests
import shutil
from itertools import compress

from pandas.plotting import register_matplotlib_converters #suppresses a future warning
register_matplotlib_converters() #suppresses a future warning

from matplotlib.backends.backend_pdf import PdfPages #for multi-pdf pages
import matplotlib.pyplot as plt 
import matplotlib.image as mpimg
import matplotlib
import matplotlib.dates as mdates
import matplotlib.path as mpath





curr_dir = os.getcwd()
os.chdir('../python_subs')
import read_ggg_misc as rgggm
import egi_subs as egs
os.chdir(curr_dir)


#%%
def egi_make_plots(oof_ffname, flag_list=['flag', 'eq', 0]):
  """
  This main function calls in the sub-functions that make the plots
  Add additional lines as other plotting/*.pdf creation routines are added.
  
  Example input:
  oof_ffname = '/home/osboxes/research/tccon/egi2020/results/cn_spec/daily/20140609/cn20140609.vav.ada.aia.oof'
  oof_ffname = '/home/osboxes/research/tccon/egi2020/results/xa/daily/20140603/xa20140603.vav.ada.aia.oof'
  """
  plot_spt_samples(oof_ffname)
  make_summary_pdf(oof_ffname, flag_list)
  plot_misc_tseries(oof_ffname, flag_list)
  
  # Other ideas for plots or subplots:
  # add 5-minute means to timeseries?
  # Something with AOD??
  # Allan variance?
  # Quantile plots?
  # A priori profiles?
  # HYSPLIT trajectories? (Not sure how this can be done...Maybe would need a central server
  #   that can be queried by FTP. I can't find one from NOAA...)
  # Later, multiple instruments (time series, one-to-one, map, etc.)
  #   Scaling factors for multiple


#%%
def plot_spt_samples(oof_ffname):
  """
  Creates a list of *.spt files to plot (2 per window), then call the subfunction
  which actually plots them and saves them to a *.pdf
  """
  xxYYYYMMDD = os.path.basename(oof_ffname)[0:10]
  plot_dir = os.path.dirname(oof_ffname) + '/plots/' 
  spt_path = os.path.dirname(oof_ffname) + '/spt/'
  pdf_name = plot_dir +  xxYYYYMMDD + '_spt_exmpls.pdf' #taking spt directory name for file name part
  
  if os.path.isdir(spt_path): #see if it's a directory
    gaswin_dirs = glob.glob(spt_path + '*_*' + '/') #gas+window directories
    spt_plt_names = [[None, None] for i1 in range(len(gaswin_dirs)) ] #to contain filenames to plot
    
    for i1, i_dir in enumerate(gaswin_dirs): #each gas+window
      gw_files = glob.glob(i_dir + 'z' + xxYYYYMMDD + '*.*')
      gw_files.sort()
      
      if len(gw_files)>0:
        spt_plt_names[i1][0] = gw_files[0] #use first spectrum, presumably aligned?
        if len(gw_files)>1:
          spt_plt_names[i1][1] = gw_files[len(gw_files)//2] #approximately median file
        
  #      pdf_name = plot_dir + os.path.basename(i_dir[0:-1]) + '_exspt.pdf' #taking spt directory name for file name part
    
    spt_plt_names = [x for x in spt_plt_names if x != [None, None] ] #remove empty lines
    if len(spt_plt_names)>0:
      plot_2spt_functions_pdf(spt_plt_names, pdf_name)



#%%
def plot_2spt_functions_pdf(names_lst, out_pdf):
  """
  This function takes a list (names_lst) where each element contains a sublist 
  of length 2, where the elements are either names to *.spt files or are "None"
  It plots 2 of the spectral fits to a *.pdf page in out_pdf.
  
  Based in part on plot_2spt_function_pdf.m
  """
  non_gas = ['Freq', 'Freq_min', 'Freq_max', 'npts', 'ASZA', 
             'obalt', 'p_alt', 'rms', 'tot_cols', 'filename', 'Tm', 'Tc', 'Cont'] #'solar', 'Cont', 
  clr_lst=plt.rcParams['axes.prop_cycle'].by_key()['color'] #default color cycle

  with PdfPages(out_pdf) as pdf:
    for i_row in names_lst:
      plt.figure(figsize=(8.5, 11))
      i_row_nemp = [x for x in i_row if x != None ]
      hsp = [None]*4 #default, 4 subplots with the 2 residuals
      if len(i_row_nemp)==1:
        hsp = [None]*2
        
      for i2, i_fname in enumerate(i_row_nemp):
        spt_raw = rgggm.read_spt(i_fname)
        if np.all(np.isnan(spt_raw['Freq'])) or np.all(np.isnan(spt_raw['Tc'])):
          break #skip NaN windows.
          
        
        fn = list(spt_raw.keys())
        fn_gas = [x for x in fn if x not in non_gas] #field names of gases
        
#        nfg = len(fn_gas) #number of fitted gases + solar + Tm & Tc + Cont
        
        maingas = fn_gas[0] #
        ctr_wn = np.floor(np.median(spt_raw['Freq'])) #center wavenumber
        
        
        ########### Plot residuals #############
        hsp[i2*2] = plt.subplot(len(hsp), 1, i2*2+1)
        resid_pct = (spt_raw['Tm']-spt_raw['Tc'])*100 #calculate the residuals in percent
        plt.plot(spt_raw['Freq'], resid_pct, 'k', linewidth=1.5)
        plt.grid(True) 
        #box on #????
        plt.ylabel('%', size=12, weight='bold')
        
        plt.title(spt_raw['filename'] + ',  ' + '%.3f' % spt_raw['ASZA'] +  '$\mathbf{^\circ}$ SZA,  ' + 
                  '%.4f' % (spt_raw['rms']*100) + '% RMS,  ' + '%.4E' % spt_raw['tot_cols'][0] + 
                  ' molec$\mathbf{\cdot cm^{-2}}$', size=12, weight='bold')
        
        plt.xlim([spt_raw['Freq_min'], spt_raw['Freq_max']])
        plt.tick_params(axis='x', which='both', labelbottom=False)
        
        ## find the largest residual, for "smart" assigning of ylimits on residuals
        resid_lim = np.abs(resid_pct);
        if len(resid_lim)>100:
          resid_lim = np.sort(np.abs(resid_lim))
          resid_lim = resid_lim[np.int32(0.995*len(resid_lim))] # takes a given percentile as maximum
        else:
          resid_lim = np.max(np.abs(resid_lim))
        
        
        if resid_lim>1:
          resid_lim=np.ceil(resid_lim)
        else: # less than 1, round to nearest 10
          resid_lim=np.ceil(resid_lim*10)/10

        plt.ylim(np.array([-1, 1])*resid_lim)
        plt.text(spt_raw['Freq_max']+(spt_raw['Freq_max']-spt_raw['Freq_min'])*0.02, 0, 
                 maingas + '\n' + '%d' % ctr_wn, fontsize=12, weight='bold') #add gas and window
        
        
        ############ Plot main species #############
        hsp[i2*2+1] = plt.subplot(len(hsp), 1, i2*2+2)
        plt.plot(spt_raw['Freq'], spt_raw['Tm'], 'k.', markersize=5, label='Tm') #plot measured points
        plt.plot(spt_raw['Freq'], spt_raw['Tc'], 'k', linewidth=1.5, label='Tc') #plot calculated points
        plt.plot(spt_raw['Freq'], spt_raw[fn_gas[0]], '-', color=clr_lst[0], linewidth=2, label=fn_gas[0]) #main gas
        for i3, i_spec in enumerate(fn_gas[1:]):
          ln_sty = '-'
          if i3>=len(clr_lst)-1:
            ln_sty = '--'
          plt.plot(spt_raw['Freq'], spt_raw[i_spec], ln_sty, color=clr_lst[(i3 % (len(clr_lst)-1))+1], 
                   linewidth=1.5, label=i_spec)
        
        plt.grid(True)
        
        
        ### set x & y-limits, add labels
        y_min = np.min(np.stack((spt_raw['Tm'], spt_raw['Tc']))) 
        y_min = np.floor(y_min*20)/20 #floor to nearest 0.05
      # y_max = np.max(np.stack((spt_raw['Tm'], spt_raw['Tc']))) 
      # if y_max>1: y_max=1.01; else y_max =1 #set upper limit
        y_max = 1.02;
        plt.ylim([y_min, y_max])
        plt.ylabel('Transmittance')
        plt.xlim([spt_raw['Freq_min'], spt_raw['Freq_max']])
        plt.xlabel('Wavenumber (cm$^{-1}$)')
        
        ### Add legend
        hsp[i2*2+1].legend(loc='center left', bbox_to_anchor=(1,0.5), borderaxespad=0,
           frameon=False)
        
        
        
        ### Manually set subplot positions
        hsp[i2*2].set_position([0.10, 0.33 + (len(hsp)*0.5-i2-1)*0.5, 0.78, 0.10])
        hsp[i2*2+1].set_position([0.10, 0.08 + (len(hsp)*0.5-i2-1)*0.5, 0.78, 0.25])
        
        
      
      ### Save current page
      # If there are all 'None' elements then don't save the page
      if not all(v is None for v in hsp):
        pdf.savefig()  # saves the current figure into a pdf page
      plt.close()


#%% plot daily summaries
def make_summary_pdf(oof_ffname, in_flag_list=['flag', 'eq', 0]):
  """ This makes a summary *.pdf page with various statistics and quick info to 
  help with diagnosis for any potential issues for the day as well as a reminder.
  Includes: number of files (including incomplete and cloudy), who ran retrievals,
  when, and on what system. Includes a map or at least a Google URL
  Includes runtimes, and SIA timeseries,
  Includes a listing of any overrides, and a histogram of flag values
  Also includes some stats of a few select species
  Example:
  oof_ffname = '/home/osboxes/research/tccon/egi2020/results/cn_spec/daily/20140609/cn20140609.vav.ada.aia.oof'
  """
  

  oof_raw = rgggm.read_oof(oof_ffname)
  oof_hdr = rgggm.parse_oof_comments(oof_raw['comments']) #parse the header comments
  
  xxYYYYMMDD = os.path.basename(oof_ffname)[0:10]
  plot_dir = os.path.dirname(oof_ffname) + '/plots/' 
  pdf_name = plot_dir +  xxYYYYMMDD + '_summary.pdf'
  
  mm_fname = os.path.dirname(oof_ffname) + '/matchup/' + xxYYYYMMDD + '_mmtch.dat'
  
  # 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_ln_nRaw = next(i1 for i1, val in enumerate(mm_raw) if val[0:5] == 'nRaw:')
  splt_line = re.split('[ ]+', mm_raw[n_ln_nRaw].strip())
  nRaw = np.int64(splt_line[1])
  inst_config = mm_raw[0].strip() #strip off the new line
  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) 
  
  
  
  with PdfPages(pdf_name) as pdf:
    ########## First page: general info  #######
    fig, ax = plt.subplots(4,2, figsize = (8.5, 11))
  #  plt.figure(figsize=(8.5, 11))
    
    
    ###### Big title #####
    plt.suptitle('Summary for ' + xxYYYYMMDD[0:2] + ' on ' + xxYYYYMMDD[2:], fontsize = 20)
    
    ##### Summary text/table #####
    ax1 = plt.subplot2grid((4,2), (0,0))
    ax1.axis('off')
    table1 = ax1.table(cellText=[['Retrievals were run by ' + oof_hdr['EGI']['name']], 
                        ['On system: ' + oof_hdr['EGI']['system']], 
                        ['Finished: ' + oof_hdr['EGI']['date']],
                        [inst_config],
                        ['%6d' % nRaw  + '   raw files were collected'],
                        ['%6d' % np.sum(mm_df['Iflag']==1)  + '   were incomplete'],
                        ['%6d' % np.sum(mm_df['Iflag']==2)  + '   were too cloudy/low signal'],
                        ['%6d' % len(oof_raw['flag'])  + '   output record lines/spectra'],
                        ['%6d' % np.sum(oof_raw['flag']==0)  + '   were flag==0']], 
                        edges='', loc='center', cellLoc='left')
    
    table1.set_fontsize(13)
    table1.scale(1.3, 1.3)
    
    
    ##### Maps #####
    # Qualities I wanted were no extra modules that were not built-in to Anaconda
    # and no required API key. I couldn't find a perfect solution, so making a URL 
    # for Google and getting a map image from Russia's Yandex...
    ## Map points, center, and zoom level
    unq_lls=np.unique( np.vstack((np.round(oof_raw['lat(deg)']*1000)/1000, np.round(oof_raw['long(deg)']*1000)/1000)).T, axis=0)
  #  if len(np.shape(unq_lls))==1: #only 1 location
  #    unq_lls = unq_lls[np.newaxis, :]
    
    url_lls = [['%.3f' % unq_lls[i1,0], '%.3f' % unq_lls[i1,1]] for i1 in range(np.shape(unq_lls)[0]) ]
    ctr_ll = ['%.3f' % np.median(unq_lls[:,0]), '%.3f' % np.median(unq_lls[:,1])]
    
    poor_maxd_km = np.array([np.max(np.sqrt((unq_lls[:,0]-np.median(unq_lls[:,0]))**2) ),
                    np.max(np.sqrt((unq_lls[:,1]-np.median(unq_lls[:,1]))**2) )])*111 #poor estimate of max lat & lon dist from median
    meters_per_pixel = np.max(poor_maxd_km*1000/np.array([180, 140])) #+/-200 or +/-150 pixels each direction, so using 180 or 140 to keep away from edges
    zoom_lvl = meters_per_pixel / np.abs( 156543 * np.cos(np.median(unq_lls[:,0]) * np.pi / 180) ) + 1E-9
    zoom_lvl = np.min([15, np.ceil(-1*np.log2(zoom_lvl)) ]) #don't zoom in past level ~15
    
    ## Google link URL
    google_url = 'https://www.google.com/maps/dir/'
    for i1_ll in url_lls:
      google_url = google_url + i1_ll[0] + ',' + i1_ll[1] + '/' #add marker point
    google_url = google_url + '/@' + ctr_ll[0] + ',' + ctr_ll[1] #"@" indicates don't actually get directions between
    # could also add zoom level here...
    
    ## Yandex link/figure URL
    #  https://staticmapmaker.com/yandex/ - Yandex is Russian, but does not require API like Google...
    # much better street level detail in Europe than NA
    yandex_clrs = {'white': 'wt', 'dark orange': 'do', 'dark blue': 'db', 'blue': 'bl', 
                   'green': 'gn', 'grey': 'gr', 'light blue': 'lb', 'night': 'nt', 
                   'orange': 'or', 'pink': 'pn', 'red': 'rd', 'violet': 'vv', 
                   'yellow': 'yw'}
    yandex_c2 = [yandex_clrs[v] for v in yandex_clrs.keys()] #just the 2 letter ids
    yandex_url = 'https://static-maps.yandex.ru/1.x/?lang=en-US&ll='
    yandex_url = yandex_url + ctr_ll[1] + ',' + ctr_ll[0] 
    
    yandex_url = yandex_url + '&z=' + '%d' % zoom_lvl
    yandex_url = yandex_url + '&l=map&size=400,300&pt='
    for i1, i1_ll in enumerate(url_lls):
      # will have problems in the case of more than 99 points
      yandex_url = yandex_url + i1_ll[1] + ',' + i1_ll[0] + ',pm' + yandex_c2[i1 % len(yandex_c2)] + 's' + '%d' % (i1+1) + '~'
    yandex_url = yandex_url[:-1] #trim off last tilde
    
    ## Try to download yandex map
    try:
      yandex_mfile = plot_dir + xxYYYYMMDD + '_yandex_map.png'
      if not os.path.isfile(yandex_mfile): #don't get map again if you already have it
        with open(yandex_mfile, 'wb') as local_map:
          resp = requests.get(yandex_url, stream=True) #set stream to True, this will return the stream content.
          
          # Set decode_content value to True, otherwise the downloaded image file's size will be zero.
          resp.raw.decode_content = True
          # Copy the response stream raw data to local image file.
          shutil.copyfileobj(resp.raw, local_map)
          # Remove the image url response object.
          del resp
      got_ymap = True
    except:
      got_ymap = False
  
    ## Add map and URL
    ax2 = plt.subplot2grid((4,2), (0,1))
    ax2.axis('off')
    if got_ymap:
      img_raw=mpimg.imread(yandex_mfile)
      ax2.imshow(img_raw)
    
    plt.text(0.5, 0, google_url, transform=ax2.transAxes, fontsize=6, verticalalignment='top',
             horizontalalignment='center')
    
    
    
    ##### Runtimes #####
    utc_sns = egs.sunrise_noon_sunset(np.median(oof_raw['day']), np.median(oof_raw['long(deg)']), np.median(oof_raw['lat(deg)']))
    dt_obj = dt.datetime(np.int64(np.median(oof_raw['year'])), 1, 1, 0, 0, 0) + dt.timedelta(np.median(oof_raw['day']) - 1) #
    dts_sr_noon_ss = [dt_obj + dt.timedelta(hours=v) for v in utc_sns]
    dts_ss = [dt_obj + dt.timedelta(hours=np.min(oof_raw['hour'])), dt_obj + dt.timedelta(hours=np.max(oof_raw['hour']))]
    dt_duration = dts_ss[1] - dts_ss[0]
    dur_HMS = [dt_duration.seconds//3600, 0, 0]
    dur_HMS[1] = (dt_duration.seconds-dur_HMS[0]*3600)//60
    dur_HMS[2] = dt_duration.seconds - dur_HMS[0]*3600 - dur_HMS[1]*60
    
    ax3 = plt.subplot2grid((4,2), (1,0))
    ax3.axis('off')
    table3 = ax3.table(cellText=[['Approx sunrise (UTC): ' + dts_sr_noon_ss[0].strftime('%H:%M:%S') ], 
                        ['Approx noon (UTC): ' + dts_sr_noon_ss[1].strftime('%H:%M:%S') ],
                        ['Approx sunset (UTC): ' + dts_sr_noon_ss[2].strftime('%H:%M:%S') ],
                        ['Approx start (UTC): ' + dts_ss[0].strftime('%H:%M:%S') ], 
                        ['Approx end (UTC): ' + dts_ss[1].strftime('%H:%M:%S') ], 
                        ['Duration (hours): ' + '%02d:%02d:%02d' % (dur_HMS[0], dur_HMS[1], dur_HMS[2])]], 
                        edges='', loc='center', cellLoc='left')
    
    table3.set_fontsize(12)
    table3.scale(1.2, 1.2)
    
    
    
    ##### Time series of signal #####
    plt.subplot2grid((4,2), (1,1))
    plt.plot(oof_raw['hour'], oof_raw['sia(AU)'], '.')
    plt.xlabel('Time (UTC)')
    plt.ylabel('Signal (AU)')
    plt.grid(True)
    
    
    ##### Overrides table #####
    ax5 = plt.subplot2grid((4,2), (2,0))
    ax5.axis('off')
    tbl5_txt = [['*.mod and *.vmr notes:'],
                [oof_hdr['EGI']['mod_notes']],
                ['Overrides:']] + \
                oof_hdr['EGI']['overrides']
  #  tbl5_clrs = ['r' if 'Warning' in s[0] else 'w' for s in tbl5_txt] #did not work. :(
    
    table5 = ax5.table(cellText=tbl5_txt, 
                        edges='', loc='center', cellLoc='left')
    table5.set_fontsize(13)
    table5.scale(1.5, 1.5)
    
    
    ##### Barchart of flags #####
    ax6 = plt.subplot2grid((4,2), (2,1))
    unq_vals, unq_cnts = np.unique(oof_raw['flag'], return_counts=True)
    
    x_idx = np.arange(len(unq_vals)) - 0.5
    ax6.bar(x_idx, unq_cnts)
    ax6.set_ylabel('Flag counts')
    ax6.set_xticks(x_idx)
    ax6.set_xticklabels(['%d' % s for s in unq_vals])
    
    
    
    ##### Summary table of some stats #####
    ax7 = plt.subplot2grid((4,2), (3,0), colspan=2)
    ax7.axis('off')
    if 'xch4(ppm)' in oof_raw:
      oof_raw['xch4(ppb)'] = oof_raw['xch4(ppm)']*1000
    
    t7_flds = ['xluft', 'xco2(ppm)', 'xch4(ppb)', 'xco(ppb)', 'xn2o(ppb)', 'xh2o(ppm)']
    t7_fmts = ['%.4f', '%.2f', '%.1f', '%.1f',  '%.1f',  '%d']
    
    # Might want to consider allowing flagging on exclusively eof fields?
#    bool_flag = oof_raw['flag']==0
    bool_flag = egs.oof_bool_flag(oof_raw, in_flag_list) 
    if np.sum(bool_flag)>0:
      tbl7_txt = [['field', 'mean', 'std', 'median']]
      for i1, iFld in enumerate(t7_flds):
        if iFld in oof_raw:
          if not np.all(np.isnan(oof_raw[iFld][bool_flag])): #prevents errors when NaN for all entries
            tbl7_txt.append([iFld, t7_fmts[i1] % np.nanmean(oof_raw[iFld][bool_flag]), t7_fmts[i1] % np.nanstd(oof_raw[iFld][bool_flag]),
                           t7_fmts[i1] % np.nanmedian(oof_raw[iFld][bool_flag]) ])
          else:
            tbl7_txt.append([iFld, t7_fmts[i1] % np.NaN, t7_fmts[i1] % np.NaN,
                           t7_fmts[i1] % np.NaN ])
            
      # Might also consider including table with all NaNs if data were flagged? or some note?
      table7 = ax7.table(cellText=tbl7_txt, 
                          edges='BT', loc='center', cellLoc='right')
      table7.set_fontsize(12)
      table7.scale(1.1, 1.1)
    
    
    ##### Add advertisment #####
    plt.figtext(0.99, 0.01, 'Created with EGI2020', horizontalalignment='right')
    
    ###### Save page 1 #######
    plt.subplots_adjust(bottom=0.04, top=0.93, hspace=0.28, wspace=0.3)
    
    
    pdf.savefig(dpi=300)  # saves the current figure into a pdf page. DPI needed for map
    plt.close()
  
  

#%% plot various time series
# later idea: possibly add ~5 minute averages

#oof_ffname = '/home/osboxes/research/tccon/egi2020/results/cn_spec/daily/20140609/cn20140609.vav.ada.aia.oof'
#oof_ffname = '/media/sf_share_mint_19_2/results/ta_oof_plotting/ta20200516.vav.ada.aia.oof'
def plot_misc_tseries(oof_ffname, in_flag_list=['flag', 'eq', 0]):
  ###### Load raw data #######
  oof_raw = rgggm.read_oof(oof_ffname)
  oof_hdr = rgggm.parse_oof_comments(oof_raw['comments']) #parse the header comments
  
  
  xxYYYYMMDD = os.path.basename(oof_ffname)[0:10]
  oof_dir = os.path.dirname(oof_ffname)
  plot_dir = os.path.dirname(oof_ffname) + '/plots/' 
  pdf_name = plot_dir +  xxYYYYMMDD + '_timeseries1.pdf'
  
#  if oof_ffname[-4:] == '.oof':
##    eof_ffname = oof_ffname[:-4] + '.eof'
#    oof_part_fname = oof_ffname[:-16]
#  elif oof_ffname[-8:] == '.oof.csv':
##    eof_ffname = oof_ffname[:-8] + '.eof.csv'
#    oof_part_fname = oof_ffname[:-20]
#  else:
#    raise ValueError("I don't understand extension on oof file: " + oof_ffname)
##  eof_raw = rgggm.read_eof(eof_ffname)
  
  gen_nc_fname = oof_dir + '/' + xxYYYYMMDD[0:2] + '*.private.nc'
  nc_ffnames = glob.glob(gen_nc_fname)
  if len(nc_ffnames) > 1:
    raise ValueError('Got ' + '%d' % len(nc_ffnames) + ' matches for ' + gen_nc_fname + ' Should be just 1')
  elif len(nc_ffnames) == 0:
    # netCDF file was not created
    raise ValueError('No matches for %r' % gen_nc_fname)
  else:
    nc_ffname = nc_ffnames[0] #just the 1 file...
  
  eof_raw = rgggm.read_prv_nc(nc_ffname)
  
  
  ####### Make listing of fields to plot timeseries of #######
  plt_flds = ['xluft', 'xco2(ppm)', 'xch4(ppb)', 'xco(ppb)', 'xn2o(ppb)', 'xh2o(ppm)']
  #plt_flds = plt_flds + ['xwco2(ppm)', 'xhdo(ppm)']
  plt_flds = plt_flds + ['fvsi(%)', 'sia(AU)', 'solzen(deg)', 'osds(ppm)']
  plt_flds = plt_flds + [['pins(mbar)', 'pout(hPa)'],  'hout(%RH)', ['tins(C)', 'tout(C)']]
  plt_flds = plt_flds + ['wind'] #special field
  plt_flds = plt_flds + ['co2_6220_fs', 'o2_7885_am_o2'] #netcdf fields are all lower case, *.eof fields had some upper case
  
  ###### make boolean flag, and time vector ######
#  bool_flag_cmb = oof_raw['flag']==0
  bool_flag_cmb = egs.oof_bool_flag(oof_raw, in_flag_list, eof_raw) 
  x_time_oof = [dt.datetime(int(oof_raw['year'][i1]), 1, 1) + dt.timedelta(days=oof_raw['day'][i1]-1) + dt.timedelta(hours=oof_raw['hour'][i1]) for i1 in range(len(bool_flag_cmb))]
  x_time_eof = [dt.datetime(int(eof_raw['year'][i1]), 1, 1) + dt.timedelta(days=np.float64(eof_raw['day'][i1]-1)) + dt.timedelta(hours=np.float64(eof_raw['hour'][i1])) for i1 in range(len(eof_raw['hour']))]
  
  clr_lst1 = matplotlib.cm.get_cmap('tab20')
  clr_lst1 = clr_lst1(np.linspace(0.0,1,20))
  mm_times = [dt.datetime(int(np.min(oof_raw['year'])), 1, 1) + dt.timedelta(days=np.min(oof_raw['day'])-1) + dt.timedelta(hours=np.min(oof_raw['hour']) ),
      dt.datetime(int(np.max(oof_raw['year'])), 1, 1) + dt.timedelta(days=np.max(oof_raw['day'])-1) + dt.timedelta(hours=np.max(oof_raw['hour']) )]
  time_diff = (mm_times[1]-mm_times[0]).seconds
  if time_diff < 300:
    xTickFmt = mdates.DateFormatter('%H:%M:%S')
  else:
    xTickFmt = mdates.DateFormatter('%H:%M')
  
  ###### Loop through fields and plot them ######
  #%%
  maxOnPlot = 4 #how many time-series per plot
  runTot = 0 #running total number of sub-plots
  with PdfPages(pdf_name) as pdf:
    nOnPlot = 0 #number on plot
    for iField in plt_flds:
      ##### Initialize new plot ######
      if nOnPlot == 0: #initialize new plot
        fig, ax = plt.subplots(maxOnPlot,1, figsize = (8.5, 11))
        hsp = [None]*maxOnPlot
        
        
      ##### Reset boolean flag and x-data (this is to allow incomplete eof data to also occasionally be plotted) #####
      plt_bool = copy.copy(bool_flag_cmb)
      x_time = copy.copy(x_time_oof)
      
      
      ##### Get data to plot from oof and eof files #####
      iFld_lst = iField
      if type(iField) is not list:
        if iField=='wind': #special case
          iFld_lst = ['wspd(m/s)', 'wdir(deg)']
        else:
          iFld_lst = [iField]
      
      #pre-allocate. Usually a vector (except for when multiple fields are on one plot or special wind case)
      dat_fields = np.zeros((len(bool_flag_cmb), len(iFld_lst))) 
      flag_lims = np.zeros((2, len(iFld_lst)))*np.NaN
      used_flds = np.ones(len(iFld_lst), dtype=bool)
      for i1, iDatF in enumerate(iFld_lst):
        if iDatF in oof_raw: #take from *.oof file first
          dat_fields[:,i1] = oof_raw[iDatF]
          flg_fldnm = [k for k,s in enumerate(list(oof_hdr['flags']['Variable'])) if s in iDatF ]
          if len(flg_fldnm)>1: #more than 1 match, like xco for xco2
            fldnm_len = np.array([len(oof_hdr['flags']['Variable'][k]) for k in flg_fldnm])
            flg_fldnm = [flg_fldnm[np.argmax(fldnm_len)]] #saying the longest match is the best...
            
          if len(flg_fldnm)==1:
            flag_lims[:,i1] = [oof_hdr['flags'].loc[flg_fldnm[0],'Vmin'], oof_hdr['flags'].loc[flg_fldnm[0],'Vmax']]
          
        elif iDatF in eof_raw: #then try *.eof
          if len(eof_raw[iDatF]) == len(bool_flag_cmb): #this is the normal case, otherwise the *.eof has a differing number of lines from *.oof
            dat_fields[:,i1] = eof_raw[iDatF]
          elif len(eof_raw[iDatF])>0 and len(iFld_lst)==1: #if it's just 1 field that's not empty, I'll let it slide that it doesn't match in length
            dat_fields = np.zeros((len(eof_raw[iDatF]), len(iFld_lst)))
            dat_fields[:,i1] = eof_raw[iDatF]
            x_time = copy.copy(x_time_eof)
            # nothing gets flagged for now...later could make separate flags for oof, eof, plus one combined
            plt_bool = np.ones(len(eof_raw[iDatF]), dtype=bool) 
            
          else: #not going to bother trying to include field of different length
            used_flds[i1] = False
        else: #just skip field
          used_flds[i1] = False
      
      if np.sum(used_flds)==0:
        continue #skip to next iteration
      
      dat_fields = dat_fields[:, used_flds]
      iFld_lst = list(compress(iFld_lst, used_flds))
      flag_lims = flag_lims[:, used_flds]
      
      ###### Rescale data if it is outside of flag limits ######
      if iField=='wind': #special case
        flag_lims = np.zeros((2, len(iFld_lst)))*np.NaN
      
      
      out_of_lim = [False, False]
      if not np.any(np.isnan(flag_lims[:,0])): #just doing the first field
        mm_diff = flag_lims[1,0] - flag_lims[0,0] #difference
        plt_mxmn = flag_lims[:,0] + np.array([-0.1, 0.1])*mm_diff
        dat_fields[dat_fields<plt_mxmn[0]] = plt_mxmn[0]
        dat_fields[dat_fields>plt_mxmn[1]] = plt_mxmn[1]
        if np.any(dat_fields<flag_lims[0,0]):
          out_of_lim[0] = True
        if np.any(dat_fields>flag_lims[1,0]): 
          out_of_lim[1] = True
      
      #### Add title if it's the first page ####
      if runTot == 0:
        plt.suptitle('Timeseries for ' + xxYYYYMMDD[0:2] + ' on ' + xxYYYYMMDD[2:], fontsize = 20)
      
      ##### Plot data #####
      hsp[nOnPlot] = plt.subplot(maxOnPlot, 1, nOnPlot+1)
      if iField=='wind': #special case. Maybe will need to plot these less often to see direction?
        if np.sum(~bool_flag_cmb)>0:
          pt_idx = np.where(~bool_flag_cmb)[0]
          plt_yet = False
          for i2,i1 in enumerate(pt_idx):
            lgd_lbl=None
            if dat_fields[i1,0]>0:
              w_angls = (90-dat_fields[i1,1]+np.array([-15, 15])) % 360 #pointer goes in direction opposite to input value
              i2_mkr = mpath.Path.wedge(w_angls[0], w_angls[1]) 
            else: #zero wind speed, plot as a dot
              i2_mkr = '.'
            if not plt_yet:
              lgd_lbl=iField+' flgd'
            plt_yet = True
            hsp[nOnPlot].plot(x_time[i1], dat_fields[i1,0], marker=i2_mkr, linestyle='None', color=clr_lst1[1,:], label=lgd_lbl)
        if np.sum(bool_flag_cmb)>0:
          pt_idx = np.where(bool_flag_cmb)[0]
          plt_yet = False
          for i2,i1 in enumerate(pt_idx):
            lgd_lbl=None
            if dat_fields[i1,0]>0:
              w_angls = (90-dat_fields[i1,1]+np.array([-15, 15])) % 360
              i2_mkr = mpath.Path.wedge(w_angls[0], w_angls[1])
            else:
              i2_mkr = '.'
            if not plt_yet:
              lgd_lbl=iField
            plt_yet = True
            hsp[nOnPlot].plot(x_time[i1], dat_fields[i1,0], marker=i2_mkr, linestyle='None', color=clr_lst1[0,:], label=lgd_lbl)
        
        hsp[nOnPlot].set_ylabel('Wind speed [m/s]')
        
      else:
        for i1, iDatF in enumerate(iFld_lst):
          if np.sum(~plt_bool)>0: #plot flagged out data
            hsp[nOnPlot].plot(list(compress(x_time, ~plt_bool)), dat_fields[~plt_bool,i1], '.', color=clr_lst1[i1*2+1,:], label=iDatF+' flgd')
          if i1==0:
            hsp[nOnPlot].set_ylabel(iDatF)
          if np.sum(plt_bool)>0:
            hsp[nOnPlot].plot(list(compress(x_time, plt_bool)), dat_fields[plt_bool,i1], '.', color=clr_lst1[i1*2,:], label=iDatF)
      
      ##### Set y-limits, add flag line if there are data out of range
      y_lims = np.array(hsp[nOnPlot].get_ylim())
      if out_of_lim[0]:
        y_lims[0] = plt_mxmn[0]
        y_lims[1] = np.max([flag_lims[0,0]+mm_diff*0.1, y_lims[1]]) 
        plt.axhline(y=flag_lims[0,0], color='b', linestyle=':')
      if out_of_lim[1]:
        y_lims[1] = plt_mxmn[1]
        y_lims[0] = np.min([flag_lims[1,0]-mm_diff*0.1, y_lims[0]]) 
        plt.axhline(y=flag_lims[1,0], color='r', linestyle=':')
  
      hsp[nOnPlot].set_ylim(y_lims)
      if runTot==0: #luft
        x_lims = hsp[nOnPlot].get_xlim()
      else: #useful for incomplete *.eofs
        hsp[nOnPlot].set_xlim(x_lims)
      
      
      hsp[nOnPlot].grid(True)
      hsp[nOnPlot].legend()
      
      hsp[nOnPlot].xaxis.set_major_formatter(xTickFmt)
      ##### Clear last x-tick labels #####
      if nOnPlot>0:
        hsp[nOnPlot-1].xaxis.set_ticklabels([])
      
      ##### Advance loop, plot current page if needed #####
      runTot = runTot+1
      nOnPlot = nOnPlot+1
      if nOnPlot==maxOnPlot: #reached capacity on current figure
        plt.setp(plt.xticks()[1], rotation=30, ha='right')
        nOnPlot=0
        pdf.savefig()  # saves the current figure into a pdf page
        plt.close()
      
      
        
    
    ###### Save final figure ######
    if nOnPlot > 0:
      plt.setp(plt.xticks()[1], rotation=30, ha='right')
      for i1 in range(len(hsp)):
        if hsp[i1] is None:
          ax_i1 = plt.subplot(maxOnPlot, 1, i1+1)   
          ax_i1.axis('off')
          
      
      pdf.savefig()  # saves the current figure into a pdf page
    plt.close()
    


  



#%%
# Other plots or subplots:
# Something with AOD??
# Allan variance?
# Quantile plots?
# A priori profiles?
# HYSPLIT trajectories? (Not sure how this can be done...Maybe would need a central server
#   that can be queried by FTP. I can't find one from NOAA...)
# Later, multiple instruments (time series, one-to-one, map, etc.)
#   Scaling factors for multiple







