Hydrolight53_Bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# HE53_Bootstrap\n",
    "# Jesse Bausell\n",
    "# 18 September 2019\n",
    "#\n",
    "# HE53_Bootstrap enables user to evaluate radiative transfer using field-measured inherent optical properties (IOPs). \n",
    "# The program accepts ac-s and hs6 data and performs bootstrapping with replacement to produce a series of depth-binned\n",
    "# absorption, attenuation, and backscattering from single spectra. The resulting output files can be uploaded into \n",
    "# Hydrolight formodeling radiative transfer. Program is specifically designed to work with Seabird Scientific ac-s \n",
    "# meter and HOBI Labs Hydroscat6 meter. See readme for more details.\n",
    "#\n",
    "# Inputs: \n",
    "    # ac-s_Data - single absorption/attenuation spectra. Program accepts ascii or hdf5/mat files.\n",
    "    # hs6_Data - single particulate backscattering spectra. Program accepts ascii or hdf5/mat files\n",
    "    # UserTemplate_ascii.txt - ascii file with user-specified data. Template provided.\n",
    "# Outputs:\n",
    "    # STATION_ac-s_bin##_itr#.txt - hydrolight compatible depth-binned absorption/attenuation\n",
    "    # STATION_hs6_bin##_itr#.txt - hydrolight compatible depth-binned particulate backscattering\n",
    "    # STATION_bin##_itr#.txt - # batch file used to run hydrolight using aforementioned depth-binned IOPs\n",
    "# Please note:\n",
    "# Bootstrapping with replacement is repeated n times\n",
    "# ##- denotes user-specified depth bin size\n",
    "#  #- denotes bootstrap iteration (e.g. 1, 2, 3 ... n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import python libraries\n",
    "import tkinter as tk # Import the proper libraries for inserting a file\n",
    "from tkinter import filedialog # Filedialog allows user to select an input file\n",
    "from datetime import datetime as dt # \n",
    "import numpy as np # Import the numpy as np\n",
    "import h5py # Import hdf library (used for mat/dfh5 files)\n",
    "import re # Import regular expression library\n",
    "import csv # Import csv reader/writer library\n",
    "import os # Import miscellaneous operating systems interface library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import user-defined functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def UserTemplate(file_path):\n",
    "    \"\"\"Opens user template and places input data into a python dictionary.\n",
    "    Inputs: \n",
    "        file_path - file pathway and name for header file (see readme)\n",
    "    Outputs:\n",
    "        User_Dict - python dictionary containing all user-supplied specifications\"\"\"\n",
    "    keY = 0 # Key to regulate while loop\n",
    "    User_Dict = {} # Create dictionary for user-input data\n",
    "    with open(file_path) as raw: # Open User template \n",
    "        while 1:\n",
    "            # This while loop cycles through the header file lines, collects important info, \n",
    "            # and puts it into a dictionary\n",
    "            title = raw.readline()\n",
    "            if 'Station Name' in title:\n",
    "                # For the first line of important information\n",
    "                keY += 1 # increase keY by one.\n",
    "            if keY > 0:\n",
    "                # If the first line of important information has been reached, do these steps \n",
    "                # for each subsequent while-loop iteration...\n",
    "                IND = title.index(':') # Find index of ':'\n",
    "                User_Dict[title[:IND]] = title[IND+2:-1] # Take information right of \":\". Exclude \"\\n\"\n",
    "                # and put it into a dictionary\n",
    "                if 'Bootstrap Iterations' in title:\n",
    "                    # Once the last important header line is reached, break the while-loop\n",
    "                    break\n",
    "        return(User_Dict)\n",
    "             "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def HDF_LOADER(hdf_name,**parameters):\n",
    "    \"\"\"This function loads processed ac-s data (individual spectra) from .mat/hdf5 files \n",
    "    that will later be bootstrapped. It is designed for data prepared by the matlab programs \n",
    "    acsPROCESS_INTERACTIVE.m or hs6PROCESS_INTERACTIVE.m. User should keep in mind that when \n",
    "    using hs6 (backscattering) data, A_CORR and C_CORR output variables are redundant.\n",
    "    Inputs:\n",
    "        hdf_name - name of hdf5/mat file\n",
    "    Outputs\n",
    "        A_CORR - absorption or particulate backscattering data\n",
    "        C_CORR - attenuation or particulate backscattering data\n",
    "        deptH - ac-s or hs6 depths\n",
    "        lambdA - ac-s or hs6 wavelengths\"\"\"\n",
    "    INST = parameters['INST'] # Define instrument as ac-s or hs6\n",
    "    with h5py.File(hdf_name) as hf:\n",
    "        # Assign variables to hdf fields\n",
    "        if 'ac-s' in INST:\n",
    "            A_CORR = hf['A_CORR'][()] # Absorption values\n",
    "            C_CORR = hf['C_CORR'][()] # Attenuation values\n",
    "        elif 'hs6' in INST:\n",
    "            A_CORR = hf['HS6_data_bbp'][()] # Absorption values\n",
    "            C_CORR = hf['HS6_data_bbp'][()] # Attenuation values\n",
    "        deptH = hf['deptH'][()] # Depth values\n",
    "        acs_wvl = hf['lambda'][()] # Wavelengths\n",
    "        # Transpose matrices \n",
    "        A_CORR = np.transpose(A_CORR) # Transpose absorption\n",
    "        C_CORR = np.transpose(C_CORR) # Transpose attenuation\n",
    "        # Simplify depth and wavelength arrays\n",
    "        deptH = deptH[0] # Make depth array one dimensional\n",
    "        lambdA = np.zeros(len(acs_wvl)) # Create one dimensional array of zeros the same size as the number of wavelengths\n",
    "        for i,j in enumerate(acs_wvl):\n",
    "            # For-loop assigns a wavelength to each zero value in the newly-created 'lambdA array (of zeros)\n",
    "            lambdA[i] = j[0] # index and assign wavelengths\n",
    "        return(A_CORR,C_CORR,deptH,lambdA)\n",
    "\n",
    "def SEABASS_TXT_LOADER(txt_name):\n",
    "    \"\"\"This function loads processed ac-s or hs6 data from Seabass-formatted .txt files.\n",
    "    It is designed for data prepared by the matlab programs acsPROCESS_SEABASS.m or hs6PROCESS_SEABASS.m.\n",
    "    It is theoretically compatible with HOBI Labs hs4 and hs2 meters, but has not been tested. User should \n",
    "    keep in mind that when using hs6 (backscattering) data, A_CORR and C_CORR output variables are redundant.\n",
    "        Inputs:\n",
    "        hdf_name - name of hdf file\n",
    "    Outputs\n",
    "        A_CORR - absorption or backscattering data\n",
    "        C_CORR - attenuation or backscattering data\n",
    "        deptH - ac-s or hs6 depths\n",
    "        lambdA - ac-s or hs6 wavelengths\"\"\"\n",
    "    with open(txt_name) as raw: # Open user-input SEABASS.txt file\n",
    "        keY = 0 # key variable to distinguish header from data\n",
    "        for title in raw:\n",
    "            # This for-loop reads every line of the SEABASS.txt file one by one\n",
    "            if 'agp' in title: \n",
    "                # SEABASS.txt is an ac-s file, find the header line with column fields\n",
    "                fielDS = title[:-1] # Remove end-of-line character\n",
    "                IND_c = fielDS.find(',cgp') # Find the start of attenuation (cgp) column headers\n",
    "                markeR1 = ',agp' # create marker for absorption (agp) column headers\n",
    "                markeR2 = ',cgp' # create marker for attenuation (cgp) column headers\n",
    "                lambdA_stra = fielDS[:IND_c] # Isolate the section of the header line with only absorpiton fields\n",
    "                lambdA_strc = fielDS[IND_c:] # Isolate the section of the header line with only attenuation fields\n",
    "            elif 'depth,bbp' in title:\n",
    "            # SEABASS.txt is an hs6 file, find the header line with column fields\n",
    "                fielDS = title[:-1] # Remove end-of-line character\n",
    "                IND_c = fielDS.find(',stimf') # Find the start of backscattering (bbp) column headers\n",
    "                markeR1 = ',bbp' # create marker for backscattering column headers\n",
    "                markeR2 = ',bbp' # make secone marker backscattering column headers\n",
    "                lambdA_stra = fielDS[:IND_c] # Isolate the section of the header line with only backscattering fields\n",
    "                lambdA_strc = fielDS[:IND_c] # Provide redundant duplicate to reduce code\n",
    "            elif 'end_header' in title:\n",
    "            # The final line of code has been reached\n",
    "                lambdA_a = lambdA_stra.split(markeR1) # Create a list of absorption (or backscttering) wavelengths \n",
    "                lambdA_a.pop(0) # Remove the first item from the list (because it's junk)\n",
    "                lambdA_a = np.asarray(lambdA_a,dtype=np.float) # Convert list into numpy array\n",
    "                lambdA_c = lambdA_strc.split(markeR2) # Create a list of attenuation (or backscttering) wavelengths\n",
    "                lambdA_c.pop(0) # Remove the first item from the list (because it's junk)\n",
    "                lambdA_c = np.asarray(lambdA_c,dtype=np.float) # Convert list into numpy array\n",
    "                keY +=1 # Because the last line of the metadata has been reached, increase keY variable by 1\n",
    "                data = dict() # Create dictionary to store data\n",
    "                fielDS = fielDS.split(',') # Create list of header line column fields            \n",
    "            elif keY == 1:\n",
    "            # Header metadata has been processed and the for-loop is examining data\n",
    "                title = title.split('\\t') # split line of data into list\n",
    "                for t in fielDS:\n",
    "                    # for-loop takes a line of SEABASS data and organizes it into python dictionary (data)\n",
    "                    IND = fielDS.index(t) \n",
    "                    if t in data.keys():\n",
    "                        data[t].append(title[IND])\n",
    "                    else:\n",
    "                        data[t] = [title[IND]]\n",
    "    deptH = np.asarray(data['depth'],dtype=np.float) # create array of IOP depth values\n",
    "    A_CORR = np.zeros([len(deptH),len(lambdA_a)])*np.nan # create nan matrix for absorption (or backscattering) values\n",
    "    C_CORR = np.zeros([len(deptH),len(lambdA_a)])*np.nan # create nan matrix for attenuation (or redundant backscattering) values\n",
    "    for i,l in enumerate(lambdA_a):\n",
    "        # for-loop fills in A_CORR and C_CORR matrices one column at a time. Wavelengths and markers (agp/cgp or bbp) are\n",
    "        # used to construct dictionary keys and index them onto matrices\n",
    "        try:\n",
    "            # If the wavelenght IS NOT an integer, put it directly into the python dictionary and index data (absorption or backscattering)\n",
    "            A_CORR[:,i] = np.asarray(data[markeR1[1:] + str(l)],dtype=np.float)\n",
    "        except:\n",
    "            # If the wavelenght IS an integer but has a decimal of 0, remove .0 and index data (absorption or backscattering)\n",
    "            A_CORR[:,i] = np.asarray(data[markeR1[1:] + str(int(l))],dtype=np.float)\n",
    "        try:\n",
    "            # If the wavelenght IS NOT an integer, put it directly into the python dictionary and index data (attenuation or backscattering)\n",
    "            C_CORR[:,i] = np.asarray(data[markeR2[1:] + str(lambdA_c[i])],dtype=np.float)\n",
    "        except:\n",
    "            # If the wavelenght IS an integer but has a decimal of 0, remove .0 and index data (attenuation or backscattering)\n",
    "            C_CORR[:,i] = np.asarray(data[markeR2[1:] + str(int(lambdA_c[i]))],dtype=np.float)\n",
    "    return(A_CORR,C_CORR,deptH,lambdA_a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def createFolder(directory):\n",
    "    \"\"\" createFolder searches for a dirctory specified by the user. If there is none, it creates one\"\"\"\n",
    "    try:\n",
    "        if not os.path.exists(directory): # If the folder doesn't exist\n",
    "            os.makedirs(directory) # Create a folder\n",
    "    except OSError: # If there is an error other than a non-existant folder\n",
    "        print ('Error: Creating directory. ' +  directory) # report error and shut down createFolder\n",
    "        \n",
    "def headER_reader(hdf5_FILE,**parameters):\n",
    "    \"\"\"Opens header template file (hdf5 file) and selects the instrument-specific header template \n",
    "    as specified by the user.\n",
    "    Inputs:\n",
    "        hdf5_FILE - hdf5 file containing generic header templates for ac-s, hs6, and hydrolight batch files\n",
    "    Keyword Arguments:\n",
    "        key - specify the header template to retrieve: \"ac-s\", \"hs6\" or \"batch\" \n",
    "    Outputs:\n",
    "        headER_final - python dictionary containing user-specified header template\"\"\"\n",
    "    key = parameters['key'] # Define user-specified keyword argument\n",
    "    with h5py.File(hdf5_FILE) as hf:\n",
    "        # Open header file\n",
    "        headER = hf[key] # Assign dictionary-like data to a variable using the keyword argument\n",
    "        headER_final = {} # Create a blank dictionary for output variable\n",
    "        for i in headER:\n",
    "            # for-loop cycles through each hdf5 dictionary key and places data into an actual dictionary\n",
    "            headER_final[i] = headER[i][()] # Create dictionary element\n",
    "        return(headER_final) # Output final instrument-secific dictionary with generic header template information\n",
    "\n",
    "def Boostrap_FileWriter(BOOTSRAP_MATRIX,User_Dict,**parameters):\n",
    "    \"\"\"Creates a Hydrolight-compatible ascii file using depth-binned IOP data. File is output in a separate folder.\n",
    "    This folder is placed in the same directory as the original file containing non-binned IOP data.\n",
    "    Inputs:\n",
    "        BOOTSRAP_MATRIX - Depth-binned IOP matrix\n",
    "        User_Dict - Python dictionary containing user-specified inputs\n",
    "    keyword arguments:\n",
    "        itR - specify bootstrap iteration. Must be an integer.\n",
    "        wvl - numpy array of IOP wavelengths.\n",
    "        INSTRUMENT - specify whether IOP matrix is 'acs' or 'hs6'\n",
    "    Outputs:\n",
    "        Hyrolight-compatible file of binned IOPs (not a variable)\"\"\"\n",
    "    ### First, import appropriate header template (instrument-dependent)\n",
    "    inST = parameters['INSTRUMENT'] # specify instrument name\n",
    "    headER = headER_reader('HE53_Bootstrap_metadata1',key=inST) # Find header metadata template for the appropriate instrument (ac-s or hs6)\n",
    "    ### Second, define user-supplied variables INDEPENDENT of instrument type\n",
    "    ## 2a. define keyword arguments\n",
    "    itR = str(parameters['itR']) # Define bootstrapping iteration\n",
    "    lambdA = parameters['wvl'] # Define wavelengths of IOP data\n",
    "    # 1b. define parameters from user-supplied template\n",
    "    biN = str(User_Dict['Depth Bin Size (m)']) # define depth bin\n",
    "    STN = User_Dict['Station Name'] # specify station name\n",
    "    now = dt.now() # Find current date and time\n",
    "    dt_string = now.strftime(\"%d/%m/%Y %H:%M:%S\") # Convert date and time into a time stamp (dd/mm/yyyy HH:MM:SS)\n",
    "    # 1c. update header metadata template with instrument-independent information\n",
    "    headER['1'] += dt_string  # Add current date and time to header metadata\n",
    "    headER['6'] += biN + ' m' # Add depth bin size to header metadata\n",
    "    headER['7'] += itR # bootstrapping iteration\n",
    "    ### Third, define user-supplied variables DEPENDENT on instrument type\n",
    "    # 3a. define parameters from user-supplied template\n",
    "    if inST == 'ac-s':\n",
    "        # If user is processing ac-s data...\n",
    "        pthWY = User_Dict['acs file*'] # Use the pathway of the original ac-s file as a place to store data\n",
    "        XTRA_proc = User_Dict['c-interpolation'] # are attenuation wavelengths interpolated to absorption using spine?\n",
    "        abs_str = '\\ta' # Absorption marker for column headers\n",
    "    elif inST == 'hs6': \n",
    "        # If user is processing hs6 data...\n",
    "        pthWY = User_Dict['hs6 file*'] # Use the pathway of the original ac-s file as a place to store data\n",
    "        XTRA_proc = User_Dict['Doxaran-Correction']\n",
    "        abs_str = '\\tbb' # Backscattering marker for column headers\n",
    "    # 3b. update header metadata template with instrument-DEPENDENT information supplied by user\n",
    "    headER['5'] += pthWY # add origianl IOP file (individual spectra)\n",
    "    headER['8'] += XTRA_proc # specify whether or not attenuation wavelengths were interpolate/Doxaran correction was applied\n",
    "    headER['10'] += str(len(lambdA)) # add number of wavelengths (channels) to the 11th header line\n",
    "    ### Fourth, specify the header fields and wavelengths for the IOP data\n",
    "    fieldS_c = '' # Empty character array to store attenuation column headers (fields) when bootstrapping ac-s data\n",
    "    for l in lambdA:\n",
    "        # for-loop takes each ac-s wavelength and converts it into an absorption and an attenuation column header\n",
    "        headER['9'] += abs_str + str(l) # Absorption column header\n",
    "        headER['10'] += '\\t' + str(l) # add wavelength onto header line 11\n",
    "        if inST == 'ac-s':\n",
    "            # If user is processign ac-s data...\n",
    "            fieldS_c += '\\tc' + str(l) # Create attenuation column headers\n",
    "    headER['9'] += fieldS_c # Add attenuation column headers onto header line 10\n",
    "    ### Fifth, specify name and location of the binned IOP file name \n",
    "    File_BS = STN + '_' + inST + '_bin' + biN + '_itr' + itR + '.txt' # Create file name using user-specified information\n",
    "    dash_IND = [w.start() for w in re.finditer('/', pthWY)] # Find forward slashes (specify folders) in the file pathway of unbinned IOPs\n",
    "    dash_IND = dash_IND[-1] # Locate the last forward slash\n",
    "    pthWY = pthWY[:dash_IND+1] # Separate pathway from file name\n",
    "    diRECTory = pthWY + inST + '/' # Add a new segment to the unbinned IOPs file pathway. This will become a new folder\n",
    "    createFolder(diRECTory) # Create a new folder for binned IOP file(s) (if it doesn't yet exist)\n",
    "    ### Sixth, write the bootstrapped IOP file!\n",
    "    with open(diRECTory+File_BS, 'w', newline='') as txtfile: # Write a new ascii file\n",
    "        for i in np.arange(len(headER)):\n",
    "            # for-loop writes the header metadata into the new file\n",
    "            txtfile.write(headER[str(i)] + '\\n')\n",
    "        for nn in np.arange(len(BOOTSRAP_MATRIX[:,0])):\n",
    "            # for-loop writes the binned IOP data into the new file\n",
    "            if ~np.isnan(BOOTSRAP_MATRIX[nn,1]):\n",
    "                # if the last line is not reached\n",
    "                datUM = str(BOOTSRAP_MATRIX[nn,:].tolist())[1:-1] # Convert row of data from numpy array into comma-separated string. Remove brackets\n",
    "                datUM = datUM.replace(', ','\\t') # replace commas with tabs\n",
    "                txtfile.write(datUM + '\\n') # Write the data into the file\n",
    "                keY = nn # this variable preserves last row of actual values e.g. (not nan values)\n",
    "        BOOTSRAP_MATRIX[keY,0] = -1 # After data has been written into the file, replace the first column in the last line (depth) with -1\n",
    "        datUM = str(BOOTSRAP_MATRIX[nn,:].tolist())[1:-1] # Re-convert last row of data to comma-separated string\n",
    "        datUM = datUM.replace(', ','\\t') # remove end of line characters\n",
    "        txtfile.write(datUM + '\\n') # re-write the last line of code into the file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def IOP_Bootstrapper(IOP_MATRIX,User_Dict,**parameters):\n",
    "    \"\"\" Performs multiple single bootstrapping (with replacement) iteration for depth-binning IOPs. Compatible with ac-s and hs6.\n",
    "    Inputs:\n",
    "        IOP_MATRIX - first column MUST be depth\n",
    "        User_Dict - python dictionary of user-supplied information (see UserTemplate)\n",
    "    Keyword arguments:\n",
    "        INST - instrument type acs/hs6\n",
    "        wvl - numpy array of wavelengths\n",
    "    Output:\n",
    "        Boostrap_FileWriter is nested in this function. Therefore, the outputs are hydrolight-compatible IOP files\"\"\"\n",
    "    INST = parameters['INST']  # Specify instrument used to collect data: acs or hs6\n",
    "    lambdA = parameters['wvl'] # Specify wavelengths (used to write the file)\n",
    "    BIN = float(User_Dict['Depth Bin Size (m)']) # Define user-input bin size for bootstrapping\n",
    "    n = int(User_Dict['Bootstrap Iterations']) # Number of bootstrap iterations to be carried out\n",
    "    Depth = IOP_MATRIX[:,0] # Find the depths for IOP spectra\n",
    "    IOP_MATRIX = IOP_MATRIX[:,1:] # Eliminate depth column from IOP_MATRIX\n",
    "    # Step 1: Create median depth bins for IOPs based on depth bin size\n",
    "    d_min = np.floor(min(Depth)) # minimum depth\n",
    "    d_max = np.ceil(max(Depth)) # maximum depth\n",
    "    Depth_edgE = np.arange(d_min,d_max+BIN,BIN) # Calculate edges (boundaries) of each depth bin\n",
    "    Depth_BINED = (Depth_edgE[1:]+Depth_edgE[:-1])/2 # # Calculate depth bin medians\n",
    "    BinHalf = BIN/2 # Calculate half of depth bin\n",
    "    # Step 2: Create bootstrapped ac-s files that are compatable with Hydrolight\n",
    "    for k in np.arange(n):\n",
    "        #  for-loop creates n bootstrap matrices, with 'n' corresponding to the number of \n",
    "        # bootstrapping iterations or \"runs\" that are specified by the user. Each iteration \n",
    "        # consists of one matrix of depth-binned IOPs (either absorption/attenuation or backscattering).\n",
    "        BOOTSRAP_MATRIX = np.zeros([len(Depth_BINED),len(IOP_MATRIX[0,:])+1])*np.nan \n",
    "        # Create empty (nan) matrix for depth-binned bootstrapped IOP spectra\n",
    "        BOOTSRAP_MATRIX[:,0] = Depth_BINED # Make the first column matrix binned depths\n",
    "        for h,i in enumerate(Depth_BINED):\n",
    "            # For each individual depth-bin, perform bootstrapping with replacement and calculate\n",
    "            # mean spectrum of bootstrapped IOP spectra. Repeat for each depth bin.\n",
    "            Depth_BIN_IND = np.where(np.logical_and(Depth >i-BinHalf, i+BinHalf>=Depth))\n",
    "            # Find index of each IOP spectrum that falls between the given depth bin\n",
    "            samplE_size = int(np.rint(len(Depth_BIN_IND[0]))) # determine indices of spectra that fall within the selected depth range\n",
    "            if samplE_size > 0:\n",
    "                # If there are IOP values at the specified depth range...\n",
    "                BootSTRAP_IND = np.random.uniform(low=np.min(Depth_BIN_IND), high=np.max(Depth_BIN_IND), size=(samplE_size,1)) \n",
    "                # (above) randomly select numbers in range of the selected indices (these are ordered according to depth)\n",
    "                BootSTRAP_IND = np.rint(BootSTRAP_IND) # Convert numbers to integers to correspond with indices\n",
    "                BOOTSRAP_TEMP = np.zeros([len(BootSTRAP_IND),len(IOP_MATRIX[0,:])]) \n",
    "                # Create an empty matrix of zeros (above) that can fit spectra (indices x wavelengths)\n",
    "                for keY,j in enumerate(BootSTRAP_IND):\n",
    "                    # This for-loop will cycle through the matrix of zeros and assign a randomly-selected IOP spectra (e.g. indices) to each row \n",
    "                    BOOTSRAP_TEMP[keY,:] = IOP_MATRIX[int(j),:] # Assign randomly-selected IOP spectrum to a row\n",
    "                BOOTSRAP_MATRIX[h,1:] = np.mean(BOOTSRAP_TEMP,axis=0) # Take mean of randomly-selected IOP spectra\n",
    "        BOOTSRAP_MATRIX = np.around(BOOTSRAP_MATRIX,decimals=6) # Round all depth-binned IOPs to six decimal places\n",
    "        Boostrap_FileWriter(BOOTSRAP_MATRIX,User_Dict,itR=k,INSTRUMENT=INST,wvl=lambdA) \n",
    "        # Write bootstrapped (depth-binned) IOP matrix into a Hydrolight-compatible data file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [],
   "source": [
    "def HE53_batch_WRITER(User_Dict):\n",
    "    \"\"\"Creates a Hydrolight-compatible ascii file using depth-binned IOP data. File is output in a separate folder.\n",
    "    This folder is placed in the same directory as the original file containing non-binned IOP data.\n",
    "    Inputs:\n",
    "        BOOTSRAP_MATRIX - Depth-binned IOP matrix\n",
    "        User_Dict - Python dictionary containing user-specified inputs\n",
    "    keyword arguments:\n",
    "        itR - specify bootstrap iteration. Must be an integer.\n",
    "        wvl - numpy array of IOP wavelengths.\n",
    "        INSTRUMENT - specify whether IOP matrix is 'acs' or 'hs6'\n",
    "    Outputs:\n",
    "        Series of Hyrolight batch files that are ready to process\"\"\"\n",
    "    ### First, import batch file template as batch_DATA.\n",
    "    batch_DATA = headER_reader('HE53_Bootstrap_metadata1',key='batch')\n",
    "    ### Second, assign station-specific information to the batch file template.\n",
    "    batch_DATA['2'] += User_Dict['Station Name'] + '_bin' + User_Dict['Depth Bin Size (m)'] + '_itr' # Construct root file. Add to 3rd line\n",
    "    batch_DATA['6'] += User_Dict['Chl Conc (mg/L)'] + ',' # Assign chlorophyll concentration to the 7th header line\n",
    "    ##2a. Compute the hour and minute in decimal form\n",
    "    timE = User_Dict['GMT (HH'].split(':') # Define time stamp (as string)\n",
    "    hr = int(timE[2]) # Find the hour (GMT)\n",
    "    hrFRACT = (int(timE[3])*60 + int(timE[4]))/3600 # Convert minutes and seconds into \"hour\" (decimal)\n",
    "    batch_DATA['a1'] += str(hr+hrFRACT) + ', 0, 0' # Add decimal hour into line a1 of header template\n",
    "    ## 2b. Compute day of year\n",
    "    d = str(dt.strptime(User_Dict['Date (mm-dd-yyyy)'],\"%m-%d-%Y\").timetuple().tm_yday) # Calculate day of year from time stamp using datetime    \n",
    "    ## 2b (continued). Add day of year, latitude, lontitude, atmospheric pressure, \"5\", % Humidity, Precipitable Water Content (cm), horizontal visibility, and windspeed to line a2 of header template line a3\n",
    "    batch_DATA['a2'] += d + ', ' + User_Dict['Latitude****'] + ', ' + User_Dict['Longitude****'] + ', ' + User_Dict[ 'Pressure (inHg)'] + ', 5, ' + User_Dict['Humidity (%)'] + ', ' + User_Dict['Precipitable Water Content (cm)'] + ', ' + User_Dict['Horizontal Visibility (km)'] + ', ' + User_Dict['Wind Speed (m/s)'] + ', -99' \n",
    "    ## 2c. Add windspeed, water temperature, and salinity to header template a4 (below)\n",
    "    batch_DATA['a3'] += User_Dict['Wind Speed (m/s)'] + ', -1.34, ' + User_Dict['Water Temperature (deg. C)'] +  ', ' + User_Dict['Salinity (PSU)']\n",
    "    ## 2d. Add ac-s and hs6 root files (not iteration specific) to header template (lines b3 & b5 respectively). Include file pathways\n",
    "    BIN = float(User_Dict['Depth Bin Size (m)'])  # Define depth bin size\n",
    "    batch_DATA['b2'] += User_Dict['acs pathway**'] + 'acs\\\\' + User_Dict['Station Name'] + '_' + 'ac-s_bin' + str(BIN) + '_itr' # ac-s root file & pathway\n",
    "    batch_DATA['b4'] += User_Dict['hs6 pathway**'] + 'hs6\\\\' + User_Dict['Station Name'] + '_' + 'hs6_bin' + str(BIN) + '_itr' # hs6 root file & pathway\n",
    "    batch_DATA['b10'] += User_Dict['acs pathway**'] + User_Dict['Solar Irradiance file'] # Add solar irradiance file & pathway to header template\n",
    "    ### Third, determine user-supplied wavelengths put them into lines (of ten) and add them to the header template\n",
    "    wvl = User_Dict['Wavelengths***'] # Define user-supllied wavelengths\n",
    "    ## 3a. Convert string to a numpy array of wavelengths\n",
    "    if ':' in wvl:\n",
    "        # If user supplied wavelengths in the format \"minimum:increment:maximum\"\n",
    "        miN_wvl,miD_wvl,maX_wvl = wvl.split(':') # Split the string into three variables: minumum wavelength, increment, & masimum wavelength\n",
    "        miN_wvl = float(miN_wvl) # Convert strings into floats (minimum)\n",
    "        miD_wvl = float(miD_wvl) # wavelength increment\n",
    "        maX_wvl = float(maX_wvl) # maximum wavelength\n",
    "        lambdA = np.arange(miN_wvl,maX_wvl+miD_wvl,miD_wvl) # Convert into numpy array\n",
    "    else:\n",
    "        # If user supplied wavelengths as a list separated by commas\n",
    "        lambdA = wvl.split(',') # split wavelengths into a list\n",
    "        lambdA = np.asarray(lambdA,dtype=np.float) # Convert list into numpy array\n",
    "    lambdA = np.sort(lambdA) # Sort wavelengths to be sure that they are in ascending order\n",
    "    batch_DATA['23'] += str(len(lambdA)-1) # Tabulate number of wavelengths and add to 24th header template line\n",
    "    ## 3b. Print order wavelengths into rows of 10 and write them into the header template. \n",
    "    ## Header template lines contain a 10 wavelengths or fewer. Wavelengths will be added between header template lines 23 and a1\n",
    "    keY = 23 # Index variable for batch_DATA dictionary (with which to create new elements)\n",
    "    for i,l in enumerate(lambdA):\n",
    "        # for-loop adds wavelengths to template header one by one. After adding 10 wavelengths, a new line is started\n",
    "        if np.remainder(i,10) == 0:\n",
    "            # if 10 wavelengths are added to header template, start a new line of wavelengths\n",
    "            keY += 1 # Increase index variable by 1\n",
    "            batch_DATA[str(keY)] = '' # Add a new dictionary element to batch_DATA\n",
    "        batch_DATA[str(keY)] += str(l) + ', ' # Add new wavelength and comma to dictionary element\n",
    "    ### Fourth, determine user-supplied depth bins and add them to template header\n",
    "    # 4a. Determine user-specified depth bins for Hydrolight to complete radiative transfer \n",
    "    # This will be the same value as used for IOP bootstrapping.\n",
    "    max_Depth = float(User_Dict['Max Depth (m)']) # Define maximum depth\n",
    "    Depth_BINS = np.arange(0,max_Depth+BIN,BIN) # Compute a numpy array for depth bins\n",
    "    keY = 5 # Reset indexing variable to 5\n",
    "    batch_DATA['a5'] += str(len(Depth_BINS)) + ', ' # Count all depth bins and add them to header template line a6\n",
    "    for i,d in enumerate(Depth_BINS):\n",
    "        # For-loop adds depth bins to the batch file template one by one. After adding 10 wavelength a new line is started\n",
    "        if d - int(d) == 0:\n",
    "            d = int(d)\n",
    "        if np.remainder(i,10) == 0 and i != 0:\n",
    "            # if 10 depth bins are added to header template, start a new line of depth bins\n",
    "            keY += 1 # Increase indes variable by one\n",
    "            batch_DATA['a' + str(keY)] = '' # Add a new dictionary element to batch_Data\n",
    "        batch_DATA['a' + str(keY)] += str(d) + ', ' # Add new depth bin and comma to dictionary element\n",
    "    ### Fifth, now that station-specific information has been added to batch file header template, generate a batch for \n",
    "    ### each bootstrap iteration\n",
    "    ## 5a. Create a separate folder for batch files\n",
    "    dash_IND = [w.start() for w in re.finditer('/', User_Dict['acs file*'])] # Find \"/'s\" for ac-s file string\n",
    "    batch_pthwy = User_Dict['acs file*'][:dash_IND[-1]] + '/batch/' # Remove ac-s file from string and add '/batch/' as new folder\n",
    "    createFolder(batch_pthwy) # if a folder named \"batch\" doesn't exist next to ac-s file, create one\n",
    "    ## 5b. Create Hydrolight batch files for each set of bootstrapped IOP files\n",
    "    n = int(User_Dict['Bootstrap Iterations']) # determine the number of bootstrapping iterations specified by user\n",
    "    lettER = ['','a','b'] # This list of characters allows the program to cycle through different dictionary keys. \n",
    "    # These keys are denoted as blank, \"a\", and \"b\", depending on where in the batch file batch_DATA elements will be placed\n",
    "    for i in np.arange(n):  \n",
    "        # for-loop creates a new batch file each time it repeats itself. Within each for-loop iteration is a while loop\n",
    "        # which is responsible for writing individual files. This while loop requires two indices (defined directly below)\n",
    "        keY = 0 # Set numerical indexing component of batch_DATA key equal to zero\n",
    "        L = 0   # Set letter index of batch_DATA lettER key equal to zero\n",
    "        with open(batch_pthwy + 'I' + batch_DATA['2'] + str(i) + '.txt', 'w', newline='') as batchFID: \n",
    "            # Write a new batch file using iteration-specific root file. Add \"I\" at the beginning \n",
    "            # of the file name per Hydrolight protocol\n",
    "            while 1:\n",
    "                # While-loop runs through batch file header template line by line and writes them into the batch file in the correct order\n",
    "                if L > 2:\n",
    "                    # If all header template lines have been written into the batch file ...\n",
    "                    break # Break the while-loop, close batch file being written, and begin a new batch file\n",
    "                else:\n",
    "                    # If not all header template lines have been written into the batch file ...\n",
    "                    # Do the following....     \n",
    "                    try: \n",
    "                        # If combined keY properly detects batch_DATA element\n",
    "                        if lettER[L] + str(keY) == '2' or lettER[L] + str(keY) == 'b2' or lettER[L] + str(keY) == 'b4':\n",
    "                            # If the header line contains a file that changes with each bootstrap iteration...\n",
    "                            batchFID.write(batch_DATA[lettER[L] + str(keY)] + str(i) + '.txt\\n') # Write a header line into batch file\n",
    "                        else:\n",
    "                            # If the header line contains a file that is the same for each bootstrap iteration...\n",
    "                            batchFID.write(batch_DATA[lettER[L] + str(keY)] + '\\n') # Write a header line into batch file\n",
    "                        keY += 1 # increase numerical indexing component by 1 to allow for the next line to be written \n",
    "                    except:                            \n",
    "                        # If the batch_DATA does not work for eather conditional statements...\n",
    "                        keY = 0 # Reset numerical indexing component to zero\n",
    "                        L += 1  # Increase letter index by 1. Move on to the next section of elements.\n",
    "                        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Execute Hydrolight53_Bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Read user-supplied template file and extract all relevant data\n",
    "template_file = filedialog.askopenfilename() # Select template file\n",
    "User_Dict = UserTemplate(template_file) # Place template information into python dictionary\n",
    "acs_file = User_Dict['acs file*'] # Define file and pathway for ac-s (single spectra)\n",
    "hs6_file = User_Dict['hs6 file*'] # Define file and pathway for hs6 (single spectra)\n",
    "\n",
    "# Step 2: Prepare ac-s data for bootstrapping by reading it into variables containing\n",
    "# absorption (Total_Abs), attenuation (Total_Atn), depth (Depth), and wavelength (lambdA)\n",
    "if 'mat' in acs_file[-4:-1]:\n",
    "    # If user specified a hdf5/mat file...\n",
    "    Total_Abs,Total_Atn,Depth,lambdA = HDF_LOADER(acs_file,INST='ac-s')\n",
    "else:\n",
    "    # If user specified a SEABASS ascii file...\n",
    "    Total_Abs,Total_Atn,Depth,lambdA = SEABASS_TXT_LOADER(acs_file)\n",
    "\n",
    "# Step 3: Prepare ac-s data for bootstrapping by ordering absorption and attenuation spectra \n",
    "# according to ascending depth\n",
    "D_IND = np.argsort(Depth) # Find indices for ascending depths\n",
    "Depth = Depth[D_IND]# Order depth according to ascending depth indices\n",
    "Total_Abs = Total_Abs[D_IND,:] # Order absorption spectra according to ascending depth indices\n",
    "Total_Atn = Total_Atn[D_IND,:] # Order attenuation spectra according to ascending depth indices\n",
    "\n",
    "# Step 4: Combine absorption and attenuation spectra into a single matrix\n",
    "l,w = Total_Abs.shape # Find dimensions of absorption and attenuation matrices (they should be the same size)\n",
    "IOP_MATRIX = np.zeros([l,(w*2)+1])*np.nan # Create a matrix of nan values large enough to combine absorption and attenuation (horizontally)\n",
    "IOP_MATRIX[:,0] = Depth # Insert depth as the first column of large matrix (IOP_MATRIX)\n",
    "IOP_MATRIX[:,1:w+1] = Total_Abs # Insert absorption into the matrix \n",
    "IOP_MATRIX[:,w+1:] = Total_Atn # Insert attenuation into the matrix\n",
    "\n",
    "# Step 5: Perform bootstapping with replacement on single ac spectra according to user-specified depth bin size. \n",
    "IOP_Bootstrapper(IOP_MATRIX,User_Dict,INST='ac-s',wvl=lambdA) # Create Hydrolight-compatible, depth-binned ac-s files for each bootstrap\n",
    "\n",
    "# Step 6: Prepare hs6 data for bootstrapping by reading it into variables containing\n",
    "# particulate backscattering (Total_hs), depth (Depth_hs), and wavelength (lambdA_hs)\n",
    "if 'mat' in hs6_file[-4:-1]:\n",
    "    # If user specified a hdf5/mat file...\n",
    "    Total_hs,Depth_hs,lambdA_hs = HDF_LOADER(hs6_file,INST = 'hs6')\n",
    "else:\n",
    "    # If user specified a SEABASS ascii file...\n",
    "    Total_hs,blank,Depth_hs,lambdA_hs = SEABASS_TXT_LOADER(hs6_file)\n",
    "\n",
    "# Step 7: Prepare hs6 data for bootstrapping by ordering particulate baskcattering spectra according to ascending depth\n",
    "D_IND_hs = np.argsort(Depth_hs) # Find indices for ascending depths\n",
    "Depth_hs = Depth_hs[D_IND_hs] # Order depth according to ascending depth indices\n",
    "Total_hs = Total_hs[D_IND_hs,:] # Order particulate backscattering spectra according to ascending depth indices\n",
    "\n",
    "# Step 8: Prepare hs6 data for bootstrapping by ordering particulate backscattering spectra according to ascending wavlength\n",
    "L_IND = np.argsort(lambdA_hs) # Find indices for ascending wavelengths\n",
    "lambdA_hs = lambdA_hs[L_IND] # Order wavelengths according to ascending wavelength indices\n",
    "Total_hs = Total_hs[:,L_IND] # Order particulate backscattering according to ascending wavelength indices\n",
    "\n",
    "# Step 9: Combine particulate backscattering specra and depths into a matrix\n",
    "l_hs, w_hs = Total_hs.shape # Find shape of particulate backscattering mtarix\n",
    "IOP_MATRIX_hs = np.zeros([l_hs,w_hs+1])*np.nan # Create nan matrix for backscattering matrix with one additional column for depth\n",
    "IOP_MATRIX_hs[:,0] = Depth_hs # Insert depth as the first column of large matrix\n",
    "IOP_MATRIX_hs[:,1:w+1] = Total_hs # Insert particulate backscttering in the remaining columns of the nan matrix\n",
    "\n",
    "# Step 10: Perform bootstapping with replacement on single ac spectra according to user-specified depth bin size. \n",
    "IOP_Bootstrapper(IOP_MATRIX_hs,User_Dict,INST='hs6',wvl=lambdA_hs) # Create Hydrolight-compatible, depth-binned ac-s files for each bootstrap\n",
    "\n",
    "# Step 11: Create Hydrolight batch files consistent with instructions\n",
    "HE53_batch_WRITER(User_Dict)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}