02_chamber_fluxes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Laser chambers (Jonathan D. Müller)\n",
    "\n",
    "Loads the raw data of soil/branch chambers, applies corrections and calculates fluxes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialization parameters\n",
    "\n",
    "# Data input\n",
    "project_path = './'\n",
    "\n",
    "project_path_amb   = project_path + '02_preprocessed_data/ambient/'\n",
    "project_path_laser = project_path + '02_preprocessed_data/laser computer/'\n",
    "project_path_par   = project_path + '02_preprocessed_data/PAR/'\n",
    "project_path_tc    = project_path + '02_preprocessed_data/thermocouples/'\n",
    "project_path_flow  = project_path + '02_preprocessed_data/flow/'\n",
    "\n",
    "# Additional data, logbook etc.\n",
    "project_path_additional = project_path + '00_additional_data/'\n",
    "cos_timing_fn       = 'COS_timing.xlsx'\n",
    "chambers_sensors_fn = 'COS_chambers_sensor_logbook.xlsx'\n",
    "\n",
    "# Output\n",
    "project_path_output = project_path + '03_processed_data/'\n",
    "output_fn = 'laser'\n",
    "\n",
    "# List of months to process\n",
    "# - If empty, all available data is processed, e.g. [  ]\n",
    "# - Otherwise, specify using a string of year-month, e.g. ['2018-06'] or ['2017-03','2018-06']\n",
    "month_list = ['2022-01']\n",
    "\n",
    "# Main parameters\n",
    "par_threshold    = 10     # Threshold by which to determine whether to keep PAR\n",
    "averaging_period = 30     # How many [s] of each measurement are used to average the data. Currently given by the COS_timing file, but here can be shortened to <25s\n",
    "keep_all_raw     = False   # If true, keep all data given by COS_timing when creating the raw output. Otherwise, keep the data of the averaging periody only\n",
    "#dew_threshold    = 98     # Percentage of RH when dew starts to form, remove data when this is the case\n",
    "\n",
    "# Choose how to find the ambient value for each measurement. There are 2 possibilities:\n",
    "# - interpolated: This interpolates ambient measurements and adds them\n",
    "# - next:         This simply finds the next ambient and adds it\n",
    "#find_ambient     = 'next'\n",
    "find_ambient     = 'interpolated'\n",
    "\n",
    "# Leaf area parameters\n",
    "proj_la_factor = np.pi # Factor used to convert total leaf area to projected leaf area\n",
    "\n",
    "# leaf area index LAI\n",
    "LAI_ctr = 1.5 # control\n",
    "LAI_irr = 2.0 # irrigation\n",
    "\n",
    "# Pump parameters: Values different from these will be deleted\n",
    "min_abs_flow_branch = 3       # Minimum allowed pump flow for branch chambers\n",
    "min_abs_flow_soil   = 0.5     # Minimum allowed pump flow for soil chambers\n",
    "#max_delta_flow = 4     # Maximum difference between ambient and sample pumps\n",
    "#max_stddev_flow = 0.5  # Maximum standard deviation within pump flows\n",
    "\n",
    "# Calibration factor for gases\n",
    "# - Factors by which the gases have to be corrected (Measured for our laser by Rafat Qubaja, Spring 2021)\n",
    "fact_h2o = 1/1           # Calibration measured\n",
    "fact_co2 = 396.5/378.52  # Calibration measured\n",
    "fact_cos = 514/583.52    # Calibration measured\n",
    "fact_co  = 1/1           # Not measured\n",
    "\n",
    "# Day/Night PAR limit\n",
    "# - Min. value of PAR when it's still day, using ambient PAR measurement\n",
    "day_par_min = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Load additionally needed data\n",
    "- Leaf area logbook\n",
    "- Chambers status (checks if the chamber was working or not)\n",
    "- Timing information (when the laser was measuring what)\n",
    "- PAR information (to which datalogger channel connected and if it's working)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read additional data\n",
    "#----------------------\n",
    "\n",
    "# Read laser timing (lt) information, used to extract data when the chamber was closed\n",
    "lt = pd.read_excel(project_path_additional + cos_timing_fn, sheet_name='timing')\n",
    "lt.rename(columns={'date':'timestamp'}, inplace=True)\n",
    "lt.drop(['comment'], axis=1, inplace=True)\n",
    "lt['timestamp'] = pd.to_datetime( lt.timestamp, format='%Y-%m-%d', utc=True, errors='raise')\n",
    "lt.dropna(subset=['timestamp'], inplace=True)\n",
    "#display(lt)\n",
    "\n",
    "# Read chamber status (cs)\n",
    "cs = pd.read_excel(project_path_additional + chambers_sensors_fn, sheet_name='chambers')\n",
    "cs.dropna(subset=['year'], inplace=True)\n",
    "cs.drop(['comment'], axis=1, inplace=True)\n",
    "# Convert dates\n",
    "cs['timestamp'] = pd.to_datetime( cs.year.map(str) + '-' + cs.month.map(str) + '-' + cs.day.map(str), format='%Y-%m-%d', utc=True, errors='raise')\n",
    "#cs['timestamp'] = pd.to_datetime( cs.year.map(str) + '-' + cs.month.map(str) + '-' + cs.day.map(str), format='%Y.0-%m.0-%d.0', utc=True, errors='raise')\n",
    "# Filter out errors (too many days in a month)\n",
    "cs = cs[cs.timestamp.notnull()]\n",
    "#display(cs)\n",
    "\n",
    "# Read PAR sensor installation and information\n",
    "par_log = pd.read_excel(project_path_additional + chambers_sensors_fn, sheet_name='par')\n",
    "par_log.dropna(subset=['date'], inplace=True)\n",
    "par_log.rename(columns={'date':'timestamp'}, inplace=True)\n",
    "par_log.drop(['comment'], axis=1, inplace=True)\n",
    "#display(par_log)\n",
    "\n",
    "# Log of flowmeter MUX channels in datalogger\n",
    "fm_log = pd.read_excel(project_path_additional + chambers_sensors_fn, sheet_name='flow')\n",
    "fm_log.dropna(subset=['date'], inplace=True)\n",
    "fm_log.rename(columns={'date':'timestamp'}, inplace=True)\n",
    "fm_log.drop(['comment'], axis=1, inplace=True)\n",
    "#display(fm_log)\n",
    "\n",
    "# Log of thermocouple MUX channels in datalogger\n",
    "tc_log = pd.read_excel(project_path_additional + chambers_sensors_fn, sheet_name='tc')\n",
    "tc_log.dropna(subset=['date'], inplace=True)\n",
    "tc_log.rename(columns={'date':'timestamp'}, inplace=True)\n",
    "tc_log.drop(['comment'], axis=1, inplace=True)\n",
    "#display(tc_log)\n",
    "\n",
    "# Log of surface area (branches and soil) in the chambers\n",
    "sa_log = pd.read_excel(project_path_additional + chambers_sensors_fn, sheet_name='surf_area')\n",
    "sa_log.dropna(subset=['date'], inplace=True)\n",
    "sa_log.rename(columns={'date':'timestamp'}, inplace=True)\n",
    "sa_log.drop(['comment'], axis=1, inplace=True)\n",
    "#display(sa_log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Debugging\n",
    "#display(par_log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Physics: Unit converstions, fluxes, etc.\n",
    "#-----------------------------------------\n",
    "\n",
    "# Calculate the dewpoint temperature in C\n",
    "def calculate_dewpointC(T_C, RH):\n",
    "    dp = 243.04*(np.log(RH/100)+((17.625*T_C)/(243.04+T_C)))/(17.625-log(RH/100)-((17.625*T_C)/(243.04+T_C)))\n",
    "    return(dp)\n",
    "\n",
    "# Calculate saturation vapour pressure from pressure and temperature\n",
    "# - 2 methods are available. Jones uses air pressure, Campbell & Norman do not\n",
    "def calculate_es(T_C, P_Pa):\n",
    "    # Jones p.348 (appendix 4)\n",
    "    #es = (1.00072+(10**(-7)*P_Pa*(0.032+5.9*10**(-6)*T_C**2))) * (611.21*np.exp( (18.678-(T_C/234.5))*T_C/(257.14+T_C) ))\n",
    "\n",
    "    # Eddypro manual: https://www.licor.com/env/support/EddyPro/topics/calculate-micromet-variables.html\n",
    "    # Campbell & Norman (1998)\n",
    "    T_K = T_C + 273.15\n",
    "    es = T_K**(-8.2) * np.exp(77.345 + 0.0057*T_K - 7235 * T_K**(-1))\n",
    "    return(es)\n",
    "\n",
    "# Calculates VPD from different environmental variables\n",
    "def calculate_VPD(T_C, h2o_mmol_mol, P_Pa):\n",
    "    # Unit conversions \n",
    "    T_K = T_C + 273.15           # Temperature in K\n",
    "    h2o_mol_mol = h2o_mmol_mol * 10**(-3) # water in [mol mol-1]\n",
    "\n",
    "    # From Eddypro manual: https://www.licor.com/env/support/EddyPro/topics/calculate-micromet-variables.html\n",
    "    e = h2o_mol_mol * P_Pa         # Water vapor partial pressure (Pa)\n",
    "    es = calculate_es(T_C, P_Pa)   # Water vapor partial pressure at saturation (Pa)\n",
    "    VPD = es - e                   # VPD (Pa)\n",
    "    return(VPD)\n",
    "\n",
    "# Converts water concentration [mmol mol] to RH [%]\n",
    "def convert_mmol_RH(T_C, h2o_mmol_mol, P_Pa):\n",
    "    # Unit conversions \n",
    "    T_K = T_C + 273.15           # Temperature in K\n",
    "    h2o_mol_mol = h2o_mmol_mol * 10**(-3) # water in [mol mol-1]\n",
    "    \n",
    "    #es = calculate_es(T_C, P_Pa)\n",
    "    #RH <- 0.263*P_Pa*((h2o_mmol_mol*18.02/28.97)/1000)*np.exp(17.67*(T_C)/(T_K-29.65))**(-1)\n",
    "    #RH = 100 if (RH > 100) else RH\n",
    "    #RH = np.nan if (RH < 5) else RH\n",
    "\n",
    "    # From Eddypro manual: https://www.licor.com/env/support/EddyPro/topics/calculate-micromet-variables.html\n",
    "    e = h2o_mol_mol * P_Pa         # Water vapor partial pressure (Pa)\n",
    "    es = calculate_es(T_C, P_Pa)   # Water vapor partial pressure at saturation (Pa)\n",
    "    RH = e/es * 100                # RH (%)\n",
    "    return(RH)\n",
    "\n",
    "# Density of dry air\n",
    "# - https://www.licor.com/env/support/EddyPro/topics/calculate-micromet-variables.html\n",
    "def calculate_rho_dry_air(T_C, h2o_mmol_mol, P_Pa):\n",
    "    # Unit conversions \n",
    "    h2o_mol_mol = h2o_mmol_mol * 10**(-3) # water in [mol mol-1]\n",
    "    \n",
    "    # Constants\n",
    "    R_dry_air = 287.058     # [J kg-1 K-1] Specific gas const dry air\n",
    "    T_K = T_C + 273.15\n",
    "    \n",
    "    # Preparations\n",
    "    R     = 8.314463             # Ideal gas constant (J K-1 mol-1)\n",
    "    M_d   = 0.02897              # molecular weights of dry air (kg mol-1)\n",
    "    M_h2o = 0.01802              # molecular weights of water vapour (kg mol-1)\n",
    "    e = h2o_mol_mol * P_Pa       # Water vapor partial pressure (Pa)\n",
    "    P_d = P_Pa - e               # Dry air partial pressure (P_d, P_a)\n",
    "    \n",
    "    rho_dry_air = P_d / (R_dry_air * T_K) # Density of dry air (use for approximation)\n",
    "    return(rho_dry_air)\n",
    "\n",
    "# Density of moist air\n",
    "def calculate_rho_moist_air(T_C, h2o_mmol_mol, P_Pa):\n",
    "    # Unit conversions \n",
    "    T_K = T_C + 273.15           # Temperature in K\n",
    "    h2o_mol_mol = h2o_mmol_mol * 10**(-3) # water in [mol mol-1]\n",
    "\n",
    "    # Preparations\n",
    "    R     = 8.314463             # Ideal gas constant (J K-1 mol-1)\n",
    "    M_d   = 0.02897              # molecular weights of dry air (kg mol-1)\n",
    "    M_h2o = 0.01802              # molecular weights of water vapour (kg mol-1)\n",
    "    e = h2o_mol_mol * P_Pa       # Water vapor partial pressure (Pa)\n",
    "    P_d = P_Pa - e               # Dry air partial pressure (P_d, P_a)\n",
    "    rho_d = P_d / (R / M_d * T_K) # Dry air mass density (rho_d, kg m-3)\n",
    "    v_d = M_d / rho_d            # Dry air molar volume (vd, m3 mol-1)\n",
    "    v_a = v_d * P_d/P_Pa         # Air molar volume (vd, m3mol-1) \n",
    "    rho_h2o = h2o_mol_mol * M_h2o / v_a # Ambient water vapor mass density (kg m-3)\n",
    "\n",
    "    # Moist air mass density (ρa, kg m-3) \n",
    "    rho_air = rho_d + rho_h2o\n",
    "\n",
    "    return(rho_air)\n",
    "\n",
    "# Dry air heat capacity at constant pressure\n",
    "# cp_d in [J kg-1 K-1]\n",
    " # https://www.licor.com/env/support/EddyPro/topics/calculate-micromet-variables.html\n",
    "def calculate_cp_dry_air(T_C):\n",
    "    cp = 1005 + ((T_C + 23.12)**2)/3364\n",
    "    return(cp)\n",
    "\n",
    "# Specific heat capacity of moist air at constant pressure\n",
    "# cp_m in [J kg-1 K-1]\n",
    "# https://www.licor.com/env/support/EddyPro/topics/calculate-micromet-variables.html\n",
    "def calculate_cp_moist_air(T_C, h2o_mmol_mol, P_Pa):\n",
    "    # Unit conversions \n",
    "    T_K = T_C + 273.15           # Temperature in K\n",
    "    h2o_mol_mol = h2o_mmol_mol * 10**(-3) # water in [mol mol-1]\n",
    "\n",
    "    # RH\n",
    "    RH = convert_mmol_RH(T_C, h2o_mmol_mol, P_Pa)\n",
    "\n",
    "    # Water vapor heat capacity at constant pressure (cp_h2o, J kg-1 K-1)\n",
    "    cp_h2o = 1859 + 0.13*RH + (0.193 + 5.6*10**(-3) * RH)*T_C + (10**(-3) + 5 * 10**(-5)*RH)*T_C**2\n",
    "\n",
    "    # Preparations\n",
    "    R     = 8.314463             # Ideal gas constant (J K-1 mol-1)\n",
    "    M_d   = 0.02897              # molecular weights of dry air (kg mol-1)\n",
    "    M_h2o = 0.01802              # molecular weights of water vapour (kg mol-1)\n",
    "    e = h2o_mol_mol * P_Pa       # Water vapor partial pressure (Pa)\n",
    "    P_d = P_Pa - e               # Dry air partial pressure (P_d, P_a)\n",
    "    rho_d = P_d / (R / M_d * T_K) # Dry air mass density (rho_d, kg m-3)\n",
    "    v_d = M_d / rho_d            # Dry air molar volume (vd, m3 mol-1)\n",
    "    v_a = v_d * P_d/P_Pa         # Air molar volume (vd, m3mol-1) \n",
    "    rho_h2o = h2o_mol_mol * M_h2o / v_a # Ambient water vapor mass density (kg m-3)\n",
    "\n",
    "    # Moist air mass density (ρa, kg m-3) \n",
    "    rho_air = rho_d + rho_h2o\n",
    "\n",
    "    # Specific humidity (Q, kg kg-1) \n",
    "    Q = rho_h2o / rho_air\n",
    "\n",
    "    # cp_moist\n",
    "    cp = calculate_cp_dry_air(T_C) * (1-Q) + cp_h2o * Q\n",
    "    return(cp)\n",
    "\n",
    "# Calculates the flux of water\n",
    "def calculate_h2o_flux(T_C, P_Pa, h2o_mmol_mol_ambient, h2o_mmol_mol_chamber, airflow_lpm, area_m2):\n",
    "    # Unit conversions\n",
    "    T_K = T_C + 273.15 # Temperature in K\n",
    "    h2o_mol_mol_ambient = h2o_mmol_mol_ambient * 10**(-3)\n",
    "    h2o_mol_mol_chamber = h2o_mmol_mol_chamber * 10**(-3)\n",
    "    \n",
    "    # Preparations\n",
    "    R     = 8.314463             # Ideal gas constant (J K-1 mol-1)\n",
    "    M_d   = 0.02897              # molecular weights of dry air (kg mol-1)\n",
    "    M_h2o = 0.01802              # molecular weights of water vapour (kg mol-1)\n",
    "    e = h2o_mol_mol_chamber * P_Pa # Water vapor partial pressure (Pa)\n",
    "    P_d = P_Pa - e               # Dry air partial pressure (P_d, P_a)\n",
    "    rho_d = P_d / (R / M_d * T_K) # Dry air mass density (rho_d, kg m-3)\n",
    "    v_d = M_d / rho_d            # Dry air molar volume (v_d, m3 mol-1)\n",
    "    v_a = v_d * P_d/P_Pa         # Air molar volume (v_a, m3mol-1)\n",
    "    \n",
    "    # Convert airflow from LPM to mol.s-1\n",
    "    airflow_m3_s = airflow_lpm / 1000 / 60 # airflow conversion to total flow [m3.s-1]\n",
    "    airflow_mol_s = airflow_m3_s / v_a     # airflow conversion to moist air [mol.s-1]\n",
    "    \n",
    "    # Note: We are measuring the outflow from the chamber, not the inflow (as opposed to the LI-6400)\n",
    "    #       Therefore, we need to change eq. 1-2 (p. 1-7) to be ue = uo - sE\n",
    "    #       This leads to eq. 1-3 changing to sE = uo*wo - (uo - sE)*we\n",
    "    #       Finally, eq. 1-4 becomes E = uo*(wo - we) / (s*(1 - we))\n",
    "    #       Where uo = airflow_mol_s\n",
    "    #             wo = h2o_mol_mol_chamber (outgoing mol fraction)\n",
    "    #             we = h2o_mol_mol_ambient (entering mol fraction)\n",
    "    #             s  = area_m2 (leaf area)\n",
    "      \n",
    "    # LI-6400 manual, eq. 1-7\n",
    "    h2o_flux = (airflow_mol_s * (h2o_mol_mol_chamber - h2o_mol_mol_ambient)) / (area_m2 * (1 - h2o_mol_mol_ambient))\n",
    "    \n",
    "    return(h2o_flux) # mol.m-2.s-1\n",
    "\n",
    "# Inspired from LI-6400 manual: Uses the water flux because water changes the air density. If stomata are open and H2o is added, the gas is more 'diluted'\n",
    "def calculate_gas_flux(T_C, P_Pa, h2o_mmol_mol_ambient, h2o_mmol_mol_chamber, gas_mol_mol_ambient, gas_mol_mol_chamber, airflow_lpm, area_m2):\n",
    "    # Unit conversions\n",
    "    T_K = T_C + 273.15 # Temperature in K\n",
    "    h2o_mol_mol_ambient = h2o_mmol_mol_ambient * 10**(-3)\n",
    "    h2o_mol_mol_chamber = h2o_mmol_mol_chamber * 10**(-3)\n",
    "    \n",
    "    # Preparations\n",
    "    R     = 8.314463             # Ideal gas constant (J K-1 mol-1)\n",
    "    M_d   = 0.02897              # molecular weights of dry air (kg mol-1)\n",
    "    M_h2o = 0.01802              # molecular weights of water vapour (kg mol-1)\n",
    "    e = h2o_mol_mol_chamber * P_Pa # Water vapor partial pressure (Pa)\n",
    "    P_d = P_Pa - e               # Dry air partial pressure (P_d, P_a)\n",
    "    rho_d = P_d / (R / M_d * T_K) # Dry air mass density (rho_d, kg m-3)\n",
    "    v_d = M_d / rho_d            # Dry air molar volume (vd, m3 mol-1)\n",
    "    v_a = v_d * P_d/P_Pa         # Air molar volume (vd, m3mol-1)\n",
    "    \n",
    "    # Convert airflow from LPM to mol.s-1\n",
    "    airflow_m3_s = airflow_lpm / 1000 / 60 # airflow conversion to total flow [m3.s-1]\n",
    "    airflow_mol_s = airflow_m3_s / v_a     # airflow conversion to moist air [mol.s-1]\n",
    "    \n",
    "    # H2O transpiration flux [mol.m-2.s-1]\n",
    "    h2o_flux_mol_m2_s = calculate_h2o_flux(T_C, P_Pa, h2o_mmol_mol_ambient, h2o_mmol_mol_chamber, airflow_lpm, area_m2)\n",
    "    \n",
    "    # Note: We are measuring the outflow from the chamber, not the inflow (as opposed to the LI-6400)\n",
    "    #       Therefore, we need to change eq. 1-2 (p. 1-7) to be ue = uo - sE\n",
    "    #       The LI-6400 uses assimilation, which is a positive flux. Previously, we used negative numbers due to CO2 uptake\n",
    "    #       The sign should reflect the direction of the flux, i.e. emission is positive, and uptake negative, so assimilation should be negative\n",
    "    #       The following calculations will therefore change eq. 1-11 (p. 1-9) to sa = uo*co - ue*ce\n",
    "    #       This leads to eq. 1-12 (p. 1-9) changing to sa = uo*co - (uo - sE)*ce\n",
    "    #       Finally, eq. 1-13 becomes a = uo*(co-ce) / s + E*ce\n",
    "    #       Where uo = airflow_mol_s\n",
    "    #             co = co2_mol_mol_chamber (outgoing mol fraction)\n",
    "    #             ce = co2_mol_mol_ambient (entering mol fraction)\n",
    "    #             s  = area_m2 (leaf area)\n",
    "    #             E  = water flux (calculated above)\n",
    "    \n",
    "    # LI-6400 manual, eq. 1-13 (p. 1-10)\n",
    "    gas_flux = (airflow_mol_s * (gas_mol_mol_chamber - gas_mol_mol_ambient)) / area_m2 + h2o_flux_mol_m2_s * gas_mol_mol_ambient\n",
    "    \n",
    "    return(gas_flux) # mol.m-2.s-1\n",
    "\n",
    "# Conductance calculations\n",
    "def calculate_cos_stomatal_conductance_ball_berry(T_C, h2o_mmol_mol, P_Pa, f_h2o_mmol_m2_s1, f_co2_umol_m2_s1, co2_umol_mol_ambient, PAR):\n",
    "    # Note: T_C is the air temperature in C\n",
    "    \n",
    "    # Calculate some initial values\n",
    "    VPD_Pa = calculate_VPD(T_C, h2o_mmol_mol, P_Pa)\n",
    "    RH = convert_mmol_RH(T_C, h2o_mmol_mol, P_Pa)\n",
    "    \n",
    "    # Pressure in hPa\n",
    "    P_kPa   = P_Pa / 1000\n",
    "    VPD_kPa = VPD_Pa / 1000\n",
    "    \n",
    "    # Stomatal conductance\n",
    "    g_s_cos = (f_h2o_mmol_m2_s1 * P_kPa * 10 / VPD_kPa) * 2 / 10000\n",
    "    \n",
    "    # Temporarily convert to df\n",
    "    temp = pd.DataFrame({'rh': RH, 'par': PAR, 'f_co2': f_co2_umol_m2_s1, 'co2_a': co2_umol_mol_ambient, 'g_s': g_s_cos})\n",
    "    \n",
    "    # Ball-Berry correction\n",
    "    temp['bbm'] = -temp['f_co2'] * temp['rh'] * 0.01 / temp['co2_a']\n",
    "    temp.loc[(temp['rh'] > 70) & (temp['par'] > day_par_min), 'g_s'] = 17.207 * temp.loc[(temp['rh'] > 70) & (temp['par'] > day_par_min), 'bbm'] / 0.0487\n",
    "    \n",
    "    return(temp['g_s'])\n",
    "\n",
    "def calculate_cos_total_conductance(f_cos_pmol_m2_s1, cos_pmol_mol_ambient):\n",
    "    # Total conductance\n",
    "    g_t_cos = - (f_cos_pmol_m2_s1 / cos_pmol_mol_ambient)\n",
    "    return(g_t_cos)\n",
    "\n",
    "def calculate_cos_internal_conductance(g_total, g_stomatal):\n",
    "    # Total conductance\n",
    "    g_i_cos = (g_total**(-1) - g_stomatal**(-1))**(-1)\n",
    "    return(g_i_cos)\n",
    "\n",
    "# Calculates leaf/soil relative uptake, i.e. LRU and SRU\n",
    "def relative_uptake(f_cos_pmol_m2_s1, f_co2_umol_m2_s1, cos_pmol_mol_ambient, co2_umol_mol_ambient):\n",
    "    ru = f_cos_pmol_m2_s1 / f_co2_umol_m2_s1 * co2_umol_mol_ambient / cos_pmol_mol_ambient\n",
    "    return(ru)\n",
    "\n",
    "def calculate_biochemical_conductance(T_leaf, LAI):\n",
    "    # Preparations\n",
    "    R  = 8.314463 # Ideal gas constant (J K-1 mol-1)\n",
    "    E0 = 40\n",
    "    ref_T_C = 20\n",
    "    ref_T_K = ref_T_C + 273.15\n",
    "    T_leaf_K = T_leaf + 273.15\n",
    "    \n",
    "    # Biochemical conductance\n",
    "    g_ca = 0.8*0.055*LAI * np.exp((E0/R) * (1/ref_T_K - 1/T_leaf_K))\n",
    "\n",
    "    return(g_ca)\n",
    "\n",
    "def calculate_mesophyll_conductance(T_leaf, LAI):\n",
    "    g_m = 0.188*LAI * np.exp(-0.5*(np.log((T_leaf/28.8)/0.61))**2)\n",
    "\n",
    "    return(g_m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main functions\n",
    "#---------------\n",
    "\n",
    "# Delete all relevant files in a folder.\n",
    "# - Used to remove 1min monthly files before running to prevent appending to existing files\n",
    "def empty_dir(directory):\n",
    "    files = glob.glob(directory + '*')\n",
    "    for f in files:\n",
    "        month_id = f[-10:-6] + '-' + f[-6:-4]\n",
    "        if(month_id in month_list):\n",
    "            print('Remove ' + f)\n",
    "            os.remove(f)\n",
    "        if(not month_list):\n",
    "            print('Remove ' + f)\n",
    "            os.remove(f)\n",
    "    pass\n",
    "\n",
    "# Read an input file\n",
    "def read_data_file(input_fn):\n",
    "    df = pd.read_csv(input_fn, na_values=['NAN'])\n",
    "    df['timestamp'] = pd.to_datetime( df.timestamp, format='%Y-%m-%d %H:%M:%S', utc=True, errors='raise')\n",
    "    df['dayid'] = df['timestamp'].apply(lambda x:(x.year*10000 + x.month*100 + x.day))\n",
    "    return(df)\n",
    "\n",
    "# This finds the file corresponding to a specific date and loads it\n",
    "def find_and_load_file(fn_list, date_idx):\n",
    "    data_fn = [i for i in fn_list if date_idx in i][0] if len([i for i in fn_list if date_idx in i]) > 0 else None\n",
    "    if(data_fn is None):\n",
    "        temp_df = read_data_file(fn_list[0])\n",
    "        temp_df = temp_df.iloc[0:0]\n",
    "    else:\n",
    "        temp_df  = read_data_file(data_fn)\n",
    "    return(temp_df)\n",
    "\n",
    "# Assign block IDs for grouping and averaging in calculations\n",
    "def create_block_id(temp):\n",
    "    temp = temp.copy()\n",
    "    blockID = 0\n",
    "    blockIDs = []\n",
    "    # Create an ID variable, used only here\n",
    "    temp['id_var'] = temp['chamber'] + temp['status']\n",
    "    x_old = temp['id_var'].tolist()[0]\n",
    "\n",
    "    for x in temp['id_var'].tolist():\n",
    "        # If same chamber, give same blockID\n",
    "        if x == x_old:\n",
    "            blockIDs.append(blockID)\n",
    "            continue\n",
    "        # If switch to new chamber, give new blockID\n",
    "        else:\n",
    "            blockID += 1\n",
    "            blockIDs.append(blockID)\n",
    "        x_old = x\n",
    "    \n",
    "    # Drop ID variable\n",
    "    temp.drop('id_var', axis=1, inplace=True)\n",
    "    \n",
    "    # Add block ID\n",
    "    temp.loc[:,'chamber_block'] = blockIDs\n",
    "    \n",
    "    return(temp)\n",
    "\n",
    "# Calculate the standard deviation of means of means\n",
    "# - provide a list of standard deviation columns of each initial mean\n",
    "# - returns a list of standard deviations, by row\n",
    "def std_means_of_means(sds):\n",
    "        # Calculate a list of variances\n",
    "        var = []\n",
    "        for i in sds:\n",
    "            var.append(i**2)\n",
    "        # Calculate the mean\n",
    "        var_mean = np.mean(var, axis=0)\n",
    "        # Calculate the standard deviation\n",
    "        sd_mean = np.sqrt(var_mean)\n",
    "        return(sd_mean)\n",
    "\n",
    "# Make means of all gases, because there are multiple measurements. Also add standard deviations\n",
    "def fix_gas_units(temp):\n",
    "    temp = temp.copy()\n",
    "    # Convert units\n",
    "    temp['OCS.1'] = temp['OCS.1'] * 1000\n",
    "    temp['OCS.2'] = temp['OCS.2'] * 1000\n",
    "    temp['CO2.1'] = temp['CO2.1'] / 1000\n",
    "    temp['CO2.2'] = temp['CO2.2'] / 1000\n",
    "    temp['CO2.3'] = temp['CO2.3'] / 1000\n",
    "    temp['CO2.4'] = temp['CO2.4'] / 1000\n",
    "    temp['H2O.1'] = temp['H2O.1'] / 10**6\n",
    "    temp['CO.1'] = temp['CO.1'] * 1\n",
    "    # drop bad columns\n",
    "    temp.drop('CO2.3', axis=1, inplace=True)\n",
    "    temp.drop('CO2.4', axis=1, inplace=True)\n",
    "    # Rename columns\n",
    "    temp.rename(columns={'OCS.1':'conc.cos.1.pmol_mol'}, inplace=True)\n",
    "    temp.rename(columns={'OCS.2':'conc.cos.2.pmol_mol'}, inplace=True)\n",
    "    temp.rename(columns={'CO2.1':'conc.co2.1.umol_mol'}, inplace=True)\n",
    "    temp.rename(columns={'CO2.2':'conc.co2.2.umol_mol'}, inplace=True)\n",
    "    temp.rename(columns={'H2O.1':'conc.h2o.mmol_mol'}, inplace=True)\n",
    "    temp.rename(columns={'CO.1': 'conc.co.nmol_mol'}, inplace=True)\n",
    "    # Apply calibration correction\n",
    "    temp['conc.cos.1.pmol_mol'] = temp['conc.cos.1.pmol_mol'] * fact_cos\n",
    "    temp['conc.cos.2.pmol_mol'] = temp['conc.cos.2.pmol_mol'] * fact_cos\n",
    "    temp['conc.co2.1.umol_mol'] = temp['conc.co2.1.umol_mol'] * fact_co2\n",
    "    temp['conc.co2.2.umol_mol'] = temp['conc.co2.2.umol_mol'] * fact_co2\n",
    "    temp['conc.h2o.mmol_mol'] = temp['conc.h2o.mmol_mol'] * fact_h2o\n",
    "    temp['conc.co.nmol_mol'] = temp['conc.co.nmol_mol'] * fact_co\n",
    "    # Make mean of gases when multiple columns are available\n",
    "    temp['conc.cos.mean.pmol_mol'] = temp[['conc.cos.1.pmol_mol','conc.cos.2.pmol_mol']].mean(axis=1)\n",
    "    temp['conc.co2.mean.umol_mol'] = temp[['conc.co2.1.umol_mol','conc.co2.2.umol_mol']].mean(axis=1)\n",
    "    # Add standard deviation of gases when multiple columns are available\n",
    "    temp['conc.cos.mean_sd.pmol_mol'] = temp[['conc.cos.1.pmol_mol','conc.cos.2.pmol_mol']].std(axis=1)\n",
    "    temp['conc.co2.mean_sd.umol_mol'] = temp[['conc.co2.1.umol_mol','conc.co2.2.umol_mol']].std(axis=1)\n",
    "    return(temp)\n",
    "\n",
    "# Identify the previous log date compared to the current date of the data\n",
    "# This log date contains the relevant information to be read and applied\n",
    "def find_log_entry(log_timestamps, this_date):\n",
    "    previous_log_date = log_timestamps[0]\n",
    "    log_entry_date = 0\n",
    "    for i in log_timestamps:\n",
    "        current_date = pd.to_datetime(str(this_date), format='%Y%m%d')\n",
    "        # Find the entry that's in between the current date and the previous\n",
    "        if((current_date > previous_log_date) and (current_date < i)):\n",
    "            log_entry_date = previous_log_date\n",
    "        # If the current date is larger than the last entry, keep the last entry\n",
    "        if(current_date > log_timestamps[-1]):\n",
    "            log_entry_date = i\n",
    "        # If the entry is on a measurement day, make it 0 so that NAs are filled in\n",
    "        if(current_date == i):\n",
    "            log_entry_date = i\n",
    "        previous_log_date = i\n",
    "    return(log_entry_date)\n",
    "\n",
    "# Read the log, and extract the relevant data\n",
    "# Requires:\n",
    "# - temp: dataframe with the data of the sensor (e.g. PAR, thermocouples, ...)\n",
    "# - log_df: log-file dataframe\n",
    "# - this_day: Current date\n",
    "# - name: Name for the new columns (e.g., 'par')\n",
    "# - unit: Unit of the measurement (e.g., 'umol_m2_s1' or 'c')\n",
    "# - old_name: The start of the name of the column in the raw data (e.g., 'par_Avg(' or similar)\n",
    "# - calibration: Factor by which the sensor needs to be multiplied to get the correct unit (standard 1, i.e. no change)\n",
    "def get_relevant_data(temp, log_df, this_day, name, unit, old_name, calibration=1):\n",
    "    raw_temp = temp.copy()\n",
    "    \n",
    "    # Identify the previous log date compared to the current date of the data\n",
    "    # This log date contains the relevant information to be read and applied\n",
    "    log_entry_date = find_log_entry(log_df['timestamp'].unique(), this_day)\n",
    "    \n",
    "    # Create empty PAR df that will be filled with PAR data\n",
    "    # List all columns in the PAR log without timestamp to have a list of chamber names\n",
    "    cols = log_df.columns.values.tolist()\n",
    "    cols.remove('timestamp')\n",
    "    # Copy all dates in the raw data that will be copied over\n",
    "    out_df = raw_temp[['timestamp']].copy()\n",
    "    # Create the empty df and fill it\n",
    "    for col in cols:\n",
    "        # Create column name\n",
    "        col_name = name + '.' + col + '.' + unit\n",
    "        # Initialise with NAs\n",
    "        out_df[col_name] = np.nan\n",
    "        # Identify MUX port and fill data in\n",
    "        port = log_df.loc[(log_df['timestamp'] == log_entry_date), col].values[0]\n",
    "        if(port > 0): # port 0 is basically bad data\n",
    "            out_df[col_name] = raw_temp[old_name + str(int(port)) + ')']\n",
    "        # If the sensor values need to be multiplied by a factor, e.g. to convert mV to umol.m-2.s-1\n",
    "        out_df[col_name] = out_df[col_name] * calibration\n",
    "    \n",
    "    return(out_df)\n",
    "\n",
    "def add_ambient(temp, amb_df):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Drop unnecessary columns\n",
    "    amb_df.drop('dayid', axis=1, inplace=True)\n",
    "    \n",
    "    # Merge PAR data\n",
    "    merged = temp.merge(amb_df, on='timestamp', how='left')\n",
    "    # Interpolate\n",
    "    merged['temp.air.ambient.c'].interpolate(method='linear', limit=4, limit_direction='both', limit_area='inside', inplace=True)\n",
    "    merged['par.ambient.umol_m2_s1'].interpolate(method='linear', limit=4, limit_direction='both', limit_area='inside', inplace=True)\n",
    "    merged['conc.h2o.irga.ambient.mmol_mol'].interpolate(method='linear', limit=4, limit_direction='both', limit_area='inside', inplace=True)\n",
    "    merged['conc.co2.irga.ambient.umol_mol'].interpolate(method='linear', limit=4, limit_direction='both', limit_area='inside', inplace=True)\n",
    "    \n",
    "    return(merged)\n",
    "\n",
    "# Adds sensor data and then interpolates it to fill the gaps\n",
    "# - temp: Data df\n",
    "# - sensor_df: Dataframe with sensor data\n",
    "# - interpolation_limit: how many values are to be interpolated\n",
    "def add_sensor_data(temp, sensor_df, interpolation_limit):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # List sensor columns only\n",
    "    cols = sensor_df.columns.values.tolist()\n",
    "    cols.remove('timestamp')\n",
    "    \n",
    "    # Merge sensor data\n",
    "    merged = temp.merge(sensor_df, on='timestamp', how='left')\n",
    "    # Interpolate sensor data\n",
    "    for col in cols:\n",
    "        merged[col].interpolate(method='linear', limit=interpolation_limit, limit_direction='both', limit_area='inside', inplace=True)\n",
    "    \n",
    "    return(merged)\n",
    "\n",
    "def fix_PAR(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # List PAR columns only\n",
    "    cols = [col for col in temp.columns if 'par.' in col]\n",
    "    \n",
    "    # Night-time corrections, by cycling through chambers: Set to 0 during the night\n",
    "    temp.loc[temp['par.ambient.umol_m2_s1'] < (par_threshold), 'par.ambient.umol_m2_s1'] = 0\n",
    "    temp.loc[(temp['timestamp'].dt.hour < 5) | (temp['timestamp'].dt.hour > 21) , 'par.ambient.umol_m2_s1'] = 0 # Define the night by the maximum of June 21\n",
    "    for col in cols:\n",
    "        temp.loc[(temp[col] < par_threshold), col] = 0 # All PAR values should be 0 at night. Small values are errors\n",
    "        temp.loc[(temp['par.ambient.umol_m2_s1'] == 0) & (temp[col] > par_threshold), col] = 0 # Remove large Apogee night-time values\n",
    "        \n",
    "    return(temp)\n",
    "\n",
    "def make_counter_and_chamber(temp, this_day):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Prepare the dataframe\n",
    "    temp['type'] = ''\n",
    "    temp['plot'] = ''\n",
    "    temp['chamber'] = '' # Better to make it empty for the \"create_block_id\" function to be faster\n",
    "    temp['status'] = ''  # Better to make it empty for the \"create_block_id\" function to be faster\n",
    "    temp['counter'] = np.nan\n",
    "    temp['chamber_block'] = np.nan\n",
    "    \n",
    "    # Extract the correct timing from the log file. If it's before the first available date, drop all data\n",
    "    if(lt['timestamp'].unique()[0] >= pd.to_datetime(str(this_day), format='%Y%m%d', utc=True)): # Data from before the log file\n",
    "        temp.drop(temp.index, inplace=True)\n",
    "        return(temp)\n",
    "    # Extract the correct timing from the log file\n",
    "    relevant_date = lt.loc[lt['timestamp'] <= this_day,'timestamp'].unique()[-1]\n",
    "    timing_dt = lt.loc[lt['timestamp'] == relevant_date]\n",
    "    \n",
    "    # In case the cycle is longer than 1h, adjust for it\n",
    "    # For 1h long cycles, the times will be divided by 1 simply\n",
    "    cycle_length = np.ceil(timing_dt['end_min'].max()/60)\n",
    "    # Mark the chambers, and add chamber NA for data that needs to be removed\n",
    "    temp['chamber'] = '' # remove data if it remains empty\n",
    "    # Save the hour in the cycle\n",
    "    temp['cycle_hour'] = temp['timestamp'].dt.strftime('%H').astype(float)\n",
    "    temp['cycle_hour'] = temp['cycle_hour'].mod(cycle_length)\n",
    "    # Create the minute relative to the hour in the cycle\n",
    "    temp['chambertiming'] = temp['timestamp'].dt.strftime('%M%S').astype(float)\n",
    "    temp['chambertiming'] = temp['chambertiming'] + temp['cycle_hour']*60*100\n",
    "    # Walk through the log and save the data\n",
    "    for index, row in timing_dt.iterrows():\n",
    "        current_start = row['start_min']*100+row['start_sec']\n",
    "        current_end   = row['end_min']*100+row['end_sec']\n",
    "        # Name of current chamber\n",
    "        if(row['status'] == 'amb'):\n",
    "            current_ch = 'amb'\n",
    "        else:\n",
    "            current_ch = 'ch' + row['type'][0] + '_' + row['plot'] + str(int(row['num']))\n",
    "        # Save the identifiers into the df\n",
    "        temp.loc[(temp['chambertiming'] >= current_start) & \\\n",
    "                 (temp['chambertiming'] <= current_end), 'type'] = row['type']\n",
    "        temp.loc[(temp['chambertiming'] >= current_start) & \\\n",
    "                 (temp['chambertiming'] <= current_end), 'plot'] = row['plot']\n",
    "        temp.loc[(temp['chambertiming'] >= current_start) & \\\n",
    "                 (temp['chambertiming'] <= current_end), 'status'] = row['status']\n",
    "        temp.loc[(temp['chambertiming'] >= current_start) & \\\n",
    "                 (temp['chambertiming'] <= current_end), 'chamber'] = current_ch\n",
    "    temp.drop(['chambertiming', 'cycle_hour'], axis=1, inplace=True)\n",
    "    \n",
    "    # Add counter\n",
    "    temp = create_block_id(temp)\n",
    "    temp['counter'] = temp.groupby(['chamber_block']).cumcount()\n",
    "    # Revert the counter to end at 0\n",
    "    for name, group in temp.groupby(['chamber_block']):\n",
    "        temp.loc[temp['chamber_block'] == name, 'counter'] = group['counter'] - group['timestamp'].count() + 1\n",
    "    \n",
    "    # Remove chamber number when not relevant\n",
    "    temp.loc[~temp['status'].isin(['oc','cc','amb','']), 'chamber'] = ''\n",
    "    # Make ambient plot into ambient\n",
    "    temp.loc[temp['status'] == 'amb', 'plot'] = 'amb'\n",
    "    \n",
    "    # Chamber with NAs are a dummy that needs to be removed\n",
    "    temp = temp.loc[temp['chamber'] != '']\n",
    "    \n",
    "    return(temp)\n",
    "\n",
    "# Define day & night by PAR\n",
    "def define_day_night(temp):\n",
    "    temp = temp.copy()\n",
    "    # Define nights as 0, days as 1\n",
    "    temp['day_night'] = 0\n",
    "    temp.loc[temp['par.ambient.umol_m2_s1'] > day_par_min, 'day_night'] = 1\n",
    "    \n",
    "    return(temp)\n",
    "\n",
    "def add_current_variable(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Add a 'current' column for par temperature and flowrate, which has the reading for the current chamber\n",
    "    temp['par.current.chamber.umol_m2_s1'] = np.nan\n",
    "    temp['temp.leaf.current.chamber.c'] = np.nan\n",
    "    temp['temp.soil.current.chamber.c'] = np.nan\n",
    "    temp['temp.air.current.chamber.c'] = np.nan\n",
    "    temp['pump.flow.current.chamber.slpm'] = np.nan\n",
    "    # Loop through chambers and fill in the gaps\n",
    "    chamber_list = temp['chamber'].dropna().unique().tolist()\n",
    "    if('amb' in chamber_list): # Ambient 'chamber' should not be treated\n",
    "        chamber_list.remove('amb')\n",
    "    for i in chamber_list:\n",
    "        if(('par.' + i + '.umol_m2_s1') in temp.columns.values):\n",
    "            temp.loc[temp['chamber'] == i, 'par.current.chamber.umol_m2_s1'] = temp.loc[temp['chamber'] == i, 'par.' + i + '.umol_m2_s1']\n",
    "        if(('temp.leaf.' + i + '.c') in temp.columns.values):\n",
    "            temp.loc[temp['chamber'] == i, 'temp.leaf.current.chamber.c'] = temp.loc[temp['chamber'] == i, 'temp.leaf.' + i + '.c']\n",
    "        if(('temp.soil.' + i + '.c') in temp.columns.values):\n",
    "            temp.loc[temp['chamber'] == i, 'temp.soil.current.chamber.c'] = temp.loc[temp['chamber'] == i, 'temp.soil.' + i + '.c']\n",
    "        if(('temp.air.' + i + '.c') in temp.columns.values):\n",
    "            temp.loc[temp['chamber'] == i, 'temp.air.current.chamber.c'] = temp.loc[temp['chamber'] == i, 'temp.air.' + i + '.c']\n",
    "        if(('pump.flow.' + i + '.slpm') in temp.columns.values):\n",
    "            temp.loc[temp['chamber'] == i, 'pump.flow.current.chamber.slpm'] = temp.loc[temp['chamber'] == i, 'pump.flow.' + i + '.slpm']\n",
    "    # When the chamber is NA, use ambient data\n",
    "    temp.loc[temp['chamber'].isna(), 'par.current.chamber.umol_m2_s1'] = temp.loc[temp['chamber'].isna(), 'par.ambient.umol_m2_s1']\n",
    "    temp.loc[temp['chamber'].isna(), 'temp.air.current.chamber.c'] = temp.loc[temp['chamber'].isna(), 'temp.air.ambient.c']\n",
    "    return(temp)\n",
    "\n",
    "def add_surface_area(temp, this_day, log_df):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Identify the previous log date compared to the current date of the data\n",
    "    # This log date contains the relevant information to be read and applied\n",
    "    log_entry_date = find_log_entry(log_df['timestamp'].unique(), this_day)\n",
    "    \n",
    "    # List sensor columns only\n",
    "    chamber_list = log_df.columns.values.tolist()\n",
    "    chamber_list.remove('timestamp')\n",
    "    \n",
    "    # Add surface area\n",
    "    temp['surface.area.m2'] = np.nan\n",
    "    for ch in chamber_list:\n",
    "        temp.loc[(temp['chamber'] == ch), 'surface.area.m2'] = log_df.loc[(log_df['timestamp'] == log_entry_date), ch].values[0]\n",
    "    \n",
    "    return(temp)\n",
    "\n",
    "# Remove data from the bad chambers, based on the logbook\n",
    "def add_log_flow(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Identify day\n",
    "    day = temp['timestamp'].dt.normalize().unique()[0]\n",
    "    \n",
    "    # Keep only relevant days, make a list of bad chambers\n",
    "    column_list = [cols for cols in cs.columns.tolist() if 'ch' in cols] # List only columns with 'ch' (i.e. chambers)\n",
    "    cs_df = cs.loc[cs['timestamp'] == day, column_list]\n",
    "    \n",
    "    # Replace the pump flow if a constant was given in the logbook, or remove values if the pump wasn't working (according to the logbook)\n",
    "    for chamber in cs_df.columns.tolist():\n",
    "        if(cs_df[chamber].values[0] > 0):\n",
    "            temp.loc[temp['chamber'] == chamber, 'pump.flow.current.chamber.slpm'] = cs_df[chamber].values[0]\n",
    "        if(cs_df[chamber].values[0] == 0): # pump not working\n",
    "            temp.loc[temp['chamber'] == chamber, 'pump.flow.current.chamber.slpm'] = np.nan\n",
    "    \n",
    "    # Add a column of unfiltered \"current\" flow data\n",
    "    temp['pump.flow.unfiltered.current.chamber.slpm'] = temp['pump.flow.current.chamber.slpm']\n",
    "    # Consider the pump as not working at a flow below the limit\n",
    "    temp.loc[(temp['type'] == 'branch') & (temp['pump.flow.current.chamber.slpm'] <= min_abs_flow_branch), 'pump.flow.current.chamber.slpm'] = np.nan\n",
    "    temp.loc[(temp['type'] == 'soil')   & (temp['pump.flow.current.chamber.slpm'] <= min_abs_flow_soil), 'pump.flow.current.chamber.slpm'] = np.nan\n",
    "    \n",
    "    return(temp)\n",
    "\n",
    "# Remove data from the bad chambers, based on the logbook\n",
    "def remove_bad_chambers(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Identify day\n",
    "    day = temp['timestamp'].dt.normalize().unique()[0]\n",
    "    \n",
    "    # Keep only relevant days, make a list of bad chambers\n",
    "    column_list = [cols for cols in cs.columns.tolist() if 'ch' in cols] # List only columns with 'ch' (i.e. chambers)\n",
    "    cs_df = cs.loc[cs['timestamp'] == day, column_list]\n",
    "    bad_chambers = cs_df.columns[(cs_df == 0).all()].tolist()\n",
    "    \n",
    "    # Show information\n",
    "    print('    Non-working chambers (from logbook): ', bad_chambers)\n",
    "    \n",
    "    # Remove data when entire chamber was bad\n",
    "    # COS (Carbonyl sulfide)\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.cos.1.pmol_mol'] = np.nan\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.cos.2.pmol_mol'] = np.nan\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.cos.mean.pmol_mol'] = np.nan\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.cos.mean_sd.pmol_mol'] = np.nan\n",
    "    # CO2 (Carbon dioxide)\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.co2.1.umol_mol'] = np.nan\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.co2.2.umol_mol'] = np.nan\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.co2.mean.umol_mol'] = np.nan\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.co2.mean_sd.umol_mol'] = np.nan\n",
    "    # H2O (Water)\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.h2o.mmol_mol'] = np.nan\n",
    "    # CO (Carbon monoxide)\n",
    "    temp.loc[temp['chamber'].isin(bad_chambers),'conc.co.nmol_mol'] = np.nan\n",
    "    \n",
    "    return(temp)\n",
    "\n",
    "# Writes output files in full temporal resolution\n",
    "def write_output_file(out_df, date_idx, out_dir):\n",
    "    # Check if output folders exist. If not, create\n",
    "    month_dir = str(date_idx)[0:4] + '-' + str(date_idx)[4:6]\n",
    "    if(not os.path.exists(out_dir + month_dir)): # make directory if it doesn't exist\n",
    "        os.makedirs(out_dir + month_dir)\n",
    "    # Create file name\n",
    "    out_fn = out_dir + month_dir + '/' + output_fn + '_' + str(date_idx) + '.csv'\n",
    "    print(str(len(out_df.dayid)), out_fn) # Shows final file size)\n",
    "    # organise data for output\n",
    "    temp_df = out_df\n",
    "    # Before saving, remove the index\n",
    "    temp_df.drop('dayid', axis=1, inplace=True)\n",
    "    #temp_df.drop('chamber_block', axis=1, inplace=True)\n",
    "    # Move timestamp column to the front\n",
    "    col = temp_df.pop('counter')\n",
    "    temp_df.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = temp_df.pop('chamber_block')\n",
    "    temp_df.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = temp_df.pop('status')\n",
    "    temp_df.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = temp_df.pop('chamber')\n",
    "    temp_df.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = temp_df.pop('type')\n",
    "    temp_df.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = temp_df.pop('plot')\n",
    "    temp_df.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = temp_df.pop('timestamp')\n",
    "    temp_df.insert(0, col.name, col, allow_duplicates=True)\n",
    "    # Remove timezone information\n",
    "    temp_df['timestamp'] = temp_df['timestamp'].dt.tz_localize(None)\n",
    "    # Write data\n",
    "    temp_df.to_csv(out_fn, index=False, encoding='utf-8', date_format='%Y-%m-%d %H:%M:%S') # Save file\n",
    "    pass\n",
    "\n",
    "# Removes all data that is not in the relevant timeframe\n",
    "def remove_obsolete_data(temp, seconds):\n",
    "    temp = temp.copy()\n",
    "    # Remove data, if 'keep_all_raw' is 'True', keep all data\n",
    "    if(keep_all_raw):\n",
    "        pass\n",
    "    else:\n",
    "        temp = temp.loc[temp['counter'] > (-seconds)]\n",
    "    return(temp)\n",
    "\n",
    "# Calculate means and standard deviations of all data columns\n",
    "def make_means(temp, seconds):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Prepare timestamps: Convert to numbers so that averaging becomes possible\n",
    "    temp['timestamp'] = temp['timestamp'].astype('int64')\n",
    "    \n",
    "    # Remove obsolete data if it wasn't removed previously\n",
    "    if(~keep_all_raw):\n",
    "        temp = temp.loc[temp['counter'] > (-seconds)]\n",
    "    \n",
    "    # Keep relevant identification columns\n",
    "    id_cols = temp[['chamber_block','type','plot','chamber','status']].drop_duplicates()\n",
    "    \n",
    "    # Make standard deviations\n",
    "    sd_df = temp.groupby(['chamber_block']).std(numeric_only=True)\n",
    "    sd_df.reset_index(inplace=True)\n",
    "    sd_df.drop('counter', axis=1, inplace=True)\n",
    "    sd_df.drop('timestamp', axis=1, inplace=True)\n",
    "    #sd_df.drop('chamber', axis=1, inplace=True)\n",
    "    #sd_df.drop('plot', axis=1, inplace=True)\n",
    "    sd_df.drop('surface.area.m2', axis=1, inplace=True)\n",
    "    sd_df = sd_df.add_suffix('.sd')\n",
    "    sd_df.rename(columns={'chamber_block.sd':'chamber_block'}, inplace=True)\n",
    "    # Make means\n",
    "    avg_df = temp.groupby(['chamber_block']).mean(numeric_only=True)\n",
    "    avg_df.drop('counter', axis=1, inplace=True)\n",
    "    avg_df.reset_index(inplace=True)\n",
    "    # Convert averaged 'timestamps stored as numbers' back to (rounded) timestamps\n",
    "    avg_df['timestamp'] = avg_df['timestamp'].astype('datetime64[ns]').dt.round('10s')\n",
    "    \n",
    "    # Merge means with standard deviations\n",
    "    summarised = avg_df.merge(sd_df, on=['chamber_block'], how='left')\n",
    "    # Remove weird indices\n",
    "    summarised.reset_index(inplace=True)\n",
    "    summarised.drop('index', axis=1, inplace=True)\n",
    "    \n",
    "    # Add identification columns (status and plot). We need to do this because text is not preserved when making averages\n",
    "    summarised = summarised.merge(id_cols, on='chamber_block', how='left')\n",
    "    # Move important columns to the front\n",
    "    col = summarised.pop('status')\n",
    "    summarised.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = summarised.pop('chamber')\n",
    "    summarised.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = summarised.pop('plot')\n",
    "    summarised.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = summarised.pop('type')\n",
    "    summarised.insert(0, col.name, col, allow_duplicates=True)\n",
    "    col = summarised.pop('timestamp')\n",
    "    summarised.insert(0, col.name, col, allow_duplicates=True)\n",
    "    \n",
    "    # Fix: Calculate mean of means, and their standard deviations, of the gases with multiple columns (cos, co2)\n",
    "    # (1) Remove previous standard deviations\n",
    "    summarised.drop('conc.cos.mean_sd.pmol_mol', axis=1, inplace=True)\n",
    "    summarised.drop('conc.co2.mean_sd.umol_mol', axis=1, inplace=True)\n",
    "    summarised.drop('conc.cos.mean_sd.pmol_mol.sd', axis=1, inplace=True)\n",
    "    summarised.drop('conc.co2.mean_sd.umol_mol.sd', axis=1, inplace=True)\n",
    "    # (2) Calculate mean of means\n",
    "    summarised['conc.cos.mean.pmol_mol'] = summarised[['conc.cos.1.pmol_mol','conc.cos.2.pmol_mol']].mean(axis=1)\n",
    "    summarised['conc.co2.mean.umol_mol'] = summarised[['conc.co2.1.umol_mol','conc.co2.2.umol_mol']].mean(axis=1)\n",
    "    # (3) Calculate standard deviation of 'means of means'\n",
    "    summarised['conc.cos.mean.pmol_mol.sd'] = std_means_of_means([summarised['conc.cos.1.pmol_mol.sd'], summarised['conc.cos.2.pmol_mol.sd']])\n",
    "    summarised['conc.co2.mean.umol_mol.sd'] = std_means_of_means([summarised['conc.co2.1.umol_mol.sd'], summarised['conc.co2.2.umol_mol.sd']])\n",
    "    \n",
    "    return(summarised)\n",
    "\n",
    "# Interpolates ambient and open chamber readings into 1min readings. Then these are taken for the ambient reading\n",
    "def find_corresponding_ambient_interpolated(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Create empty dataframe with a 1min frequency\n",
    "    idx = pd.date_range(start=temp['timestamp'].tolist()[0],\n",
    "                        end=temp['timestamp'].tolist()[-1],\n",
    "                        freq='s')\n",
    "    time_df = pd.DataFrame(idx, index=None, columns=['timestamp'])\n",
    "    \n",
    "    # Create a list of relevant columns\n",
    "    column_list = [cols for cols in temp.columns.tolist() if 'conc' in cols] # List only columns with concentrations\n",
    "    column_list = [cols for cols in column_list if not cols.endswith('.sd')] # Ignore standard deviation columns\n",
    "    column_list = [cols for cols in column_list if not ('.irga.' in cols) and not ('.calibrated.' in cols)] # Ignore calibrated and IRGA columns\n",
    "    column_list = column_list + ['temp.leaf.current.chamber.c', 'temp.soil.current.chamber.c', 'temp.air.current.chamber.c']\n",
    "    \n",
    "    # Extract ambient & interpolate\n",
    "    ambient_df = time_df.merge(temp.loc[(temp['status'] == 'amb'),['timestamp'] + column_list], on='timestamp', how='left')\n",
    "    col = ambient_df.pop('timestamp') # Interpolation doesn't work on timestamps, temporarily take it out\n",
    "    ambient_df = ambient_df.interpolate(method='linear', limit=3600, limit_direction='both')\n",
    "    ambient_df.insert(0, col.name, col, allow_duplicates=True) # And put the timestamp back\n",
    "    # Rename columns\n",
    "    ambient_df = ambient_df.add_suffix('.ambient')\n",
    "    ambient_df.rename(columns={'timestamp.ambient':'timestamp'}, inplace=True)\n",
    "    \n",
    "    # Extract open chamber & interpolate, if it exists. Create an empty df if not\n",
    "    oc_data_list = []\n",
    "    chamber_list = temp['chamber'].unique().tolist()\n",
    "    # Interpolate chamber by chamber\n",
    "    for ch in chamber_list:\n",
    "        if('oc' in temp['status'].unique()):\n",
    "            oc_temp = time_df.merge(temp.loc[(temp['chamber'] == ch) & (temp['status'] == 'oc'),['timestamp'] + column_list], on='timestamp', how='left')\n",
    "            col = oc_temp.pop('timestamp') # Interpolation doesn't work on timestamps, temporarily take it out\n",
    "            oc_temp = oc_temp.interpolate(method='linear', limit=3600, limit_direction='both')\n",
    "            oc_temp.insert(0, col.name, col, allow_duplicates=True) # And put the timestamp back\n",
    "        else: # Add multiple columns to the df, but leave them empty\n",
    "            oc_temp = time_df.reindex(columns=['timestamp']+column_list)\n",
    "        # Add chamber label\n",
    "        oc_temp['chamber'] = ch\n",
    "        oc_data_list.append(oc_temp)\n",
    "    oc_df = pd.concat(oc_data_list) # Now merge all\n",
    "    # Rename columns\n",
    "    oc_df = oc_df.add_suffix('.oc')\n",
    "    oc_df.rename(columns={'timestamp.oc':'timestamp'}, inplace=True)\n",
    "    oc_df.rename(columns={'chamber.oc':'chamber'}, inplace=True)\n",
    "    \n",
    "    # Merge with data. Note: Removes the ambient from the temp df before merging, so that no ambient-to-ambient flux is calculated\n",
    "    merged = temp.loc[temp['chamber'] != 'amb'].merge(ambient_df, on='timestamp', how='left')\n",
    "    merged = merged.merge(oc_df, on=['timestamp','chamber'], how='left')\n",
    "    \n",
    "    return(merged)\n",
    "\n",
    "# Find the next ambient or open chamber after the current measurement, mark it, then append it\n",
    "def find_corresponding_ambient_next(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Create a list of relevant columns\n",
    "    column_list = [cols for cols in temp.columns.tolist() if 'conc' in cols] # List only columns with concentrations\n",
    "    column_list = [cols for cols in column_list if not cols.endswith('.sd')] # Ignore standard deviation columns\n",
    "    column_list = [cols for cols in column_list if not ('.irga.' in cols) and not ('.calibrated.' in cols)] # Ignore calibrated and IRGA columns\n",
    "    column_list = column_list + ['temp.leaf.current.chamber.c', 'temp.soil.current.chamber.c', 'temp.air.current.chamber.c']\n",
    "    \n",
    "    # Extract ambient & prepare\n",
    "    amb_df = temp.loc[(temp['status'] == 'amb'),['chamber_block'] + column_list].copy()\n",
    "    amb_df = amb_df.add_suffix('.ambient') # Rename columns\n",
    "    amb_df.rename(columns={'chamber_block.ambient':'chamber_block'}, inplace=True)\n",
    "    # Extract open chamber & prepare\n",
    "    oc_df = temp.loc[(temp['status'] == 'oc'),['chamber_block'] + column_list].copy()\n",
    "    oc_df = oc_df.add_suffix('.oc') # Rename columns\n",
    "    oc_df.rename(columns={'chamber_block.oc':'chamber_block'}, inplace=True)\n",
    "    # Ambient column to hold the chamber block number\n",
    "    amb_df['amb_block'] = np.nan\n",
    "    oc_df['amb_block'] = np.nan\n",
    "    \n",
    "    # Search for the corresponding next ambient or open chamber, do it chamber by chamber\n",
    "    chamber_list = temp['chamber'].unique().tolist()\n",
    "    amb_list = [ ]\n",
    "    for ch in chamber_list:\n",
    "        temp_amb = amb_df.copy() # Ambient for each chamber is treated separately, and then concatenated (below)\n",
    "        # Cycle through closed chambers\n",
    "        cc_block_list = temp.loc[(temp['status'] == 'cc') & (temp['chamber'] == ch), 'chamber_block'].unique().tolist()\n",
    "        # Make sure the open chamber is the same chamber\n",
    "        oc_block_list = temp.loc[(temp['status'] == 'oc') & (temp['chamber'] == ch), 'chamber_block'].unique().tolist()\n",
    "        # For ambient, we don't have to take the same chamber\n",
    "        amb_block_list = temp.loc[(temp['status'] == 'amb'), 'chamber_block'].unique().tolist()\n",
    "        \n",
    "        for block in cc_block_list:\n",
    "            # Find next block with amb or oc. Avoid ValueError if the list is empty\n",
    "            if(len(oc_block_list) > 0):\n",
    "                if(block < np.max(oc_block_list)):\n",
    "                    oc_block =  np.min([x for x in oc_block_list if x > block])\n",
    "            else:\n",
    "                oc_block =  np.nan\n",
    "            if(len(amb_block_list) > 0):\n",
    "                if(block < np.max(amb_block_list)):\n",
    "                    amb_block =  np.min([x for x in amb_block_list if x > block])\n",
    "            else:\n",
    "                amb_block =  np.nan\n",
    "            # In the ambient and OC dfs, replace the chamber_block with the current block\n",
    "            temp_amb.loc[temp_amb['chamber_block'] == amb_block, 'amb_block'] = block\n",
    "            oc_df.loc[oc_df['chamber_block'] == oc_block, 'amb_block'] = block\n",
    "        # If the current chamber is not ambient, add this to the list of ambient blocks that were found\n",
    "        if(ch != 'amb'):\n",
    "            amb_list.append(temp_amb.loc[~temp_amb['amb_block'].isna()])\n",
    "    amb_df = pd.concat(amb_list)\n",
    "    \n",
    "    # Sometimes, the is another OC or ambient measured before a chamber (which will be ignored, or NaN). Remove it\n",
    "    amb_df['chamber_block'] = amb_df['amb_block']\n",
    "    amb_df = amb_df.loc[~amb_df['amb_block'].isna()]\n",
    "    amb_df.drop('amb_block', axis=1, inplace=True)\n",
    "    oc_df['chamber_block'] = oc_df['amb_block']\n",
    "    oc_df = oc_df.loc[~oc_df['amb_block'].isna()]\n",
    "    oc_df.drop('amb_block', axis=1, inplace=True)\n",
    "    # Make sure the dtype is the same as the chamber block\n",
    "    amb_df['chamber_block'] = amb_df['chamber_block'].astype(int)\n",
    "    oc_df['chamber_block'] = oc_df['chamber_block'].astype(int)\n",
    "    \n",
    "    # Merge with data. Note: Removes the ambient from the temp df before merging, so that no ambient-to-ambient flux is calculated\n",
    "    merged = temp.loc[temp['status'] != 'amb'].merge(amb_df, on='chamber_block', how='left')\n",
    "    merged = merged.merge(oc_df, on='chamber_block', how='left')\n",
    "    \n",
    "    return(merged)\n",
    "\n",
    "# Calculate all fluxes for all gases\n",
    "def calculate_fluxes(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # The TSI flow meters use SLPM with the following standard values for air flow:\n",
    "    air_temperature_tsi_c = 21.11 # [°C]\n",
    "    air_pressure_tsi_kpa  = 101.3 # [kPa]\n",
    "    \n",
    "    # Fluxes with total surface area, between chamber and ambient measurements\n",
    "    temp['flux.h2o.ch_amb.mmol_m2_s'] = 10**3 * calculate_h2o_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    temp['flux.co2.ch_amb.umol_m2_s1'] = 10**6 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co2.mean.umol_mol.ambient'] * 10**(-6),\n",
    "                                                            temp['conc.co2.mean.umol_mol'] * 10**(-6),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    temp['flux.cos.ch_amb.pmol_m2_s'] = 10**12 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.cos.mean.pmol_mol.ambient'] * 10**(-12),\n",
    "                                                            temp['conc.cos.mean.pmol_mol'] * 10**(-12),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    temp['flux.co.ch_amb.nmol_m2_s'] = 10**9 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co.nmol_mol.ambient'] * 10**(-9),\n",
    "                                                            temp['conc.co.nmol_mol'] * 10**(-9),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    \n",
    "    # Fluxes with projected leaf area, between chamber and ambient measurements\n",
    "    temp['flux.h2o.ch_amb.proj_la.mmol_m2_s'] = 10**3 * calculate_h2o_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    temp['flux.co2.ch_amb.proj_la.umol_m2_s1'] = 10**6 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co2.mean.umol_mol.ambient'] * 10**(-6),\n",
    "                                                            temp['conc.co2.mean.umol_mol'] * 10**(-6),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    temp['flux.cos.ch_amb.proj_la.pmol_m2_s'] = 10**12 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.cos.mean.pmol_mol.ambient'] * 10**(-12),\n",
    "                                                            temp['conc.cos.mean.pmol_mol'] * 10**(-12),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    temp['flux.co.ch_amb.proj_la.nmol_m2_s'] = 10**9 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co.nmol_mol.ambient'] * 10**(-9),\n",
    "                                                            temp['conc.co.nmol_mol'] * 10**(-9),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    # For soil chambers, there is no projected leaf area\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.h2o.ch_amb.proj_la.mmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.co2.ch_amb.proj_la.umol_m2_s1'] = np.nan\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.cos.ch_amb.proj_la.pmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.co.ch_amb.proj_la.nmol_m2_s'] = np.nan\n",
    "    \n",
    "    # Fluxes with total surface area, between chamber and open chamber measurements\n",
    "    temp['flux.h2o.ch_oc.mmol_m2_s'] = 10**3 * calculate_h2o_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    temp['flux.co2.ch_oc.umol_m2_s1'] = 10**6 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co2.mean.umol_mol.oc'] * 10**(-6),\n",
    "                                                            temp['conc.co2.mean.umol_mol'] * 10**(-6),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    temp['flux.cos.ch_oc.pmol_m2_s'] = 10**12 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.cos.mean.pmol_mol.oc'] * 10**(-12),\n",
    "                                                            temp['conc.cos.mean.pmol_mol'] * 10**(-12),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    temp['flux.co.ch_oc.nmol_m2_s'] = 10**9 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co.nmol_mol.oc'] * 10**(-9),\n",
    "                                                            temp['conc.co.nmol_mol'] * 10**(-9),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2'])\n",
    "    \n",
    "    # Fluxes with projected leaf area, between chamber and open chamber measurements\n",
    "    temp['flux.h2o.ch_oc.proj_la.mmol_m2_s'] = 10**3 * calculate_h2o_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    temp['flux.co2.ch_oc.proj_la.umol_m2_s1'] = 10**6 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co2.mean.umol_mol.oc'] * 10**(-6),\n",
    "                                                            temp['conc.co2.mean.umol_mol'] * 10**(-6),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    temp['flux.cos.ch_oc.proj_la.pmol_m2_s'] = 10**12 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.cos.mean.pmol_mol.oc'] * 10**(-12),\n",
    "                                                            temp['conc.cos.mean.pmol_mol'] * 10**(-12),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    temp['flux.co.ch_oc.proj_la.nmol_m2_s'] = 10**9 * calculate_gas_flux(air_temperature_tsi_c,\n",
    "                                                            air_pressure_tsi_kpa*1000,\n",
    "                                                            temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                            temp['conc.h2o.mmol_mol'],\n",
    "                                                            temp['conc.co.nmol_mol.oc'] * 10**(-9),\n",
    "                                                            temp['conc.co.nmol_mol'] * 10**(-9),\n",
    "                                                            temp['pump.flow.current.chamber.slpm'],\n",
    "                                                            temp['surface.area.m2']/proj_la_factor)\n",
    "    # For soil chambers, there is no projected leaf area\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.h2o.ch_oc.proj_la.mmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.co2.ch_oc.proj_la.umol_m2_s1'] = np.nan\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.cos.ch_oc.proj_la.pmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['type'] == 'soil', 'flux.co.ch_oc.proj_la.nmol_m2_s'] = np.nan\n",
    "    \n",
    "    # There is no flux between open chamber and open chamber\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.h2o.ch_oc.mmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.co2.ch_oc.umol_m2_s1'] = np.nan\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.cos.ch_oc.pmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.co.ch_oc.nmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.h2o.ch_oc.proj_la.mmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.co2.ch_oc.proj_la.umol_m2_s1'] = np.nan\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.cos.ch_oc.proj_la.pmol_m2_s'] = np.nan\n",
    "    temp.loc[temp['status'] == 'oc', 'flux.co.ch_oc.proj_la.nmol_m2_s'] = np.nan\n",
    "    \n",
    "    return(temp)\n",
    "\n",
    "# Calculate all fluxes for all gases\n",
    "def calculate_conductances(temp):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Calculate VPD\n",
    "    temp['VPD.Pa.ambient'] = calculate_VPD(temp['temp.air.current.chamber.c'], temp['conc.h2o.mmol_mol.ambient'], temp['P.irga.kPa']*1000)\n",
    "    temp['VPD.Pa.oc'] = calculate_VPD(temp['temp.air.current.chamber.c.oc'], temp['conc.h2o.mmol_mol.oc'], temp['P.irga.kPa']*1000)\n",
    "    temp['VPD.Pa'] = calculate_VPD(temp['temp.air.current.chamber.c'], temp['conc.h2o.mmol_mol'], temp['P.irga.kPa']*1000)\n",
    "    \n",
    "    # Conductances for COS with total leaf area, between chamber and ambient measurements\n",
    "    temp['cond.stomatal.ch_amb.cos'] = calculate_cos_stomatal_conductance_ball_berry(temp['temp.air.current.chamber.c'],\n",
    "                                                                    temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                                    temp['P.irga.kPa']*1000,\n",
    "                                                                    temp['flux.h2o.ch_amb.mmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_amb.umol_m2_s1'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.ambient'],\n",
    "                                                                    temp['par.ambient.umol_m2_s1'])\n",
    "    temp['cond.total.ch_amb.cos'] = calculate_cos_total_conductance(temp['flux.cos.ch_amb.pmol_m2_s'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.ambient'])\n",
    "    temp['cond.internal.ch_amb.cos'] = calculate_cos_internal_conductance(temp['cond.total.ch_amb.cos'],\n",
    "                                                                    temp['cond.stomatal.ch_amb.cos'])\n",
    "    temp['rel_uptake.ch_amb'] = relative_uptake(temp['flux.cos.ch_amb.pmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_amb.umol_m2_s1'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.ambient'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.ambient'])\n",
    "    \n",
    "    # Conductances for COS with projected leaf area, between chamber and ambient measurements\n",
    "    temp['cond.stomatal.ch_amb.proj_la.cos'] = calculate_cos_stomatal_conductance_ball_berry(temp['temp.air.current.chamber.c'],\n",
    "                                                                    temp['conc.h2o.mmol_mol.ambient'],\n",
    "                                                                    temp['P.irga.kPa']*1000,\n",
    "                                                                    temp['flux.h2o.ch_amb.proj_la.mmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_amb.proj_la.umol_m2_s1'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.ambient'],\n",
    "                                                                    temp['par.ambient.umol_m2_s1'])\n",
    "    temp['cond.total.ch_amb.proj_la.cos'] = calculate_cos_total_conductance(temp['flux.cos.ch_amb.proj_la.pmol_m2_s'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.ambient'])\n",
    "    temp['cond.internal.ch_amb.proj_la.cos'] = calculate_cos_internal_conductance(temp['cond.total.ch_amb.proj_la.cos'],\n",
    "                                                                    temp['cond.stomatal.ch_amb.proj_la.cos'])\n",
    "    temp['rel_uptake.ch_amb.proj_la'] = relative_uptake(temp['flux.cos.ch_amb.proj_la.pmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_amb.proj_la.umol_m2_s1'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.ambient'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.ambient'])\n",
    "    \n",
    "    # Conductances for COS with total leaf area, between chamber and open chamber measurements\n",
    "    temp['cond.stomatal.ch_oc.cos'] = calculate_cos_stomatal_conductance_ball_berry(temp['temp.air.current.chamber.c'],\n",
    "                                                                    temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                                    temp['P.irga.kPa']*1000,\n",
    "                                                                    temp['flux.h2o.ch_oc.mmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_oc.umol_m2_s1'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.oc'],\n",
    "                                                                    temp['par.ambient.umol_m2_s1'])\n",
    "    temp['cond.total.ch_oc.cos'] = calculate_cos_total_conductance(temp['flux.cos.ch_oc.pmol_m2_s'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.oc'])\n",
    "    temp['cond.internal.ch_oc.cos'] = calculate_cos_internal_conductance(temp['cond.total.ch_oc.cos'],\n",
    "                                                                    temp['cond.stomatal.ch_oc.cos'])\n",
    "    temp['rel_uptake.ch_oc'] = relative_uptake(temp['flux.cos.ch_oc.pmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_oc.umol_m2_s1'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.oc'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.oc'])\n",
    "    \n",
    "    # Conductances for COS with projected leaf area, between chamber and open chamber measurements\n",
    "    temp['cond.stomatal.ch_oc.proj_la.cos'] = calculate_cos_stomatal_conductance_ball_berry(temp['temp.air.current.chamber.c'],\n",
    "                                                                    temp['conc.h2o.mmol_mol.oc'],\n",
    "                                                                    temp['P.irga.kPa']*1000,\n",
    "                                                                    temp['flux.h2o.ch_oc.proj_la.mmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_oc.proj_la.umol_m2_s1'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.oc'],\n",
    "                                                                    temp['par.ambient.umol_m2_s1'])\n",
    "    temp['cond.total.ch_oc.proj_la.cos'] = calculate_cos_total_conductance(temp['flux.cos.ch_oc.proj_la.pmol_m2_s'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.oc'])\n",
    "    temp['cond.internal.ch_oc.proj_la.cos'] = calculate_cos_internal_conductance(temp['cond.total.ch_oc.proj_la.cos'],\n",
    "                                                                    temp['cond.stomatal.ch_oc.proj_la.cos'])\n",
    "    temp['rel_uptake.ch_oc.proj_la'] = relative_uptake(temp['flux.cos.ch_oc.proj_la.pmol_m2_s'],\n",
    "                                                                    temp['flux.co2.ch_oc.proj_la.umol_m2_s1'],\n",
    "                                                                    temp['conc.cos.mean.pmol_mol.oc'],\n",
    "                                                                    temp['conc.co2.mean.umol_mol.oc'])\n",
    "    \n",
    "    # Biochemical conductance\n",
    "    temp['cond.biochemical'] = np.nan\n",
    "    temp.loc[temp['plot'] == 'irr', 'cond.biochemical'] = calculate_biochemical_conductance(temp['temp.leaf.current.chamber.c'], LAI_irr)\n",
    "    temp.loc[temp['plot'] == 'ctr', 'cond.biochemical'] = calculate_biochemical_conductance(temp['temp.leaf.current.chamber.c'], LAI_ctr)\n",
    "    \n",
    "    # Mesophyll conductance\n",
    "    temp['cond.mesophyll'] = np.nan\n",
    "    temp.loc[temp['plot'] == 'irr', 'cond.mesophyll'] = calculate_mesophyll_conductance(temp['temp.leaf.current.chamber.c'], LAI_irr)\n",
    "    temp.loc[temp['plot'] == 'ctr', 'cond.mesophyll'] = calculate_mesophyll_conductance(temp['temp.leaf.current.chamber.c'], LAI_ctr)\n",
    "\n",
    "    return(temp)\n",
    "\n",
    "def calibrate_co2_h2o(temp, cal_stats):\n",
    "    temp = temp.copy()\n",
    "    \n",
    "    # Correlation of H2O\n",
    "    print('    Correlation between IRGA ambient and laser')\n",
    "    h2o_mask = (temp['status'] == 'amb') & (temp['conc.h2o.irga.ambient.mmol_mol'].notna()) & (temp['conc.h2o.mmol_mol'].notna())\n",
    "    n_h2o = len(temp.loc[h2o_mask,'conc.h2o.mmol_mol'])\n",
    "    if(n_h2o >= 10):\n",
    "        h2o_corr = stats.linregress(temp.loc[h2o_mask,'conc.h2o.mmol_mol'], temp.loc[h2o_mask,'conc.h2o.irga.ambient.mmol_mol'])\n",
    "        print('        H2O: ', h2o_corr, n_h2o)\n",
    "    # Correlation of CO2\n",
    "    co2_mask = (temp['status'] == 'amb') & (temp['conc.co2.irga.ambient.umol_mol'].notna()) & (temp['conc.co2.mean.umol_mol'].notna())\n",
    "    n_co2 = len(temp.loc[co2_mask,'conc.co2.mean.umol_mol'])\n",
    "    if(n_co2 >= 10):\n",
    "        co2_corr = stats.linregress(temp.loc[co2_mask,'conc.co2.mean.umol_mol'], temp.loc[co2_mask,'conc.co2.irga.ambient.umol_mol'])\n",
    "        print('        CO2: ', co2_corr, n_co2)\n",
    "    \n",
    "    # Flag the data:\n",
    "    # H2O bad  = 1\n",
    "    # CO2 bad  = 2\n",
    "    # Both bad = 3\n",
    "    temp['flag'] = np.nan\n",
    "    if((n_h2o < 10) | (n_co2 < 10)):\n",
    "        temp['flag'] = 0\n",
    "    else:\n",
    "        if(h2o_corr.rvalue < 0.90):\n",
    "            temp['flag'] = 1\n",
    "        if(co2_corr.rvalue < 0.90):\n",
    "            temp['flag'] = 2\n",
    "        if((h2o_corr.rvalue < 0.90) & (co2_corr.rvalue < 0.90)):\n",
    "            temp['flag'] = 3\n",
    "        \n",
    "    # Add information to the calibration statistics df\n",
    "    temp_cal_df = pd.DataFrame(columns=['date','conc.h2o.slope','conc.h2o.intercept','conc.h2o.r2','conc.h2o.p','conc.h2o.sd','conc.h2o.n',\\\n",
    "                                        'conc.co2.slope','conc.co2.intercept','conc.co2.r2','conc.co2.p','conc.co2.sd','conc.co2.n'], index=[0])\n",
    "    temp_cal_df['date'] = temp['timestamp'].dt.date.unique()\n",
    "    if(n_h2o >= 10):\n",
    "        temp_cal_df['conc.h2o.slope']     = h2o_corr.slope\n",
    "        temp_cal_df['conc.h2o.intercept'] = h2o_corr.intercept\n",
    "        temp_cal_df['conc.h2o.r2']        = h2o_corr.rvalue\n",
    "        temp_cal_df['conc.h2o.p']         = h2o_corr.pvalue\n",
    "        temp_cal_df['conc.h2o.sd']        = h2o_corr.stderr\n",
    "    temp_cal_df['conc.h2o.n'] = n_h2o\n",
    "    if(n_co2 >= 10):\n",
    "        temp_cal_df['conc.co2.slope']     = co2_corr.slope\n",
    "        temp_cal_df['conc.co2.intercept'] = co2_corr.intercept\n",
    "        temp_cal_df['conc.co2.r2']        = co2_corr.rvalue\n",
    "        temp_cal_df['conc.co2.p']         = co2_corr.pvalue\n",
    "        temp_cal_df['conc.co2.sd']        = co2_corr.stderr\n",
    "    temp_cal_df['conc.co2.n'] = n_co2\n",
    "    cal_stats = pd.concat([cal_stats, temp_cal_df])\n",
    "        \n",
    "    # Apply corrections\n",
    "    temp['conc.h2o.calibrated.mmol_mol'] = np.nan\n",
    "    temp['conc.co2.calibrated.umol_mol'] = np.nan\n",
    "    if(n_h2o >= 10):\n",
    "        if(h2o_corr.rvalue > 0.75):\n",
    "            temp['conc.h2o.calibrated.mmol_mol'] = h2o_corr.slope * temp['conc.h2o.mmol_mol'] + h2o_corr.intercept\n",
    "    if(n_co2 >= 10):\n",
    "        if(co2_corr.rvalue > 0.50):\n",
    "            temp['conc.co2.calibrated.umol_mol'] = co2_corr.slope * temp['conc.co2.mean.umol_mol'] + co2_corr.intercept\n",
    "    \n",
    "    return(temp, cal_stats)\n",
    "\n",
    "# Writes output files averaged by chamber and time\n",
    "def write_stats_file(out_df, out_dir):\n",
    "    # Check if output folders exist. If not, create\n",
    "    if(not os.path.exists(out_dir)):\n",
    "        os.makedirs(out_dir)\n",
    "    # Create file name\n",
    "    out_fn = out_dir + '/calibration_stats_co2_h2o.csv'\n",
    "    # Remove timezone information\n",
    "    #out_df['timestamp'] = out_df['timestamp'].dt.tz_localize(None)\n",
    "    # If the file does not exist, create it with the header\n",
    "    out_df.to_csv(out_fn, index=False, encoding='utf-8', date_format='%Y-%m-%d')\n",
    "\n",
    "# Writes output files averaged by chamber and time\n",
    "def write_summarised_file(out_df, date_idx, out_dir):\n",
    "    # Check if output folders exist. If not, create\n",
    "    this_month = str(date_idx)[0:6]\n",
    "    if(not os.path.exists(out_dir)): # make directory if it doesn't exist\n",
    "        os.makedirs(out_dir)\n",
    "    # Create file name\n",
    "    out_fn = out_dir + '/' + output_fn + '_' + str(this_month) + '.csv'\n",
    "    # Remove timezone information\n",
    "    #out_df['timestamp'] = out_df['timestamp'].dt.tz_localize(None)\n",
    "    # If the monthly file does not exist, create it with the header\n",
    "    if not os.path.isfile(out_fn):\n",
    "        out_df.to_csv(out_fn, index=False, encoding='utf-8', date_format='%Y-%m-%d %H:%M:%S')\n",
    "    else: # otherwise append without writing the header\n",
    "        out_df.to_csv(out_fn, mode='a', header=False, index=False, encoding='utf-8', date_format='%Y-%m-%d %H:%M:%S')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Full processing pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Applies all the corrections in order. This gets called by the main loops that are in control of reading and writing files\n",
    "def apply_corrections(laser, amb, par, tc, flow, cal_stats, date_idx):\n",
    "    # Show info\n",
    "    print('- Applying corrections on ', np.unique(laser['dayid'])[0], '...')\n",
    "    # Add average columns of the different gases, and standard deviation\n",
    "    laser = fix_gas_units(laser)\n",
    "    # get the right PAR, flow and temperature data\n",
    "    par_df  = get_relevant_data(par, par_log, date_idx, name='par', unit='umol_m2_s1', old_name='par_Avg(', calibration=100)\n",
    "    tc_df   = get_relevant_data(tc, tc_log, date_idx, name='temp', unit='c', old_name='Tc(')\n",
    "    flow_df = get_relevant_data(flow, fm_log, date_idx, name='pump.flow', unit='slpm', old_name='FlowMeter(')\n",
    "    # Merge all data\n",
    "    temp = add_ambient(laser, amb)\n",
    "    temp = add_sensor_data(temp, par_df, interpolation_limit=120) # Interpolation for max. 2min (=120s)\n",
    "    temp = add_sensor_data(temp, tc_df, interpolation_limit=120)    # Interpolation for max. 4 values (=4s)\n",
    "    temp = add_sensor_data(temp, flow_df, interpolation_limit=120)  # Interpolation for max. 4 values (=4s)\n",
    "    temp = fix_PAR(temp) # Fix nighttime PAR data\n",
    "    # Temporarily add counter and chamber\n",
    "    temp = make_counter_and_chamber(temp, date_idx)\n",
    "    # Define day/night\n",
    "    temp = define_day_night(temp)\n",
    "    # Checks if there is enough data left to properly do all the calculations\n",
    "    if(len(temp['timestamp']) > 240):\n",
    "        # Add a column for current PAR, temperature and current flow\n",
    "        temp = add_current_variable(temp)\n",
    "        # Add leaf area\n",
    "        temp = add_surface_area(temp, date_idx, sa_log)\n",
    "        # Add the pump flow rate from the log file, if not previously given by measurements\n",
    "        temp = add_log_flow(temp)\n",
    "        # Remove bad chambers\n",
    "        temp = remove_bad_chambers(temp)\n",
    "        # Remove the data outside the averaging period, if specified\n",
    "        temp = remove_obsolete_data(temp, averaging_period)\n",
    "        # Daily correction of CO2 and H2O by correlating to IRGA values, if available\n",
    "        temp, cal_stats = calibrate_co2_h2o(temp, cal_stats)\n",
    "    else:\n",
    "        print('- Data of', date_idx, 'too short for corrections')\n",
    "        \n",
    "    return(temp, cal_stats)\n",
    "\n",
    "def calculate_values(temp, date_idx):\n",
    "    # Make means, standard deviations, and\n",
    "    # remove all except last number of seconds, provided in the variable 'averaging_period'\n",
    "    temp = make_means(temp, averaging_period)\n",
    "    # Find the corresponding ambient measurements\n",
    "    if(find_ambient == 'interpolated'):\n",
    "        temp = find_corresponding_ambient_interpolated(temp)\n",
    "    elif(find_ambient == 'next'):\n",
    "        temp = find_corresponding_ambient_next(temp)\n",
    "    else:\n",
    "        raise ValueError('Invalid choice: Please choose which method is used to find ambient data')\n",
    "    # calculate fluxes\n",
    "    temp = calculate_fluxes(temp)\n",
    "    # Calculate conductances\n",
    "    temp = calculate_conductances(temp)\n",
    "    \n",
    "    return(temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run calculations\n",
    "#-----------------\n",
    "\n",
    "# List all files in the directories\n",
    "laser_fn_list = sorted(glob.glob(project_path_laser + '*/*', recursive=True))\n",
    "amb_fn_list   = sorted(glob.glob(project_path_amb + '*/*', recursive=True))\n",
    "par_fn_list   = sorted(glob.glob(project_path_par + '*/*', recursive=True))\n",
    "tc_fn_list    = sorted(glob.glob(project_path_tc + '*/*', recursive=True))\n",
    "flow_fn_list  = sorted(glob.glob(project_path_flow + '*/*', recursive=True))\n",
    "saved = []\n",
    "\n",
    "# Prepare 1h directory. Because files aren't overwritten typically, we need to empty it first\n",
    "empty_dir(project_path_output + find_ambient + '/' + '1h/')\n",
    "empty_dir(project_path_output + find_ambient + '/' + '1min/')\n",
    "empty_dir(project_path_output + find_ambient + '/' + 'stats_calibration/')\n",
    "\n",
    "# Prepare empty dataframe for calibration statistics\n",
    "cal_stats_df = pd.DataFrame(columns=['date','conc.h2o.slope','conc.h2o.intercept','conc.h2o.r2','conc.h2o.p','conc.h2o.sd','conc.h2o.n',\\\n",
    "                                     'conc.co2.slope','conc.co2.intercept','conc.co2.r2','conc.co2.p','conc.co2.sd','conc.co2.n'])\n",
    "\n",
    "for fn_i, fn in enumerate(laser_fn_list):\n",
    "    # Only run data in the month list\n",
    "    current_month = fn.replace(project_path_laser[:-1], '')[1:8]\n",
    "    if((current_month in month_list) or (len(month_list) == 0)):\n",
    "        pass\n",
    "    else:\n",
    "        continue\n",
    "    \n",
    "    # Determine current day\n",
    "    this_day = fn[-12:-4]\n",
    "    \n",
    "    # Debugging message\n",
    "    if(fn_i % 1 == 0): # % 20 to show every 20th file being loaded\n",
    "        print( '{:<07}'.format(str(round(fn_i * 100 / len(laser_fn_list), 4))) + '%\\t' + str(this_day[0:4]) + '-' + str(this_day[4:6]) + '-' + str(this_day[6:8])) # Show status\n",
    "        \n",
    "    # Now find corresponding files and load all the data\n",
    "    laser_df = read_data_file(fn)\n",
    "    amb_df  = find_and_load_file(amb_fn_list, this_day)\n",
    "    par_df  = find_and_load_file(par_fn_list, this_day)\n",
    "    tc_df   = find_and_load_file(tc_fn_list, this_day)\n",
    "    flow_df = find_and_load_file(flow_fn_list, this_day)\n",
    "    \n",
    "    # Now apply corrections\n",
    "    out_df, cal_stats_df = apply_corrections(laser_df, amb_df, par_df, tc_df, flow_df, cal_stats_df, this_day)\n",
    "    \n",
    "    # If there is any data left after corrections, i.e. it wasn't all deleted due to bad data\n",
    "    if(len(out_df) > 0):\n",
    "        # Write output in full temporal resolution\n",
    "        write_output_file(out_df, this_day, project_path_output + find_ambient + '/' + 'raw_qc/')\n",
    "        \n",
    "        # Create means and calculate fluxes for the averaging period\n",
    "        out_df = calculate_values(out_df, this_day)\n",
    "        \n",
    "        # Write output averaged output file\n",
    "        write_summarised_file(out_df, this_day, project_path_output + find_ambient + '/' + '1min/')\n",
    "    pass\n",
    "\n",
    "write_stats_file(cal_stats_df, project_path_output + find_ambient + '/' + 'stats_calibration/')\n",
    "print('Done...')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagnostics graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from plotnine import *\n",
    "from mizani.breaks import date_breaks\n",
    "from mizani.formatters import date_format\n",
    "# Colours\n",
    "cbPalette = ['#000000', '#E69F00', '#0072B2', '#CC00CC', '#009E73', '#D55E00', '#CC79A7', '#FF3300', '#F0E442', '#56B4E9']\n",
    "\n",
    "plt = ggplot(out_df)\n",
    "plt = plt + geom_point(aes(x='conc.h2o.mmol_mol.ambient', y='conc.h2o.irga.ambient.mmol_mol', colour='plot'), size=0.5)\n",
    "plt = plt + labs(x='Laser $H_2 O$', y='IRGA $H_2 O$', colour='Type', parse=True)\n",
    "plt = plt + scale_colour_manual(values=cbPalette)\n",
    "plt = plt + theme_bw()\n",
    "plt = plt + theme(text=element_text(family='serif'), axis_text_x=element_text(rotation=90, hjust=0.5))\n",
    "print(plt)\n",
    "\n",
    "plt = ggplot(out_df)\n",
    "plt = plt + geom_point(aes(x='conc.co2.mean.umol_mol.ambient', y='conc.co2.irga.ambient.umol_mol', colour='plot'), size=0.5)\n",
    "plt = plt + labs(x='Laser $CO_2$', y='IRGA $CO_2$', colour='Type', parse=True)\n",
    "plt = plt + scale_colour_manual(values=cbPalette)\n",
    "plt = plt + theme_bw()\n",
    "plt = plt + theme(text=element_text(family='serif'), axis_text_x=element_text(rotation=90, hjust=0.5))\n",
    "print(plt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dp",
   "language": "python",
   "name": "dp"
  },
  "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}