CRNS_COSMOS_Europe_Inversion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79783aca",
   "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",
    "Authors:\n",
    "Roland Baatz roland.baatz @ zalf.de\n",
    "\n",
    "Maintainers and contact:\n",
    "Currently maintained by the authors.\n",
    "\n",
    "This file requires the data set of the COSMOS Europe publication:\n",
    "https://essd.copernicus.org/articles/14/1125/2022/essd-14-1125-2022.html\n",
    "\n",
    "Copyright (C) Leibniz Centre for Agricultural Landscape Research (ZALF)\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 warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\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",
    "\n",
    "# Function to mask values outside the 1st and 99th percentiles\n",
    "def mask_outliers(series):\n",
    "    lower = np.percentile(series, 1)\n",
    "    upper = np.percentile(series, 99)\n",
    "    return series.where((series >= lower) & (series <= upper), np.nan)\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 cosmic ray 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 cosmic ray 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 cosmic ray 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 cosmic ray 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 cosmic 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 cosmic 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\")) \n"
   ]
  },
  {
   "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",
    "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['mean_N'] = np.nan\n",
    "dfi['mean_P'] = np.nan\n",
    "dfi['number_of_days_for_regression'] = np.nan\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95f2ecf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Iterate over EU Study sites:\n",
    "plog_flag_id=0\n",
    "suffix_csv = \"rmse_u\"\n",
    "#csv_filename='my_dataframe_moderated_nelder_rmse.csv'\n",
    "csv_filename=f'my_dataframe_moderated_{suffix_csv}-.csv'\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",
    "    for ii,nn in enumerate(sites_names):\n",
    "        #plog_flag_id=1\n",
    "        #for ii, nn in enumerate(sites_names[3:5], start=3):\n",
    "        print(\"ID:\",ii, \"Reading:\",nn)\n",
    "        with open(\"your_file.txt\", \"a\") as file:\n",
    "            file.write(f'ID: {ii} Site: {nn}\\n')\n",
    "\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",
    "        DF.index.name = 'Date'\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",
    "            \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",
    "        #print(DF.index)\n",
    "        #print(DF.head(2))\n",
    "        #print(Absolute_conv(80,25))\n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Incoming\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric=suffix_csv)\n",
    "        dfi.loc[ii, 'beta'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi'] = my_pars[2]\n",
    "        dfi.loc[ii, 'beta_sd'] = my_pars[3]\n",
    "        dfi.loc[ii, 'omega_sd'] = my_pars[4]\n",
    "        dfi.loc[ii, 'psi_sd'] = my_pars[5]\n",
    "        print(my_pars[8])\n",
    "\n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Apty\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric=suffix_csv)\n",
    "        dfi.loc[ii, 'beta_apty'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega_apty'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi_apty'] = my_pars[2]\n",
    "\n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Oulu\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric=suffix_csv)\n",
    "        dfi.loc[ii, 'beta_oulu'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega_oulu'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi_oulu'] = my_pars[2]\n",
    "\n",
    "        \n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Mxco\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric=suffix_csv)\n",
    "        dfi.loc[ii, 'beta_mxco'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega_mxco'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi_mxco'] = my_pars[2]\n",
    "\n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Lmks\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric=suffix_csv)\n",
    "        dfi.loc[ii, 'beta_lmks'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega_lmks'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi_lmks'] = my_pars[2]\n",
    "\n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Athens\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric=suffix_csv)\n",
    "        dfi.loc[ii, 'beta_athens'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega_athens'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi_athens'] = my_pars[2]\n",
    "\n",
    "\n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Incoming\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric='mae')\n",
    "        dfi.loc[ii, 'beta_mae'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega_mae'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi_mae'] = my_pars[2]\n",
    "        #dfi.loc[ii, 'beta_sd'] = my_pars[3]\n",
    "        #dfi.loc[ii, 'omega_sd'] = my_pars[4]\n",
    "        #dfi.loc[ii, 'psi_sd'] = my_pars[5]\n",
    "\n",
    "        \n",
    "        my_pars=CRNS_opt.parameter_estimator(DF['Pressure'],DF[\"MOD\"],DF[\"Incoming\"],DF[\"Abs_h\"],\"exp\",plot_flag=plog_flag_id,site_name=nn,error_metric=\"rmse\")\n",
    "        dfi.loc[ii, 'beta_rmse'] = my_pars[0]\n",
    "        dfi.loc[ii, 'omega_rmse'] = my_pars[1]\n",
    "        dfi.loc[ii, 'psi_rmse'] = my_pars[2]\n",
    "\n",
    "        \n",
    "        mean_P_McJD = my_pars[6]\n",
    "        number_of_days_for_regression = my_pars[7]\n",
    "        dfi.loc[ii, 'number_of_days_for_regression'] = number_of_days_for_regression\n",
    "        dfi.loc[ii, 'mean_N'] = my_pars[8]\n",
    "        #print(my_pars[8])\n",
    "\n",
    "        g, x, f_lat, f_bar, F_scale, beta_McJD_eff, beta_McJD = CRNS_lib.beta_e(mean_P_McJD,CutoffRigidity[ii],sites_lat[ii])\n",
    "        \n",
    "        dfi.loc[ii, 'beta_McJD_eff'] = beta_McJD_eff\n",
    "        dfi.loc[ii, 'beta_McJD'] = beta_McJD\n",
    "        dfi.loc[ii, 'mean_P'] = mean_P_McJD\n",
    "\n",
    "    dfi.to_csv(csv_filename, sep=\";\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef3c1db",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (my_env_for_swn_kernel)",
   "language": "python",
   "name": "my_env_for_swn_kernel"
  },
  "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}