Figure_7.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef5cb19d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "This Source Code Form is subject to the terms of the Mozilla Public\n",
    "License, v. 2.0. If a copy of the MPL was not distributed with this\n",
    "file, You can obtain one at http://mozilla.org/MPL/2.0/.\n",
    "\n",
    "\n",
    "Authors:\n",
    "Roland Baatz roland.baatz @ zalf.de\n",
    "Data:\n",
    "Paolo Nasta paolo.nasta @ unina.it\n",
    "\n",
    "Maintainers and contact:\n",
    "Currently maintained by the authors.\n",
    "\n",
    "This file requires output of the CRNS_*.ipynb\n",
    "\n",
    "Copyright (C) Leibniz Centre for Agricultural Landscape Research (ZALF)\n",
    "\"\"\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e45a361d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "current_dir = os.getcwd()\n",
    "lib_dir = os.path.join(current_dir, 'lib')\n",
    "sys.path.append(lib_dir)\n",
    "\n",
    "import numpy as np\n",
    "import CRNS_library as CRNS_lib\n",
    "import time\n",
    "import scipy.io as sio\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import CRNS_optimizer as CRNS_opt\n",
    "from scipy.optimize import curve_fit\n",
    "import matplotlib.dates as mdates\n",
    "\n",
    "# Function to calculate RMSE\n",
    "def calculate_rmse(predictions, targets):\n",
    "    return np.sqrt(((predictions - targets) ** 2).mean())\n",
    "\n",
    "\n",
    "\n",
    "file_path_15cm = './data/Alento/Alento2_15cm_precip.txt'\n",
    "file_path_5cm = './data/Alento/Alento2_5cm.txt'\n",
    "alento_site_start=4\n",
    "alento_site_end=5\n",
    "#alento_site_start=12\n",
    "#alento_site_end=13\n",
    "\n",
    "def func(x, p1, p2):\n",
    "    return p1 * np.cos(p2 * x) + p2 * np.sin(p1 * x)\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75c1888",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "src_path = './data'\n",
    "##Read in data for one site as example:\n",
    "#read incoming neutron intensity\n",
    "file_path=src_path+'/Scaling/JUNG1_Data_UTC.TXT'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Incoming'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Jung = new_data.set_index('DateTime')\n",
    "\n",
    "##Read in data for one site as example:\n",
    "#read incoming neutron intensity\n",
    "file_path=src_path+'/Scaling/OULU_1h_UTC.txt'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Oulu'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Oulu = new_data.set_index('DateTime')\n",
    "\n",
    "\n",
    "##Read in data for one site as example:\n",
    "#read incoming neutron intensity\n",
    "file_path=src_path+'/Scaling/LMKS.txt'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Lmks'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Lmks = new_data.set_index('DateTime')\n",
    "\n",
    "##Read in data for one site as example:\n",
    "#read incoming neutron intensity\n",
    "file_path=src_path+'/Scaling/APTY_1h_UTC.txt'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Apty'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Apty = new_data.set_index('DateTime')\n",
    "\n",
    "\n",
    "\n",
    "##Read in data for one site as example:\n",
    "#read incoming neutron ray intensity\n",
    "file_path=src_path+'/Scaling/AthensR8.53Alt260 m.txt'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Athens'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Athens = new_data.set_index('DateTime')\n",
    "\n",
    "\n",
    "##Read in data for one site as example:\n",
    "#read incoming neutron ray intensity\n",
    "file_path=src_path+'/Scaling/MXCO_1h_UTC8.28Alt2274m.txt'\n",
    "new_data = pd.read_csv(file_path, sep=';', names=['DateTime', 'Mxco'], parse_dates=['DateTime'])\n",
    "new_data['DateTime'] = pd.to_datetime(new_data['DateTime'])\n",
    "#print(new_data.head())\n",
    "Mxco = new_data.set_index('DateTime')\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fefe8fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Schäfertal ID 3 removed; and columns renamed; asterisk remarks removed; BUC001 data gaps were removed\n",
    "fname='EU_Study_sites.csv'\n",
    "fname=src_path+\"/Scaling/COSMOS_Europe_Data/Additional_information.csv\"\n",
    "site_info = pd.read_csv(fname)\n",
    "fname=src_path+\"/Scaling/COSMOS_Europe_Data/General_information.csv\"\n",
    "\n",
    "site_general_info = pd.read_csv(fname)\n",
    "sites_names = site_info.StationID.to_list()\n",
    "CutoffRigidity = site_info.CutoffRigidity.to_list()\n",
    "sites_Bulkdensity = site_info.Bulkdensity.to_list()\n",
    "sites_Altitude = site_general_info.Altitude.to_list()\n",
    "sites_Precip = site_general_info.Meanannualprecipitation.to_list()\n",
    "sites_Temperature = site_general_info.Meanairtemperature.to_list()\n",
    "sites_lat = site_general_info.Latitude.to_list()\n",
    "dfi = pd.DataFrame({\n",
    "    'StationID': sites_names,\n",
    "    'CutoffRigidity': CutoffRigidity,\n",
    "    'sites_Precip': sites_Precip,\n",
    "    'sites_Temperature': sites_Temperature,\n",
    "    'Bulkdensity': sites_Bulkdensity,\n",
    "    'Altitude': sites_Altitude\n",
    "})\n",
    "\n",
    "# Add new columns (beta, omega, psi) with initial values (you can update these later)\n",
    "dfi['beta'] = np.nan\n",
    "dfi['beta_sd'] = np.nan\n",
    "dfi['beta_mae'] = np.nan\n",
    "dfi['beta_rmse'] = np.nan\n",
    "dfi['beta_oulu'] = np.nan\n",
    "dfi['beta_apty'] = np.nan\n",
    "dfi['beta_lmks'] = np.nan\n",
    "dfi['beta_mxco'] = np.nan\n",
    "dfi['beta_athens'] = np.nan\n",
    "dfi['omega'] = np.nan\n",
    "dfi['omega_sd'] = np.nan\n",
    "dfi['omega_oulu'] = np.nan\n",
    "dfi['omega_apty'] = np.nan\n",
    "dfi['omega_mxco'] = np.nan\n",
    "dfi['omega_athens'] = np.nan\n",
    "dfi['omega_lmks'] = np.nan\n",
    "dfi['omega_mae'] = np.nan\n",
    "dfi['omega_rmse'] = np.nan\n",
    "dfi['psi'] = np.nan\n",
    "dfi['psi_sd'] = np.nan\n",
    "dfi['psi_oulu'] = np.nan\n",
    "dfi['psi_apty'] = np.nan\n",
    "dfi['psi_athens'] = np.nan\n",
    "dfi['psi_mxco'] = np.nan\n",
    "dfi['psi_lmks'] = np.nan\n",
    "dfi['psi_mae'] = np.nan\n",
    "dfi['psi_rmse'] = np.nan\n",
    "dfi['BetaDesilets'] = np.nan\n",
    "dfi['mean_P'] = np.nan\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95f2ecf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "### Iterate over EU Study sites\n",
    "#Read CRNS and calibrate\n",
    "plog_flag_id=0\n",
    "src_path = './data'\n",
    "suffix_csv = \"rmse_u\"\n",
    "csv_filename=src_path+\"my_dataframe_moderated_\"+suffix_csv+'.csv'\n",
    "\n",
    "\n",
    "if os.path.exists(csv_filename):\n",
    "    dfi = pd.read_csv(csv_filename, sep=';', index_col=0)  # Set index_col to 0 for row names\n",
    "    #LOAD DATAFRAME FROM CSV\n",
    "else:\n",
    "    print(\"ERROR File not present\")\n",
    "    exithere()\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5210b17",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "beta = dfi.loc[alento_site_start, 'beta']\n",
    "omega = dfi.loc[alento_site_start, 'omega']\n",
    "psi = dfi.loc[alento_site_start, 'psi']\n",
    "SiteName_Alento = dfi.loc[alento_site_start,'StationID']\n",
    "print(\"Beta=\", beta)\n",
    "print(\"omega=\",omega)\n",
    "print(\"psi=\",psi)\n",
    "print(\"Name=\",SiteName_Alento)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b4102bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "for ii, nn in enumerate(sites_names[alento_site_start:alento_site_end], start=alento_site_start):\n",
    "    print(\"ID:\",ii, \"Reading:\",nn)\n",
    "    #if nn in ['BUC001']:\n",
    "    #    pass\n",
    "    #else:\n",
    "    #print(ii, nn,sites[ii]) #ID, site and fname, same\n",
    "    DF = CRNS_lib.load_EU_data(src_path+'/Scaling/COSMOS_Europe_Data/',nn)\n",
    "    DF['Abs_h'] = CRNS_opt.Absolute_conv(DF.RH, DF.TEMP)\n",
    "    #resample to 1 hour exactly, for merging with 1h inc. radiation\n",
    "    DF = DF.resample('1H').mean()\n",
    "\n",
    "    #print(DF.tail(3))\n",
    "    DF = DF.join(Jung)\n",
    "    DF = DF.join(Apty)\n",
    "    DF = DF.join(Oulu)\n",
    "    DF = DF.join(Athens)\n",
    "    DF = DF.join(Mxco)\n",
    "    DF = DF.join(Lmks)\n",
    "    #print(DF.head(3))\n",
    "    #print(\"\")\n",
    "    #print(DF.tail(3))\n",
    "    #calculate relative scaling for inc. radiation\n",
    "    \n",
    "    # Normalizing and masking outliers for each column\n",
    "    prctl_values = [0.1, 0.1, 0.1, 0.4, 0.1, 0.1]\n",
    "    columns = [\"Incoming\", \"Apty\", \"Oulu\", \"Mxco\", \"Lmks\", \"Athens\"]\n",
    "    for idcol, col in enumerate(columns):\n",
    "        lower = 0+prctl_values[idcol]\n",
    "        upper=100-prctl_values[idcol]\n",
    "        lower, upper = np.percentile(DF[col].dropna(), [lower,upper])\n",
    "        DF[col] = DF[col].mask((DF[col] < lower) | (DF[col] > upper))\n",
    "    \n",
    "    DF[\"Incoming\"] = DF[\"Incoming\"]/np.nanmean(DF[\"Incoming\"])\n",
    "    DF[\"Apty\"] = DF[\"Apty\"]/np.nanmean(DF[\"Apty\"])\n",
    "    DF[\"Oulu\"] = DF[\"Oulu\"]/np.nanmean(DF[\"Oulu\"])\n",
    "    DF[\"Mxco\"] = DF[\"Mxco\"]/np.nanmean(DF[\"Mxco\"])\n",
    "    DF[\"Lmks\"] = DF[\"Lmks\"]/np.nanmean(DF[\"Lmks\"])\n",
    "    DF[\"Athens\"] = DF[\"Athens\"]/np.nanmean(DF[\"Athens\"])\n",
    "    \n",
    "    mydf = pd.concat([DF['Pressure'],DF[\"MOD\"],DF[\"Incoming\"],DF[\"Abs_h\"]], axis=1)\n",
    "    mydf.columns = ['Pressure', 'MOD', 'Incoming', 'Abs_h']\n",
    "\n",
    "    \"\"\"\n",
    "    #check for ranges and remove outliers\n",
    "    \"\"\"\n",
    "\n",
    "    print(\"Len before removing values out of range:\", len(mydf))\n",
    "    #print(mydf.tail(48))\n",
    "    mydf.loc[mydf['MOD'] < 50, 'MOD'] = np.nan\n",
    "    mydf.loc[mydf['MOD'] > 6000, 'MOD'] = np.nan\n",
    "    mydf.loc[mydf['Pressure'] < 400, 'Pressure'] = np.nan\n",
    "    mydf.loc[mydf['Pressure'] > 1200, 'Pressure'] = np.nan\n",
    "    mydf.loc[mydf['Incoming'] < 0.5, 'Incoming'] = np.nan\n",
    "    mydf.loc[mydf['Incoming'] > 1.5, 'Incoming'] = np.nan\n",
    "    mydf.loc[mydf['Abs_h'] < 0, 'Abs_h'] = np.nan\n",
    "    mydf.loc[mydf['Abs_h'] > 35, 'Abs_h'] = np.nan\n",
    "    \n",
    "    mydf=mydf.dropna(subset=['MOD'])\n",
    "    mydf=mydf.dropna(subset=['Pressure'])\n",
    "    mydf=mydf.dropna(subset=['Incoming'])\n",
    "    mydf=mydf.dropna(subset=['Abs_h'])\n",
    "\n",
    "    #print(mydf.tail(48))\n",
    "    print(\"Len after removing values out of range:\", len(mydf))\n",
    "\n",
    "    \"\"\"\n",
    "    #check for nearest neighboor and remove those with neutron count difference larger 10% to nn+/-1\n",
    "    \"\"\"\n",
    "    mydf['MOD_diff_minus_1'] = abs(mydf.MOD.pct_change(periods=-1)) # difference in neutron count\n",
    "    mydf['MOD_diff_plus_1'] = abs(mydf.MOD.pct_change(periods=1)) # difference in neutron count\n",
    "    mydf = mydf[(mydf['MOD_diff_minus_1']  <= 0.10)&(mydf['MOD_diff_plus_1']  <= 0.10)]\n",
    "\n",
    "\n",
    "    mydf['MOD_diff_minus_1'] = abs(mydf.Pressure.pct_change(periods=-1)) # difference in air pressure\n",
    "    mydf['MOD_diff_plus_1'] = abs(mydf.Pressure.pct_change(periods=1)) # difference in air pressure\n",
    "    mydf = mydf[(mydf['MOD_diff_minus_1']  <= 0.10)&(mydf['MOD_diff_plus_1']  <= 0.10)]\n",
    "    print(\"Len after removing values MOD_diff_plus_1 and MOD_diff_minus_1:\", len(mydf))\n",
    "\n",
    "    mydf['MOD_diff_minus_1'] = abs(mydf.Abs_h.pct_change(periods=-1)) # difference in air h\n",
    "    mydf['MOD_diff_plus_1'] = abs(mydf.Abs_h.pct_change(periods=1)) # difference in air h\n",
    "    mydf = mydf[(mydf['MOD_diff_minus_1']  <= 0.10)&(mydf['MOD_diff_plus_1']  <= 0.10)]\n",
    "    \n",
    "    print(\"Len after removing values MOD_diff_plus_1 and MOD_diff_minus_1:\", len(mydf))\n",
    "    #print(mydf.tail(10))\n",
    "    #print(mydf.head(3))\n",
    "    \n",
    "    daily_df = mydf.resample('D').mean()\n",
    "    print(\"daily_df Min Max and SD and NA.ct are       :\", np.nanmin(daily_df), np.nanmax(daily_df), np.nanstd(daily_df),np.count_nonzero(np.isnan(daily_df)))\n",
    "    \n",
    "    daily_df=daily_df.dropna(subset=['MOD'])\n",
    "    daily_df=daily_df.dropna(subset=['Pressure'])\n",
    "    daily_df=daily_df.dropna(subset=['Incoming'])\n",
    "    daily_df=daily_df.dropna(subset=['Abs_h'])\n",
    "    \n",
    "    #print(daily_df.tail(10))\n",
    "    print(\"daily_df Min Max and SD and NA.ct NA.drop   :\", np.nanmin(daily_df), np.nanmax(daily_df), np.nanstd(daily_df),np.count_nonzero(np.isnan(daily_df)))\n",
    "    # Step 2: Calculate daily differences\n",
    "    # Prepare data for optimization\n",
    "    mod_values = daily_df['MOD'].values\n",
    "    pressure_values = daily_df['Pressure'].values\n",
    "    incoming_values = daily_df['Incoming'].values\n",
    "    abs_h_values = daily_df['Abs_h'].values\n",
    "print(\"\\n daily_df head:\")\n",
    "print(daily_df.head(3))\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e61c97a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Read the CSV file\n",
    "data = pd.read_csv(file_path_15cm, delim_whitespace=True)\n",
    "# Step 2: Group by Year, Month, and Day\n",
    "grouped = data.groupby(['Year', 'Month', 'Day'])\n",
    "# Step 3: Filter out days with 2 or more missing entries\n",
    "valid_days = grouped.filter(lambda x: x.shape[0] >= 22)\n",
    "# Step 4: Compute daily averages\n",
    "daily_avg_swc = valid_days.groupby(['Year', 'Month', 'Day']).mean().reset_index()\n",
    "\n",
    "# Step 1: Read the CSV file\n",
    "data = pd.read_csv(file_path_5cm, delim_whitespace=True)\n",
    "# Step 2: Group by Year, Month, and Day\n",
    "grouped = data.groupby(['Year', 'Month', 'Day'])\n",
    "# Step 3: Filter out days with 2 or more missing entries\n",
    "valid_days = grouped.filter(lambda x: x.shape[0] >= 22)\n",
    "# Step 4: Compute daily averages\n",
    "daily_avg_swc_5cm = valid_days.groupby(['Year', 'Month', 'Day']).mean().reset_index()\n",
    "\n",
    "###Merge two dataframes for where data is present\n",
    "df_daily_avg_swc = pd.DataFrame(daily_avg_swc)\n",
    "df_daily_avg_swc['Date'] = pd.to_datetime(df_daily_avg_swc[['Year', 'Month', 'Day']])\n",
    "df_daily_avg_swc.drop(['Year', 'Month', 'Day','Hour'], axis=1, inplace=True)\n",
    "\n",
    "df_daily_avg_swc_5cm = pd.DataFrame(daily_avg_swc_5cm)\n",
    "df_daily_avg_swc_5cm['Date'] = pd.to_datetime(df_daily_avg_swc_5cm[['Year', 'Month', 'Day']])\n",
    "df_daily_avg_swc_5cm.drop(['Year', 'Month', 'Day','Hour'], axis=1, inplace=True)\n",
    "\n",
    "# Merging the DataFrames\n",
    "merged_df = pd.merge(daily_df, df_daily_avg_swc, on='Date', how='inner')\n",
    "merged_df = pd.merge(merged_df, df_daily_avg_swc_5cm, on='Date', how='inner')\n",
    "print(\" HEAD \")\n",
    "print(merged_df.head(3))\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65772c97",
   "metadata": {},
   "outputs": [],
   "source": [
    "### for the merged_df: Identify a 5 day period with 5 days of measurement, where Incoming changes strongly, and Rainfall is 0.\n",
    "#save this period as index.\n",
    "# Plot MOD, Incoming, Abs_h, and average_SM for this period only.\n",
    "\n",
    "myvars = ['Incoming','Abs_h','Pressure']\n",
    "for iiid, myvar in enumerate(myvars):\n",
    "    print(myvar, iiid)\n",
    "    if myvar=='Abs_h':\n",
    "        my_vars_diff = \"Abs_diff\"\n",
    "    else:\n",
    "        my_vars_diff = f'{myvar}_diff'\n",
    "    print(my_vars_diff)\n",
    "    merged_df[my_vars_diff] = merged_df[myvar].diff().abs()\n",
    "    \"\"\"\n",
    "    filtered_df = merged_df.copy()\n",
    "    temp_df = filtered_df[my_vars_diff].fillna(-np.inf)\n",
    "    \n",
    "    # Find the top 3 maximum values and their indices\n",
    "    top_3_values = temp_df.nlargest(3).values\n",
    "    top_3_indices = temp_df.nlargest(3).index\n",
    "    \n",
    "    print(f\"Top 3 change values: {top_3_values}\")\n",
    "    print(f\"Indices of top 3 changes: {top_3_indices}\")\n",
    "\n",
    "    for i, (max_value, max_index) in enumerate(zip(top_3_values, top_3_indices), 1):\n",
    "        print(f\"\\nMax change value {i}: {max_value}\")\n",
    "        print(f\"Index of max change {i}: {max_index}\")\n",
    "\n",
    "        index = filtered_df.index[filtered_df[my_vars_diff] == max_value].tolist()\n",
    "        if index:\n",
    "            extracted_index = index[0]  # Assuming you want the first match if duplicates\n",
    "        else:\n",
    "            extracted_index = None  # Handle the case where no match is found\n",
    "\n",
    "        print(f\"Extracted Index {i}: {extracted_index}\")\n",
    "        print(f\"Extracted Index {i}: {extracted_index}, {filtered_df.loc[extracted_index, my_vars_diff]}\")\n",
    "        print(\" \")\n",
    "\n",
    "        # Display surrounding rows for the extracted index\n",
    "        print(merged_df.loc[(extracted_index-2):(extracted_index+2), ['Date', myvar, my_vars_diff, 'Rainfall', 'SWC_15cm', 'SWC_5cm']])\n",
    "        print(\" \")\n",
    "    \"\"\"\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0f62c56",
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_df = merged_df.copy()\n",
    "printflag=0\n",
    "\n",
    "# Filter for periods where over THREE days: Rainfall == 0\n",
    "inc_df = filtered_df['Incoming_diff'].fillna(-np.inf)\n",
    "absh_df = filtered_df['Abs_diff'].fillna(-np.inf)\n",
    "pressure_df = filtered_df['Pressure_diff'].fillna(-np.inf)\n",
    "\n",
    "# Find the maximum value and its index\n",
    "for ii in range(1,10):\n",
    "    # Get the second largest value in column 'A'\n",
    "    my_n_largest = filtered_df['Incoming_diff'].nlargest(ii).iloc[-1]\n",
    "    my_index=filtered_df['Incoming_diff']==my_n_largest\n",
    "    max_index = inc_df.idxmax()\n",
    "    \n",
    "    index = filtered_df.index[filtered_df['Incoming_diff'] == my_n_largest].tolist()\n",
    "    if index:\n",
    "        extracted_index = index[0]  # Assuming you want the first match\n",
    "    else:\n",
    "        extracted_index = None  # Handle the case where no match is found\n",
    "    tmp_val=filtered_df.loc[extracted_index, 'Incoming_diff']*100\n",
    "    if printflag:\n",
    "        print(\"Extracted Index:\", extracted_index, f'Value = {tmp_val:.1f} %')\n",
    "        print(\" \")\n",
    "        print(merged_df.loc[(extracted_index-2):(extracted_index+2),['Date','Incoming','Incoming_diff','Rainfall','SWC_5cm']])\n",
    "        #print(merged_df.loc[(extracted_index-2):(extracted_index+2),['Date','Incoming','Incoming_diff','Rainfall','SWC_15cm','SWC_5cm']])\n",
    "        print(\" \")\n",
    "\n",
    "    \n",
    "    # Get the second largest value in column 'A'\n",
    "    my_n_largest = filtered_df['Pressure_diff'].nlargest(ii).iloc[-1]\n",
    "    #my_n_largest = filtered_df['Incoming_diff'].nsmallest(ii).iloc[-1]\n",
    "    my_index=filtered_df['Pressure_diff']==my_n_largest\n",
    "    max_index = pressure_df.idxmax()\n",
    "\n",
    "    if printflag: print(\" \")\n",
    "    index = filtered_df.index[filtered_df['Pressure_diff'] == my_n_largest].tolist()\n",
    "    if index:\n",
    "        extracted_index = index[0]  # Assuming you want the first match\n",
    "    else:\n",
    "        extracted_index = None  # Handle the case where no match is found\n",
    "    tmp_val=filtered_df.loc[extracted_index, 'Pressure_diff']\n",
    "    if printflag:\n",
    "        print(\"Extracted Index:\", extracted_index, f'Value = {tmp_val:.1f} hPa')\n",
    "        print(\" \")\n",
    "        print(merged_df.loc[(extracted_index-2):(extracted_index+2),['Date','Pressure','Pressure_diff','Rainfall','SWC_5cm']])\n",
    "        print(\" \")\n",
    "\n",
    "    \n",
    "    # Get the second largest value in column 'A'\n",
    "    my_n_largest = filtered_df['Abs_diff'].nlargest(ii).iloc[-1]\n",
    "    #my_n_largest = filtered_df['Incoming_diff'].nsmallest(ii).iloc[-1]\n",
    "    my_index=filtered_df['Abs_diff']==my_n_largest\n",
    "    max_index = absh_df.idxmax()\n",
    "\n",
    "    if printflag: print(\" \")\n",
    "    index = filtered_df.index[filtered_df['Abs_diff'] == my_n_largest].tolist()\n",
    "    if index:\n",
    "        extracted_index = index[0]  # Assuming you want the first match\n",
    "    else:\n",
    "        extracted_index = None  # Handle the case where no match is found\n",
    "    tmp_val=filtered_df.loc[extracted_index, 'Abs_diff']\n",
    "    if printflag:\n",
    "        print(\"Extracted Index:\", extracted_index, f'Value = {tmp_val:.1f} Abs. Humidity')\n",
    "        print(\" \")\n",
    "        print(merged_df.loc[(extracted_index-2):(extracted_index+2),['Date','Abs_h','Abs_diff','Rainfall','SWC_5cm']])\n",
    "        #print(merged_df.loc[(extracted_index-2):(extracted_index+2),['Date','Incoming','Incoming_diff','Rainfall','SWC_15cm','SWC_5cm']])\n",
    "        print(\" \")\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31879047",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de703a44",
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_df = merged_df.copy()\n",
    "\n",
    "#extended range\n",
    "inc_indices=np.arange(523,532,1)\n",
    "abs_indices=np.arange(502,512,1)\n",
    "pres_indices=np.arange(432,445,1)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "filtered_df = merged_df.copy()\n",
    "\n",
    "beta_ref=-0.0076\n",
    "omega_ref=-0.0054\n",
    "psi_ref=0.8\n",
    "bd = 1.126\n",
    "lw=0.055\n",
    "\n",
    "p_ref=np.nanmean(filtered_df['Pressure'])\n",
    "i_ref=np.nanmean(filtered_df['Incoming'])\n",
    "h_ref=np.nanmean(filtered_df['Abs_h'])\n",
    "\n",
    "beta = dfi.loc[alento_site_start, 'beta']\n",
    "omega = dfi.loc[alento_site_start, 'omega']\n",
    "psi = dfi.loc[alento_site_start, 'psi']\n",
    "\n",
    "SWC_calib = filtered_df.loc[inc_indices[0],['SWC_5cm']].values\n",
    "\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = filtered_df.loc[inc_indices[0],\"MOD_new\"]\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[inc_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = filtered_df.loc[inc_indices[0],\"MOD_standard\"]\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[inc_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "SWC_calib = np.mean(filtered_df.loc[inc_indices,['SWC_5cm']].values)\n",
    "\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = np.mean(filtered_df.loc[inc_indices,\"MOD_new\"])\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[inc_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = np.mean(filtered_df.loc[inc_indices,\"MOD_standard\"])\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[inc_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "# Create a figure with three columns\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "ax1 = axes[0]\n",
    "# Plot the four SWC lines\n",
    "ax1.plot(filtered_df.loc[inc_indices,'Date'], SWC_new,       label='SWC (new approach)', color='red', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[inc_indices,'Date'], SWC_standard,  label='SWC (ref approach)', color='green', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[inc_indices,'Date'], filtered_df.loc[inc_indices, 'SWC_5cm'],  label='SWC ( 5cm TDT)', color='blue', marker='*',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[inc_indices,'Date'], filtered_df.loc[inc_indices, 'SWC_5cm'], label='SWC (observed 5cm)', color='blue', marker='',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[inc_indices,'Date'], filtered_df.loc[inc_indices, 'SWC_15cm'], label='SWC (15cm TDT)', color='purple', linestyle=':')\n",
    "# Add a legend\n",
    "ax1.legend()\n",
    "# Add labels\n",
    "ax1.set_xlabel('Date')\n",
    "ax1.set_ylabel('SWC')\n",
    "# Create a second y-axis for 'Incoming'\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(filtered_df.loc[inc_indices,'Date'], filtered_df.loc[inc_indices, 'Incoming'], label='Incoming', color='grey', linestyle='--')\n",
    "ax2.set_ylabel('Incoming')\n",
    "# Add a legend for the secondary y-axis\n",
    "\n",
    "\n",
    "the_indices = abs_indices\n",
    "SWC_calib = filtered_df.loc[the_indices[0],['SWC_5cm']].values\n",
    "\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = filtered_df.loc[the_indices[0],\"MOD_new\"]\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[the_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = filtered_df.loc[the_indices[0],\"MOD_standard\"]\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[the_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "SWC_calib = np.mean(filtered_df.loc[the_indices,['SWC_5cm']].values)\n",
    "\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_new\"])\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[the_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_standard\"])\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[the_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "ax1 = axes[1]\n",
    "# Plot the four SWC lines\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_new,       label='SWC (new approach)', color='red', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_standard,  label='SWC (ref approach)', color='green', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_5cm'], label='SWC (observed 5cm)', color='blue', marker='*',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_15cm'], label='SWC (15cm TDT)', color='purple', linestyle=':')\n",
    "ax1.legend()\n",
    "# Add labels\n",
    "ax1.set_xlabel('Date')\n",
    "ax1.set_ylabel('SWC')\n",
    "# Create a second y-axis for 'Incoming'\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'Abs_h'], color='grey', linestyle='--')\n",
    "ax2.set_ylabel('Absolute Humidity')\n",
    "# Add a legend for the secondary y-axis\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "the_indices = pres_indices\n",
    "SWC_calib = filtered_df.loc[the_indices[0],['SWC_5cm']].values\n",
    "\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = filtered_df.loc[the_indices[0],\"MOD_new\"]\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[the_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = filtered_df.loc[the_indices[0],\"MOD_standard\"]\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[the_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "SWC_calib = np.mean(filtered_df.loc[the_indices,['SWC_5cm']].values)\n",
    "\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_new\"])\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[the_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_standard\"])\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[the_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "ax1 = axes[2]\n",
    "# Plot the four SWC lines\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_new,       label='SWC (new approach)', color='red', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_standard,  label='SWC (ref approach)', color='green', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_5cm'], label='SWC (observed 5cm)', color='blue', marker='*',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_15cm'], label='SWC (15cm TDT)', color='purple', linestyle=':')\n",
    "ax1.legend()\n",
    "# Add labels\n",
    "ax1.set_xlabel('Date')\n",
    "ax1.set_ylabel('SWC')\n",
    "# Create a second y-axis for 'Incoming'\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'Pressure'], color='grey', linestyle='--',marker='x')\n",
    "ax2.set_ylabel('Air Pressure')\n",
    "# Add a legend for the secondary y-axis\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "fig.tight_layout()\n",
    "\n",
    "# Show plot\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Alento2.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f687b39",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import matplotlib.dates as mdates\n",
    "filtered_df = merged_df.copy()\n",
    "\n",
    "\n",
    "# Convert 'Date' column to datetime format if not already\n",
    "#filtered_df['Date'] = pd.to_datetime(filtered_df['Date'])\n",
    "\n",
    "# Set 'Date' column as the index\n",
    "#filtered_df.set_index('Date', inplace=True)\n",
    "\n",
    "\n",
    "#print(filtered_df.index)\n",
    "# Ensure the index is in datetime format\n",
    "#filtered_df.index = pd.to_datetime(filtered_df.index)#\n",
    "\n",
    "#print(filtered_df.index)\n",
    "\n",
    "beta_ref=-0.0076\n",
    "omega_ref=-0.0054\n",
    "psi_ref=0.8\n",
    "bd = 1.126\n",
    "lw=0.055\n",
    "\n",
    "filtered_df[\"SWC_Weighted\"] = filtered_df[\"SWC_15cm\"].copy()\n",
    "depths = [0.05, 0.15]  # depths in cm\n",
    "depths = [0.15, 0.30]  # depths in cm\n",
    "radial_distance = 0  # assuming no radial distribution\n",
    "print(\"LENGTH:\")\n",
    "print(len(filtered_df))\n",
    "##Calculate weighting:\n",
    "for ii in range(len(filtered_df)):\n",
    "    #print(ii)\n",
    "    swc_values = [filtered_df.loc[ii,\"SWC_5cm\"]+lw, filtered_df.loc[ii,\"SWC_15cm\"]+lw]\n",
    "    # Calculate weighted SWC\n",
    "    weighted_swcs = []\n",
    "    for d, swc in zip(depths, swc_values):\n",
    "        weight = CRNS_lib.Wd(d, radial_distance, bd, swc)\n",
    "        weighted_swcs.append(swc * weight)\n",
    "\n",
    "    # Calculate average weighted SWC\n",
    "    average_weighted_swc = sum(weighted_swcs) / len(weighted_swcs)\n",
    "    #print(f\"Average Weighted SWC: {average_weighted_swc}\")\n",
    "    filtered_df.loc[ii,\"SWC_Weighted\"] = average_weighted_swc-lw\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "p_ref=np.nanmean(filtered_df['Pressure'])\n",
    "i_ref=np.nanmean(filtered_df['Incoming'])\n",
    "h_ref=np.nanmean(filtered_df['Abs_h'])\n",
    "\n",
    "\n",
    "beta = dfi.loc[alento_site_start, 'beta']\n",
    "omega = dfi.loc[alento_site_start, 'omega']\n",
    "psi = dfi.loc[alento_site_start, 'psi']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ffe2fb2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec7f0e70",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Function to calculate RMSE\n",
    "def calculate_rmse(predictions, targets):\n",
    "    return np.sqrt(((predictions - targets) ** 2).mean())\n",
    "\n",
    "\n",
    "filtered_df = merged_df.copy()\n",
    "#print(merged_df.index)\n",
    "#print(merged_df.columns)\n",
    "\n",
    "# Convert 'Date' column to datetime format if not already\n",
    "#filtered_df['Date'] = pd.to_datetime(filtered_df['Date'])\n",
    "\n",
    "# Set 'Date' column as the index\n",
    "#filtered_df.set_index('Date', inplace=True)\n",
    "\n",
    "\n",
    "#print(filtered_df.index)\n",
    "# Ensure the index is in datetime format\n",
    "#filtered_df.index = pd.to_datetime(filtered_df.index)#\n",
    "\n",
    "#print(filtered_df.index)\n",
    "\n",
    "beta_ref=-0.0076\n",
    "omega_ref=-0.0054\n",
    "psi_ref=0.8\n",
    "bd = 1.126\n",
    "lw=0.055\n",
    "\n",
    "filtered_df[\"SWC_Weighted\"] = filtered_df[\"SWC_15cm\"].copy()\n",
    "depths = [0.05, 0.15]  # depths in cm\n",
    "depths = [0.15, 0.30]  # depths in cm\n",
    "radial_distance = 0  # assuming no radial distribution\n",
    "\n",
    "print(len(filtered_df))\n",
    "##Calculate weighting:\n",
    "for ii in range(len(filtered_df)):\n",
    "    #print(ii)\n",
    "    swc_values = [filtered_df.loc[ii,\"SWC_5cm\"]+lw, filtered_df.loc[ii,\"SWC_15cm\"]+lw]\n",
    "    # Calculate weighted SWC\n",
    "    weighted_swcs = []\n",
    "    for d, swc in zip(depths, swc_values):\n",
    "        weight = CRNS_lib.Wd(d, radial_distance, bd, swc)\n",
    "        weighted_swcs.append(swc * weight)\n",
    "\n",
    "    # Calculate average weighted SWC\n",
    "    average_weighted_swc = sum(weighted_swcs) / len(weighted_swcs)\n",
    "    #print(f\"Average Weighted SWC: {average_weighted_swc}\")\n",
    "    filtered_df.loc[ii,\"SWC_Weighted\"] = average_weighted_swc-lw\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "p_ref=np.nanmean(filtered_df['Pressure'])\n",
    "i_ref=np.nanmean(filtered_df['Incoming'])\n",
    "h_ref=np.nanmean(filtered_df['Abs_h'])\n",
    "\n",
    "\n",
    "beta = dfi.loc[alento_site_start, 'beta']\n",
    "omega = dfi.loc[alento_site_start, 'omega']\n",
    "psi = dfi.loc[alento_site_start, 'psi']\n",
    "\n",
    "\n",
    "the_indices = inc_indices.copy()\n",
    "\n",
    "print(\"the_indices\")\n",
    "print(the_indices)\n",
    "print(\"filtered_df.loc[the_indices[0],['SWC_Weighted']].values\")\n",
    "print(filtered_df.loc[the_indices[0],['SWC_Weighted']].values)\n",
    "\n",
    "SWC_calib = np.mean(filtered_df.loc[inc_indices,['SWC_Weighted']].values)\n",
    "\n",
    "print(\"CALIB SWC:\")\n",
    "print(SWC_calib)\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = np.mean(filtered_df.loc[inc_indices,\"MOD_new\"])\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[inc_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = np.mean(filtered_df.loc[inc_indices,\"MOD_standard\"])\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[inc_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "print(np.mean(SWC_new))\n",
    "print(np.mean(SWC_standard))\n",
    "\n",
    "print(the_indices)\n",
    "print(filtered_df.loc[the_indices,'Date'])\n",
    "print(\" \")\n",
    "print(\" \")\n",
    "\n",
    "# Create a figure with three columns\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "ax1 = axes[0]\n",
    "# Assuming 'filtered_df' is your DataFrame and 'the_indices' are the indices of interest\n",
    "# Compute RMSE values\n",
    "rmse_new = calculate_rmse(filtered_df.loc[the_indices, 'SWC_Weighted'], SWC_new)\n",
    "rmse_standard = calculate_rmse(filtered_df.loc[the_indices, 'SWC_Weighted'], SWC_standard)\n",
    "label_new = f'SWC (new approach) RMSE={rmse_new:.3f}'\n",
    "label_standard = f'SWC (ref approach) RMSE={rmse_standard:.3f}'\n",
    "\n",
    "improvement_percent = ((rmse_standard - rmse_new) / rmse_standard) * 100\n",
    "print(\"Improvement incoming neutrons:\", improvement_percent, \"%\")\n",
    "print(\" \")\n",
    "\n",
    "\n",
    "# Plot the four SWC lines\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_new,       label=label_new, color='red', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_standard,  label=label_standard, color='green', linestyle='--')\n",
    "#ax1.plot(filtered_df.loc[inc_indices,'Date'], filtered_df.loc[inc_indices, 'SWC_5cm'],  label='SWC ( 5cm TDT)', color='blue', marker='*',linestyle='-')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_Weighted'],  label='SWC (weighted TDT)', color='blue', marker='*',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[inc_indices,'Date'], filtered_df.loc[inc_indices, 'SWC_5cm'], label='SWC (observed 5cm)', color='blue', marker='',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[inc_indices,'Date'], filtered_df.loc[inc_indices, 'SWC_15cm'], label='SWC (15cm TDT)', color='purple', linestyle=':')\n",
    "ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "# Add a legend\n",
    "ax1.legend()\n",
    "# Add labels\n",
    "ax1.set_xlabel('Date')\n",
    "ax1.set_ylabel('SWC')\n",
    "# Create a second y-axis for 'Incoming'\n",
    "ax2 = ax1.twinx()\n",
    "\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_p[the_indices], color='grey', linestyle=':')\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_i[the_indices], color='grey', linestyle='-')\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_h[the_indices], color='grey', linestyle=':')\n",
    "ax2.set_ylabel('Scaling factors')\n",
    "ax2.set_title('- Incoming neutron intensity dynamics')\n",
    "# Add a legend for the secondary y-axis\n",
    "\n",
    "\n",
    "the_indices = abs_indices.copy()\n",
    "\n",
    "SWC_calib = np.mean(filtered_df.loc[the_indices,['SWC_Weighted']].values)\n",
    "\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_new\"])\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[the_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_standard\"])\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[the_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "ax1 = axes[1]\n",
    "# Assuming 'filtered_df' is your DataFrame and 'the_indices' are the indices of interest\n",
    "# Compute RMSE values\n",
    "rmse_new = calculate_rmse(filtered_df.loc[the_indices, 'SWC_Weighted'], SWC_new)\n",
    "rmse_standard = calculate_rmse(filtered_df.loc[the_indices, 'SWC_Weighted'], SWC_standard)\n",
    "label_new = f'SWC (new approach) RMSE={rmse_new:.3f}'\n",
    "label_standard = f'SWC (ref approach) RMSE={rmse_standard:.3f}'\n",
    "\n",
    "improvement_percent = ((rmse_standard - rmse_new) / rmse_standard) * 100\n",
    "print(\"Improvement absolute air humidity:\", improvement_percent, \"%\")\n",
    "print(\" \")\n",
    "\n",
    "# Plot the four SWC lines\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_new,       label=label_new, color='red', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_standard,  label=label_standard, color='green', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_Weighted'], label='SWC (weighted TDT)', color='blue', marker='*',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_15cm'], label='SWC (15cm TDT)', color='purple', linestyle=':')\n",
    "ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "ax1.legend()\n",
    "# Add labels\n",
    "ax1.set_ylabel('SWC')\n",
    "# Create a second y-axis for 'Incoming'\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_p[the_indices], color='grey', linestyle=':')\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_i[the_indices], color='grey', linestyle=':')\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_h[the_indices], color='grey', linestyle='-')\n",
    "ax2.set_ylabel('Scaling factors')\n",
    "ax2.set_title('- Absolute humidity dynamics')\n",
    "# Add a legend for the secondary y-axis\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "the_indices = pres_indices.copy()\n",
    "\n",
    "SWC_calib = np.mean(filtered_df.loc[the_indices,['SWC_Weighted']].values)\n",
    "cor_p_new = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta)\n",
    "cor_h_new = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega)\n",
    "cor_i_new = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi)\n",
    "\n",
    "filtered_df[\"MOD_new\"] = filtered_df['MOD'] / cor_p_new / cor_i_new / cor_h_new\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_new\"])\n",
    "N0_new = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_new = ((0.0808/((filtered_df.loc[the_indices,\"MOD_new\"]/N0_new)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "cor_p = CRNS_opt.corr_exp_pressure(filtered_df['Pressure'], p_ref, beta_ref)\n",
    "cor_h = CRNS_opt.corr_lin_humidity(filtered_df['Abs_h'], h_ref, omega_ref)\n",
    "cor_i = CRNS_opt.corr_lin_incoming_correction(filtered_df['Incoming'],i_ref, psi_ref)\n",
    "\n",
    "filtered_df[\"MOD_standard\"] = filtered_df['MOD'] / cor_p / cor_i / cor_h\n",
    "Npih_calibrated = np.mean(filtered_df.loc[the_indices,\"MOD_standard\"])\n",
    "N0_standard = Npih_calibrated / (0.0808/((SWC_calib+lw) / bd + 0.115) + 0.372)\n",
    "SWC_standard = ((0.0808/((filtered_df.loc[the_indices,\"MOD_standard\"]/N0_standard)-0.372)-0.115)*bd)-lw\n",
    "\n",
    "\n",
    "\n",
    "ax1 = axes[2]\n",
    "# Assuming 'filtered_df' is your DataFrame and 'the_indices' are the indices of interest\n",
    "# Compute RMSE values\n",
    "rmse_new = calculate_rmse(filtered_df.loc[the_indices, 'SWC_Weighted'], SWC_new)\n",
    "rmse_standard = calculate_rmse(filtered_df.loc[the_indices, 'SWC_Weighted'], SWC_standard)\n",
    "label_new = f'SWC (new approach) RMSE={rmse_new:.3f}'\n",
    "label_standard = f'SWC (ref approach) RMSE={rmse_standard:.3f}'\n",
    "\n",
    "improvement_percent = ((rmse_standard - rmse_new) / rmse_standard) * 100\n",
    "print(\"Improvement atmospheric pressure:\", improvement_percent, \"%\")\n",
    "print(\" \")\n",
    "\n",
    "# Plot the four SWC lines\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_new,       label=label_new, color='red', linestyle='--')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], SWC_standard,  label=label_standard, color='green', linestyle='--')\n",
    "#ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_5cm'], label='SWC (5cm TDT)', color='purple', linestyle='-.')\n",
    "ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_Weighted'], label='SWC (weighted TDT)', color='blue', marker='*',linestyle='-')\n",
    "#ax1.plot(filtered_df.loc[the_indices,'Date'], filtered_df.loc[the_indices, 'SWC_15cm'], label='SWC (15cm TDT)', color='purple', linestyle=':')\n",
    "ax1.legend()\n",
    "# Add labels\n",
    "ax1.set_ylabel('SWC')\n",
    "ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))\n",
    "ax1.tick_params(axis='x', rotation=45)\n",
    "\n",
    "# Create a second y-axis for 'Incoming'\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_p[the_indices], color='grey', linestyle='-')\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_i[the_indices], color='grey', linestyle=':')\n",
    "ax2.plot(filtered_df.loc[the_indices,'Date'], cor_h[the_indices], color='grey', linestyle=':')\n",
    "ax2.set_ylabel('Scaling factors')\n",
    "ax2.set_title('- Air pressure dynamics')\n",
    "# Add a legend for the secondary y-axis\n",
    "# Adjust layout to prevent overlap\n",
    "fig.tight_layout()\n",
    "\n",
    "# Show plot\n",
    "plt.show()\n",
    "\n",
    "\n",
    "pngname=\"./plots/Alento2_scaling_plotted.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "500b2336",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdac0d3c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fcab386",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9 (py_env_lda)",
   "language": "python",
   "name": "py_env_lda"
  },
  "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": 5
}