Data-Analisys.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a934b755-b362-42a9-87c7-e3e1674720d9",
   "metadata": {},
   "source": [
    "# **Data analisys for validation MODIS vs VESPA-22 IWV measurements**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "488471b6-adbf-4a2e-a31c-29417919d63c",
   "metadata": {},
   "source": [
    "## **Import libraries**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4cc8665-3aa3-454c-90be-832818ee0c32",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from patsy import dmatrices\n",
    "\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7ac9551-6957-459e-b2e8-04ed088d935c",
   "metadata": {},
   "source": [
    "## **Selected data reception from Thule-mix-LAADS.ipynb**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "419221e3-eafc-467e-a331-0c4bdcb2df25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# passaggio parametri di filtraggio dati\n",
    "%store -r max_IWV\n",
    "%store -r max_dist_from_THAAO\n",
    "%store -r deltatime\n",
    "# passaggio dati MODIS\n",
    "%store -r new_mod_vespa\n",
    "%store -r new_myd_vespa\n",
    "%store -r new_mxd_vespa\n",
    "# passaggio plot_defaults\n",
    "%store -r plot_defaults\n",
    "# passaggio dati completi non filtrati con deltatime\n",
    "%store -r vespa_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "378d74dc-d1e8-4f45-85af-cc95a68ed763",
   "metadata": {},
   "source": [
    "### **Controllo dati ricevuti**"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7b2cb90d-9048-49a3-922e-e441bdbdb87e",
   "metadata": {},
   "source": [
    "new_mod_vespa.head(1)     # ok: raw per test"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d92b8b61-ed6e-4e2c-97c6-aac0313bb4de",
   "metadata": {},
   "source": [
    "new_myd_vespa.head(1)     # ok: raw per test"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9ac2da4a-b101-42b9-a65a-533cb377aee7",
   "metadata": {
    "tags": []
   },
   "source": [
    "new_mxd_vespa.head(1)     # ok: raw per test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29d9f64b-15a4-4434-afca-f85a70002f40",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_data.head(1)     # ok: raw per test\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a81d19e9-ee76-4e71-bdac-44da35876c11",
   "metadata": {
    "tags": []
   },
   "source": [
    "# converto il time range selezionato da datetime a timestamp (float in seconds, default=ns)\n",
    "vespa_time_range['timestamp'] = pd.to_numeric(vespa_time_range['Time'].values)/ 10 ** 9  # to timestamp in seconds (default=ns)\n",
    "type(vespa_time_range['timestamp'][2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3064a27c-efe6-4d26-8c68-277dfd4a337a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# converto da datetime a timestamp (float in seconds, default=ns)\n",
    "vespa_data['timestamp'] = pd.to_numeric(vespa_data['Time'].values)/ 10 ** 9  # to timestamp in seconds (default=ns)\n",
    "type(vespa_data['timestamp'][2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad34435a-1811-409c-bb9c-b8cea548e7d7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# ok: raw per test\n",
    "print(type(vespa_data['Time'][0]))\n",
    "vespa_data.head(1)    "
   ]
  },
  {
   "cell_type": "raw",
   "id": "2c40fde1-d035-45d7-9d57-3a726080df0a",
   "metadata": {
    "tags": []
   },
   "source": [
    "from datetime import datetime\n",
    "#dt = datetime.today()  # Get timezone naive now\n",
    "#seconds = dt.timestamp()\n",
    "vespa_data['time'] = datetime.timestamp(vespa_data['Time'])\n",
    "vespa_data.head(1)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf01ad9a-5ab4-4a8f-aa7e-aa1957dfeb91",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_mod_data_p_Dt = len(new_mod_vespa.IWV_MODIS) \n",
    "num_myd_data_p_Dt = len(new_myd_vespa.IWV_MODIS)\n",
    "num_mxd_data_p_Dt = len(new_mxd_vespa.IWV_MODIS)\n",
    "num_modmyd_data_p_Dt = num_mod_data_p_Dt + num_myd_data_p_Dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d78ab45-ae27-4a41-a7f8-3b693c4997f9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(f'-- Selection parameters --')\n",
    "print(f'max IWV = {max_IWV} [mm]')\n",
    "print(f'max distance from THAOO = {max_dist_from_THAAO} [km]')\n",
    "print(f'Deltatime chosen for data points selection = {deltatime}\\n')\n",
    "print(f'-- Data points selected --')\n",
    "print(f'Selected number of MOD data points =          {num_mod_data_p_Dt}')\n",
    "print(f'Selected number of MYD data points =          {num_myd_data_p_Dt}\\n')\n",
    "print(f'Total selected number of MXD data points     =  {num_mxd_data_p_Dt}')\n",
    "print(f'Total selected number of MOD+MYD data points =  {num_modmyd_data_p_Dt}\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8eb5b2a-175b-4494-add3-7ec961c70014",
   "metadata": {},
   "source": [
    "## **Analisi spettrale**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "051dab1e-b2b0-4c93-9041-accdb483eae3",
   "metadata": {
    "tags": []
   },
   "source": [
    "### **Selezione e plot dei dati in un intervallo temporale** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f99676fd-63a6-4e93-a05d-b6d9689a4f8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.pandas\n",
    "\n",
    "color_THAAO = 'green'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d11f1df9-42d2-4859-a4cb-e1ae460ea289",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# seleziono un intervallo temporale che voglio analizzare \n",
    "# NOTA: .copy() è necessario per evitare una shallow copy ed il settingwithcopywarning\n",
    "vespa_time_range = vespa_data.loc[vespa_data[\"Time\"].between('2019-01-01 00:00:00','2021-01-01 00:00:00'),:].copy()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "47d1a34b-b412-4d8c-9068-b7a2eeb092e4",
   "metadata": {},
   "source": [
    "vespa_time_range.tail(2)  #head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb3a9ea8-06ee-4921-83d2-25668a417e15",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# test selezione\n",
    "print(len(vespa_time_range))\n",
    "vespa_time_range.head(3)  #.tail(3)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "787a5cef-401f-4ed1-b3d2-a2a167d09c1b",
   "metadata": {
    "tags": []
   },
   "source": [
    "# riconverto il time range selezionato da timestamp (float in seconds) a datetime \n",
    "vespa_time_range['datetime'] = pd.to_datetime(vespa_time_range['timestamp'], \n",
    "                                              utc=True, unit='s')\n",
    "type(vespa_time_range['datetime'][2].value)\n",
    "vespa_time_range.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c21d8aa3-b130-46f1-bce9-60f5725b5de5",
   "metadata": {},
   "outputs": [],
   "source": [
    "vespa_ts_scatter = vespa_time_range.IWV_THAAO.hvplot.scatter(color=color_THAAO, label='VESPA',\n",
    "                                                       fields={'date': 'Date'},\n",
    "                                                       title='IWV THAAO - VESPA', \n",
    "                                                       **plot_defaults)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "24169a6f-8ab9-43cf-895e-af4245188f4c",
   "metadata": {},
   "source": [
    "vespa_ts_scatter_rec = vespa_time_range.IWV_THAAO.hvplot.scatter(color='red', label='VESPA_rec',\n",
    "                                                       fields={'date': 'Date'},\n",
    "                                                       title='IWV THAAO - VESPA', \n",
    "                                                       **plot_defaults)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22532776-1ef6-4cac-b2ac-7aff117bce11",
   "metadata": {},
   "outputs": [],
   "source": [
    "#(vespa_ts_scatter*vespa_ts_scatter_rec).opts(show_grid=True)\n",
    "(vespa_ts_scatter).opts(show_grid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5766bab9-1b36-4fd2-a249-a5a07d3d43ed",
   "metadata": {},
   "source": [
    "### **Interpolazione**\n",
    "Poiché i dati non sono equispaziati e contengono buchi provo ad interpolare prima di fare lo spettro <br>\n",
    "\n",
    "vedi [gaps](http://mres.uni-potsdam.de/index.php/2017/08/22/data-voids-and-spectral-analysis-dont-be-afraid-of-gaps/) <br>\n",
    "**NOTA**: Since all these (and other) interpolation techniques might introduce artifacts into the data, <br> \n",
    "it is always advisable to (1) keep the total number of data points constant before and after <br>\n",
    "interpolation, (2) report the method employed for estimating the evenly-spaced data sequence, <br>\n",
    "and (3) explore the effect of interpolation on the variance of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b95d72a6-8fa4-4886-96ab-22f043fbbe45",
   "metadata": {},
   "outputs": [],
   "source": [
    "# creo la time series di timestamp con dati temporalmente equispaziati"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a608b3e3-baf5-4058-b22c-ea0524a72686",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(vespa_time_range['Time'][0])\n",
    "print(vespa_time_range['Time'].iloc[-1])\n",
    "print(len(vespa_time_range['Time']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc9144c6-1767-4cd2-b074-dae8bb31da5e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "date_rng = pd.date_range(start=vespa_time_range['Time'][0], \n",
    "                         end=vespa_time_range['Time'].iloc[-1], \n",
    "                         periods=len(vespa_time_range['Time']))\n",
    "print(len(date_rng))\n",
    "print(date_rng[0])\n",
    "print(date_rng[-1])\n",
    "print(type(date_rng))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "681b9067-e656-40e1-8423-d90543139dbd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_time_range['new_date_rng'] = date_rng"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f79fe086-f522-484c-90aa-c062aea08226",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_time_range.tail(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e3c0333-d08e-4ae0-8707-d0e49bd005a5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# converto il time range selezionato da dateime a timestamp (float in seconds, default=ns)\n",
    "vespa_time_range['new_timestamp_rng'] = pd.to_numeric(vespa_time_range['new_date_rng'].values)/ 10 ** 9  # to timestamp in seconds (default=ns)\n",
    "vespa_time_range.tail(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f45f892-5f5a-4367-80b7-cbfd6b26f5e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_time_range = vespa_time_range.set_index(['Time','new_date_rng'], \n",
    "                            drop=False).sort_index()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "711522ff-32e7-4415-8773-ef282a924165",
   "metadata": {
    "tags": []
   },
   "source": [
    "vespa_time_range"
   ]
  },
  {
   "cell_type": "raw",
   "id": "04e7bee1-9fc9-41ef-a7de-c28f65a84a86",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "raw",
   "id": "c08e5534-f271-4463-8aca-96d37f3b3679",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b7edc87a-4a5e-485c-b9bf-03a6716b69b0",
   "metadata": {},
   "source": [
    "#### **Linear, Cubic spline interpolation, Piecewise Cubic Hermite Interpolating Polynomial interpolation, Akima1D**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c023493-7774-4ac4-a54e-285ddf0b855e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d, CubicSpline, PchipInterpolator, Akima1DInterpolator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be4d0d6e-4d4b-4b4b-b2a3-d7e1747d8728",
   "metadata": {},
   "outputs": [],
   "source": [
    "linear = interp1d(vespa_time_range['timestamp'], vespa_time_range['IWV_THAAO'], kind='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87da62a0-6443-40a1-8e74-efab894850b7",
   "metadata": {},
   "source": [
    "Interpolo sui dati reali con le date datetime convertite in timestamp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a166c841-a8a7-4639-b3c2-00e71527848e",
   "metadata": {},
   "outputs": [],
   "source": [
    "CS = CubicSpline(vespa_time_range['timestamp'],\n",
    "                 vespa_time_range['IWV_THAAO'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcfba4f4-8912-4cf5-acd4-35fb0eca478e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PCHIP = PchipInterpolator(vespa_time_range['timestamp'], \n",
    "                          vespa_time_range['IWV_THAAO'], \n",
    "                          axis=0, extrapolate=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50452af7-cd0f-49aa-ad17-bcbf4b9f1a62",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Akima1D = Akima1DInterpolator(vespa_time_range['timestamp'],\n",
    "                              vespa_time_range['IWV_THAAO'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7f97024-ac00-4e93-8658-16d1ccabf797",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_time_range['linear'] =  linear(vespa_time_range['new_timestamp_rng'])\n",
    "vespa_time_range['CS'] =  CS(vespa_time_range['new_timestamp_rng'])\n",
    "vespa_time_range['PCHIP'] =  PCHIP(vespa_time_range['new_timestamp_rng'])\n",
    "vespa_time_range['Akima1D'] =  Akima1D(vespa_time_range['new_timestamp_rng'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35b3ca17-1f4b-4b73-b611-aee997a5f6cf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_time_range.loc[vespa_time_range['Time']=='2020-01-08 15:56:28']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "165a4553-5b20-4cdc-b063-59b090fd6766",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_time_range['new_timestamp_rng'][1]-vespa_time_range['new_timestamp_rng'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f52b9fc5-2fff-464a-aec7-d7ef7158d6d9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "type(vespa_time_range.new_date_rng[0])\n",
    "#type(vespa_time_range.Time[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93ab5aed-d29a-4407-8382-d2a5ff541686",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vespa_linear = vespa_time_range.linear.hvplot.scatter(x='new_date_rng',\n",
    "                                                      color='blue',\n",
    "                                                      label='linear',\n",
    "                                                      title='IWV THAAO - VESPA',\n",
    "                                                      **plot_defaults)\n",
    "vespa_CS = vespa_time_range.CS.hvplot.scatter(x='new_date_rng',\n",
    "                                              color='red',\n",
    "                                              label='CS',\n",
    "                                              title='IWV THAAO - VESPA',\n",
    "                                              **plot_defaults)\n",
    "vespa_PCHIP = vespa_time_range.PCHIP.hvplot.scatter(x='new_date_rng', \n",
    "                                                    color='cyan', \n",
    "                                                    label='PCHIP',\n",
    "                                                    title='IWV THAAO - VESPA', \n",
    "                                                    **plot_defaults)\n",
    "vespa_Akima1D = vespa_time_range.Akima1D.hvplot.scatter(x='new_date_rng', \n",
    "                                                        color='yellow', \n",
    "                                                        label='Akima1D',\n",
    "                                                        title='IWV THAAO - VESPA', \n",
    "                                                        **plot_defaults)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f13ceac2-a625-4705-b91d-4788907540c7",
   "metadata": {
    "tags": []
   },
   "source": [
    "# grafico vespa nel range temporale definito\n",
    "vespa_ts_scatter"
   ]
  },
  {
   "cell_type": "raw",
   "id": "8f33cc3d-cfa6-4c7b-a9e3-9d4f5acaf3d4",
   "metadata": {
    "tags": []
   },
   "source": [
    "# grafico interpolazione\n",
    "vespa_PCHIP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5bd8d86-d0d4-4086-b613-6193383c0d4a",
   "metadata": {},
   "source": [
    "#### **Plot of interpolations and original data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f3845af-57e3-421b-8fa8-f8d504a59483",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# grafico vespa nel range temporale definito ed interpolazione\n",
    "(vespa_ts_scatter\n",
    "*vespa_linear\n",
    "#*vespa_CS\n",
    "*vespa_PCHIP\n",
    "#*vespa_Akima1D\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7030f9c-ef8b-4bf9-8ccb-1928509a457c",
   "metadata": {
    "tags": []
   },
   "source": [
    "### **Periodogrammi**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d792ecf-8fea-48d4-bf73-db6029d343bf",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### **Dati equispaziati (interpolati)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f87f7a96-f01c-40ce-8d77-2420953b9357",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Passo di campionamento dt (sample spacing, secondi tra un dato ricampionato e il successivo)\n",
    "dt = vespa_time_range['new_timestamp_rng'][10002]-vespa_time_range['new_timestamp_rng'][10001]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c83e9ff5-1ca3-49c9-a05b-3d4a874fcd36",
   "metadata": {},
   "source": [
    "##### **FFT**"
   ]
  },
  {
   "cell_type": "raw",
   "id": "0db86ba3-8718-41ce-927a-0fc619268086",
   "metadata": {
    "tags": []
   },
   "source": [
    "#Fourier Analysis in Python\n",
    "import numpy as np\n",
    "from scipy.fft import fft\n",
    "dt = 4502.40588\n",
    "# compute PSD using simple FFT\n",
    "N = len(vespa_time_range.PCHIP)\n",
    "print(N)\n",
    "df = 1. / (N * dt)\n",
    "PSD = abs(dt*fft((vespa_time_range.PCHIP).values)[:N / 2])**2\n",
    "f = df * np.arange(N / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f8854cf-fa63-4986-a5f5-4b25f4eb2ecf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Fourier Analysis in Python\n",
    "from scipy.fft import fft, fftfreq\n",
    "import numpy as np\n",
    "\n",
    "# Number of sample points\n",
    "N = len(vespa_time_range)\n",
    "\n",
    "x = (vespa_time_range['new_timestamp_rng']).values #np.linspace(0.0, N*T, N, endpoint=False)\n",
    "y = (vespa_time_range['PCHIP']).values #np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)\n",
    "yf = fft(y)\n",
    "xf = fftfreq(N, dt)[0:N//2]\n",
    "\n",
    "print(len(xf), len(2.0/N * np.abs(yf[0:N//2])), len(yf), ' N =',N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a7901f1-c6a3-43af-a06e-41f1661eda9b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Nota unità di misura: xf_s [s**(-1)], xf_d [days**(-1)] \n",
    "#print(xf[1], 2.0/N * np.abs(yf[1//2]))\n",
    "print(len(xf), len(2.0/N * np.abs(yf[0:N//2])),  ' N = ',N )\n",
    "print(type(xf), type(2.0/N * np.abs(yf[0:N//2])))\n",
    "FFT_all_f = pd.DataFrame({'xf_s':xf, 'xf_d':xf*86400, 'yf':2.0/N * np.abs(yf[0:N//2])}) \n",
    "#FFT_all_f.set_index('xf_s',drop=True,inplace=True)\n",
    "FFT_all_f.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fccf1c0-e1d6-4a29-9ad7-fdb863f5cc9f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# drop first row to plot log(x)\n",
    "FFT = FFT_all_f.iloc[1:].copy()\n",
    "print(len(FFT))\n",
    "FFT['xt_d'] = (1./FFT.xf_d)  # unità di misura: xt_d [days] \n",
    "FFT.set_index(['xf_s','xf_d'],drop=True,inplace=True)\n",
    "f_max_gg = FFT.yf*86400\n",
    "FFT.head(2)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "67d2abed-2c91-4129-afcb-305499cdb9f2",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "source": [
    "hvplot.help('line')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0fbca2f-eecc-4370-a913-bd8f968c8dea",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import hvplot.pandas\n",
    "from bokeh.plotting import figure, output_file, show, ColumnDataSource\n",
    "from bokeh.models import HoverTool,WheelZoomTool, PanTool, ResetTool\n",
    "#hvplot.extension('bokeh')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "47b9da8a-8ef6-4dd9-afd2-f689bb916a63",
   "metadata": {
    "tags": []
   },
   "source": [
    "# %opts Line [tools=[hover] hover_cols=['xf_s','yf','xt_d'] grid=True]\n",
    "# %opts Line [grid=True]\n",
    "#tools = [PanTool(), WheelZoomTool(), ResetTool(), hover]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88bbedb3-9876-489c-8125-f48f2014d91b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# tooltips cambia solo il nome nel hovertool\n",
    "hover = HoverTool(tooltips=[('f_seconds','@xf_s'),('t_days','@xt_d'),('FFT','@yf')])\n",
    "FFT.hvplot.line(x='xf_s', y='yf', logx=True, grid=True, hover_cols=['xf_s','yf','xt_d']).opts(tools=[hover])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "dfe97e8a-54c6-444e-8205-2594c565bcfe",
   "metadata": {
    "tags": []
   },
   "source": [
    "FFT.hvplot.line(x='xf_d', y='yf', logx=True, grid=True, hover_cols=['xf_s','yf','xt_d']).opts(tools=[hover])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "876a912f-70b0-4178-a500-bdeefcd7b9ea",
   "metadata": {
    "tags": []
   },
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "raw",
   "id": "ca6e10cf-eedd-4c26-8435-1ec6ee176d3d",
   "metadata": {
    "tags": []
   },
   "source": [
    "#plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))\n",
    "plt.semilogx(xf, 2.0/N * np.abs(yf[0:N//2]))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab8899f7-f5fe-4b2c-9e44-0d1ff81d418b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(1/dt)\n",
    "f_max = 3*(10**-8)     # Hz\n",
    "f_max_gg = f_max*86400      # 1/giorni\n",
    "f_max_2 = 6.1*(10**-8)     # Hz\n",
    "f_max_gg_2 = f_max_2*86400      # 1/giorni\n",
    "print(f'f_max={f_max} [Hz], f_max_gg={f_max} [day-1]')\n",
    "TT=(1./(f_max_gg))\n",
    "print(f'TT = {TT} [days]')\n",
    "\n",
    "print(f'f_max_2={f_max_2} [Hz], f_max_gg_2={f_max_2} [day-1]')\n",
    "TT_2=(1./(f_max_gg_2))\n",
    "print(f'TT_2 = {TT_2} [days]')\n",
    "#TT/86400 # a okkio 579 giorni il picco più grande"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fc108aa-6aae-4963-9d5e-54e452d1fa2c",
   "metadata": {},
   "source": [
    "##### **PS and PSD** \n",
    "Power Spectrum e Power Spectral Density base <br> \n",
    "(scipy.signal.periodogram(data,fs,scaling='spectrum') <br>\n",
    "unità di misura: <br>\n",
    "IWV: $[mm^2]$ <br>\n",
    "frequency: $[Hz]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "416b346e-0711-4238-b114-a631e8bf4de7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "fs = 1/dt \n",
    "T = len(vespa_time_range)*dt  # T=N*dt [s]\n",
    "df = 1/T  # df=1/T=fs/N , dim bin in freq. df = (1/(N*dt))[Hz]=fatt. di conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e6fb6ba-0d1b-4c3d-9246-fdceca834ab3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f_ps, PS_f = signal.periodogram(vespa_time_range.PCHIP,fs, \n",
    "                                window='blackman', #'boxcar', #'hann', \n",
    "                                scaling='spectrum')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee059f54-56b5-4a05-9f04-fdc6936af8d6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f_psd, PSD_f = signal.periodogram(vespa_time_range.PCHIP,fs, scaling='density')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "514a93d7-bf78-4b2f-b906-3055f2729712",
   "metadata": {
    "tags": []
   },
   "source": [
    "freq = f_ps-f_psd\n",
    "print(min(freq),max(freq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5ec58c7-d555-4118-863d-566ca8114959",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# print di controllo\n",
    "print(f'len(vespa_time_range) = {len(vespa_time_range)}')\n",
    "print(f'len(PS)  = {len(PS_f)}, len(f_ps)  = {len(f_ps)}')\n",
    "print(f'len(PSD) = {len(PSD_f)}, len(f_psd) = {len(f_ps)}\\n')\n",
    "\n",
    "print(f'dt = {dt} [s], fs = {fs} [s^-1]')\n",
    "print(f'dt_o = {dt/3600} [hours]')\n",
    "print(f'T   = {T} [s], df = {df} 1/(N*dt))[Hz]')\n",
    "print(f'T_a = {T/(365*86400)} [years]\\n')\n",
    "\n",
    "print(f'min(PS) = {min(PS_f)}, max(PS) = {max(PS_f)}')\n",
    "print(f'min(PSD) = {min(PSD_f)}, max(PSD) = {max(PSD_f)}')\n",
    "print(f'min(f_ps)  = {min(f_ps)}, max(f_ps)  = {max(f_ps)}')\n",
    "print(f'min(f_psd) = {min(f_psd)}, max(f_psd) = {max(f_psd)}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e567b697-2cca-4616-a0f1-ca4ca34b5c0d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PS_PSD_f = pd.DataFrame({'f_ps':f_ps, 'f_ps_d':f_ps*86400, 'PS':PS_f,\n",
    "                        'f_psd':f_psd, 'f_psd_d':f_psd*86400, 'PSD':PSD_f}) #2.0/N * np.abs(yf[0:N//2])}) \n",
    "PS_PSD_f.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7afb7bba-ba9e-4173-825e-9275007e1744",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# drop first row to plot log(x)\n",
    "PS_PSD = PS_PSD_f.iloc[1:].copy()\n",
    "print(len(f_ps))\n",
    "PS_PSD['t_ps_d'] = (1./PS_PSD.f_ps_d)  # unità di misura: xt_d [days] \n",
    "PS_PSD['t_psd_d'] = (1./PS_PSD.f_psd_d)  # unità di misura: xt_d [days] \n",
    "#PS_PSD.set_index(['f_ps','f_ps_d','f_psd','f_psd_d'],drop=True,inplace=True)\n",
    "PS_PSD.set_index(['f_ps','f_psd'],drop=True,inplace=True)\n",
    "PS_PSD.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c33ebfe3-946f-4a09-b285-cf1899c227b6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# tooltips cambia il nome nel hovertool\n",
    "hover = HoverTool(tooltips=[('f_ps','@f_ps'),('t_ps_d','@t_ps_d'),('PS','@PS')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b7e41cc-a324-4094-9c50-5e7eda01700b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PS_PSD.hvplot.line(x='f_ps', y='PS', logx=True, logy=True, grid=True, \n",
    "                   xlabel='frequency [Hz]', ylabel='PS [mm**2]', \n",
    "                   hover_cols=['f_ps','PS','t_ps_d']).opts(tools=[hover])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf1d2440-6065-47ac-a521-df3b53ba6af2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# tooltips cambia il nome nel hovertool\n",
    "hover = HoverTool(tooltips=[('f_psd','@f_psd'),('t_psd_d','@t_psd_d'),('PSD','@PSD')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f45e58f4-0f81-49e5-b926-9234f2ddfd7c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PS_PSD.hvplot.line(x='f_psd', y='PSD', logx=True, logy=True, grid=True, \n",
    "                   xlabel='frequency [Hz]', ylabel='PSD [mm**2/Hz]', \n",
    "                   hover_cols=['f_psd','PSD','t_psd_d']).opts(tools=[hover])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "924f564d-1ce9-48a6-809d-3d99e4e76f0f",
   "metadata": {},
   "source": [
    "##### **PSD Welch** \n",
    "Power Spectral Density Welch <br>\n",
    "(scipy.signal.welch(data, fs, nperseg=1024, scaling='spectrum') <br>\n",
    "unità di misura: <br>\n",
    "IWV: $[mm^2]$ <br>\n",
    "frequency: $[Hz]$"
   ]
  },
  {
   "cell_type": "raw",
   "id": "59178d15-9350-4a42-b279-22de805598ad",
   "metadata": {
    "tags": []
   },
   "source": [
    "f_w, PSD_welch = signal.welch(vespa_time_range.PCHIP, fs)#, nperseg=4096)\n",
    "len(PSD_welch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc5f6d22-793e-45db-8ef5-958bfa4f84de",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f_w, PSD_welch = signal.welch(vespa_time_range.PCHIP, fs,\n",
    "                              #scaling='spectrum',\n",
    "#                              window='boxcar', \n",
    "#                              window='hann', \n",
    "#                              window='blackman', \n",
    "                              window=('kaiser', 14), \n",
    "                              nperseg=len(vespa_time_range.PCHIP))\n",
    "len(PSD_welch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9de2f35a-53b5-4063-aeaa-d7edf780faaa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PSD_wel = pd.DataFrame({'f_w':f_w, 'f_w_d':f_w*86400, 'PSD_welch':PSD_welch}) \n",
    "print(PSD_wel.head(2))\n",
    "# drop first row to plot log(x)\n",
    "PSD_w = PSD_wel.iloc[1:].copy()\n",
    "print(len(f_ps))\n",
    "PSD_w['t_w_d'] = (1./PSD_w.f_w_d)  # unità di misura: xt_d [days] \n",
    "PSD_w.set_index(['f_w'],drop=True,inplace=True)\n",
    "PSD_w.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7539c572-d743-4ce3-ac81-4c1ee266bd24",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# tooltips cambia il nome nel hovertool\n",
    "hover = HoverTool(tooltips=[('f_w','@f_w'),('t_w_d','@t_w_d'),('PSD_welch','@PSD_welch')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdb7f90d-c08a-416c-8d74-dd634f1b0760",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PSD_w.hvplot.line(x='f_w', y='PSD_welch', logx=True, logy=True, grid=True, \n",
    "                   xlabel='frequency [Hz]', ylabel='PSD Welch [mm**2/Hz]', \n",
    "                   hover_cols=['f_w','PSD_welch','t_w_d']).opts(tools=[hover])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "2b7cbf98-2538-4f34-b23d-ae3fc9e7cd28",
   "metadata": {
    "tags": []
   },
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#plt.semilogy(f_p, Pxx_den)\n",
    "plt.loglog(f_ps, PS_f)\n",
    "plt.xlim([1e-8, 2*1e-4]) #0.00014])\n",
    "plt.ylim([1e-9, 1e1])\n",
    "plt.xlabel('frequency [Hz]')\n",
    "plt.ylabel('PS [mm**2]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "0c8e56bc-44ae-4c76-b454-4e3ab09dd305",
   "metadata": {
    "tags": []
   },
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#plt.semilogy(f_p, Pxx_den)\n",
    "plt.loglog(f_psd, PSD)\n",
    "plt.xlim([1e-8, 2*1e-4]) #0.00014])\n",
    "plt.ylim([1e-2, 1e10])\n",
    "plt.xlabel('frequency [Hz]')\n",
    "plt.ylabel('PSD [mm**2/Hz]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "ee081b6d-09a7-4fe7-8edd-7412b9fad283",
   "metadata": {
    "tags": []
   },
   "source": [
    "vespa_remove_mean = vespa_time_range.PCHIP #- vespa_time_range.PCHIP.mean()\n",
    "f_welch, Pwelch_spec = signal.welch(vespa_remove_mean, fs, nperseg=8192)#, scaling='spectrum')\n",
    "#f_welch, Pwelch_spec = signal.welch(vespa_time_range.PCHIP, fs, scaling='spectrum')\n",
    "#plt.loglog(f_welch, Pwelch_spec)\n",
    "\n",
    "plt.loglog(f_w, PSD_welch)\n",
    "plt.xlabel('frequency [Hz]')\n",
    "plt.ylabel('PSD')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d36d5a84-08bb-4f83-86dd-90dd9a54dbbb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Pxx_den.hvplot.line(x='f', y='Pxx_den')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6b37fe-0c0b-42e8-acea-2e58d54ac311",
   "metadata": {},
   "source": [
    "#### **Error/Trend/Seasonality (ETS) model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc54f88e-0ed4-49d3-9097-3cfce6991bb3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from statsmodels.tsa.seasonal import seasonal_decompose\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d95fbcb4-40a0-4de3-a8a5-b9c01095b8c7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "new_mxd_vespa_1 = new_mxd_vespa.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8aba119b-30eb-4ea2-be67-9b4d1650403a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "new_mxd_vespa_1 = new_mxd_vespa_1.set_index('Time')\n",
    "new_mxd_vespa_1 = new_mxd_vespa_1.sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d97ad28b-26d7-42a0-9493-290d22ad6e5d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#new_mxd_vespa_1 = new_mxd_vespa_1.asfreq(pd.infer_freq(new_mxd_vespa_1.index),method='bfill')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ff11cad-976a-46f8-8638-fb20695672f7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#new_mxd_vespa_1.asfreq(method='bfill')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "50b07c6f-14f4-474e-99fa-04884e453ad3",
   "metadata": {
    "tags": []
   },
   "source": [
    "result_add = seasonal_decompose(new_mxd_vespa_1['IWV_MODIS'],\n",
    "                                #period=pd.to_datetime(new_mxd_vespa['Time_IWV_MODIS']).astype(np.int64),\n",
    "                                period=500,\n",
    "                                model = 'add')\n",
    "\n",
    "result_add.plot();\n",
    "#result_mul = seasonal_decompose(df['Sales'],model = 'mul')\n",
    "#result_mul.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "920842ac-59d3-4ed0-b010-289f73c17f7c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "len(vespa_time_range.PCHIP)\n",
    "vespa_time_range.head(2)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "89db939d-7d04-43e5-b3bb-f83317db3e67",
   "metadata": {
    "tags": []
   },
   "source": [
    "fs = 1/dt \n",
    "T = len(vespa_time_range)*dt  # T=N*dt [s]\n",
    "df = 1/T  # df=1/T=fs/N , dim bin in freq. df = (1/(N*dt))[Hz]=fatt. di conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb3f9b66-2768-4d6c-9974-ee4f2b8b0b7d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# print di controllo\n",
    "print(f'len(vespa_time_range) = {len(vespa_time_range)}')\n",
    "print(f'len(PS)  = {len(PS_f)}, len(f_ps)  = {len(f_ps)}')\n",
    "print(f'len(PSD) = {len(PSD_f)}, len(f_psd) = {len(f_ps)}\\n')\n",
    "\n",
    "print(f'dt = {dt} [s], fs = {fs} [s^-1]')\n",
    "print(f'dt_o = {dt/3600} [hours]')\n",
    "print(f'T   = {T} [s], df = {df} 1/(N*dt))[Hz]')\n",
    "print(f'T_a = {T/(365*86400)} [years]\\n')\n",
    "\n",
    "print(f'min(PS) = {min(PS_f)}, max(PS) = {max(PS_f)}')\n",
    "print(f'min(PSD) = {min(PSD_f)}, max(PSD) = {max(PSD_f)}')\n",
    "print(f'min(f_ps)  = {min(f_ps)}, max(f_ps)  = {max(f_ps)}')\n",
    "print(f'min(f_psd) = {min(f_psd)}, max(f_psd) = {max(f_psd)}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "217c0baa-4226-420e-9b30-74467bce0f9c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# ricreo un df pandas con i dataset che mi interessano \n",
    "vespa_ETS_timerange = (vespa_time_range[['new_date_rng', 'new_timestamp_rng','PCHIP']].copy()).droplevel(0).drop('new_date_rng', axis=1)\n",
    "#vespa_ETS_tr.index[100:200].mean()\n",
    "vespa_ETS_timerange.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2553da4-0c50-494a-96ee-17065e81c9eb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "T_delta = vespa_ETS_timerange.new_timestamp_rng[101]-vespa_ETS_timerange.new_timestamp_rng[100]\n",
    "print(T_delta)\n",
    "periodo = int(365*86400/T_delta)\n",
    "periodo "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56ee9dd5-3384-41d8-92d1-6240be3facfb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "result_add = seasonal_decompose(vespa_ETS_timerange.PCHIP,\n",
    "                                #period=T_delta.seconds,\n",
    "                                period=periodo,\n",
    "                                model = 'add', extrapolate_trend='freq')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e7735a8-c3a8-45bd-942f-c6491ea074c3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#print(f'trend={result_add.trend},\\n\\n seasonal={result_add.seasonal},\\n\\n resid={result_add.resid}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f742c4f4-4d3a-4d3d-b3cd-5084d238a00a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "result_add.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6713de59-23d5-4802-860a-2d05917a33a8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from statsmodels.tsa.seasonal import STL\n",
    "\n",
    "stl = STL(vespa_ETS_timerange.PCHIP, period=periodo, seasonal=13)\n",
    "res = stl.fit()\n",
    "fig = res.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd61b677-d64c-409d-8b1b-d1f676d709bc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(vespa_ETS_timerange.index.freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6b9a80c-81f7-4509-98ae-3df14bf65991",
   "metadata": {},
   "source": [
    "#### **Lomb-Scargle periodogram for unevenly spaced data**"
   ]
  },
  {
   "cell_type": "raw",
   "id": "20c0dd1e-095a-4512-ae4a-2bf828d70f07",
   "metadata": {},
   "source": [
    "import scipy.signal as signal\n",
    "pgram = signal.lombscargle(x, y, w, normalize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "768ca197-0970-4e90-93e6-581f7e9ac160",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from astropy.timeseries import LombScargle\n",
    "import astropy.units as u\n",
    "\n",
    "dy = 0.05          # 5% error on VESPA measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff6c4aa6-7cd3-4895-9c0b-b11de0bdaf08",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "frequency, power = LombScargle(vespa_time_range['timestamp'],\n",
    "                               vespa_time_range['IWV_THAAO'], \n",
    "                               dy#, \n",
    "                               # normalization='psd'\n",
    "                              ).autopower()\n",
    "minimum_frequency = min(frequency)\n",
    "print(f'min(frequency) = {min(frequency)}, max(frequency) = {max(frequency)}')\n",
    "print(f'min(power) = {min(power)}, max(power) = {max(power)}')\n",
    "print(f'min(vespa_time_range.IWV_THAAO) = {min(vespa_time_range.IWV_THAAO.values)}')\n",
    "print(f'max(vespa_time_range.IWV_THAAO) = {max(vespa_time_range.IWV_THAAO.values)}')\n",
    "print(f'min(vespa_time_range.timestamp) = {min(vespa_time_range.timestamp.values)}')\n",
    "print(f'max(vespa_time_range.timestamp) = {max(vespa_time_range.timestamp.values)}\\n')\n",
    "print(f'len(vespa_time_range) = {len(vespa_time_range)}')\n",
    "print(f'len(power) = {len(power)}')\n",
    "print(f'len(frequency) = {len(frequency)}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35f88b28-d95c-4a74-bfe6-541b94c58c7b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# tolgo i valori nan dalla colonna IWV\n",
    "print(vespa_data['IWV_THAAO'].isnull().values.any())\n",
    "vespa_dropna = vespa_data.dropna(subset = ['IWV_THAAO'])\n",
    "print(vespa_dropna['IWV_THAAO'].isnull().values.any())"
   ]
  },
  {
   "cell_type": "raw",
   "id": "01aab9f1-e0ac-405c-81a7-617ee293b429",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "source": [
    "frequency, power = LombScargle(vespa_time_range['timestamp'], \n",
    "                               vespa_time_range['IWV_THAAO'], \n",
    "                               dy#, \n",
    "                               # normalization='psd'\n",
    "                              ).autopower(minimum_frequency=minimum_frequency,\n",
    "                                          maximum_frequency=1.e-6,\n",
    "                                          samples_per_peak=100)\n",
    "print(f'min(frequency) = {min(frequency)}, max(frequency) = {max(frequency)}')\n",
    "print(f'min(power) = {min(power)}, max(power) = {max(power)}')\n",
    "print(f'min(vespa_time_range.IWV_THAAO) = {min(vespa_time_range.IWV_THAAO.values)}')\n",
    "print(f'max(vespa_time_range.IWV_THAAO) = {max(vespa_time_range.IWV_THAAO.values)}')\n",
    "print(f'len(frequency) = {len(frequency)}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94c87a49-09e9-4240-962c-428ede87d0fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "frequency, power = LombScargle(vespa_dropna['timestamp'], \n",
    "                               vespa_dropna['IWV_THAAO'], \n",
    "                               dy#,\n",
    "                               #fit_mean=True#,\n",
    "                               # normalization='psd'\n",
    "                              ).autopower() #nyquist_factor=10)#) #method='fast')\n",
    "\n",
    "minimum_frequency = min(frequency)\n",
    "maximum_frequency = max(frequency)\n",
    "print(f'min(frequency) = {min(frequency)}, max(frequency) = {max(frequency)}')\n",
    "print(f'min(power) = {min(power)}, max(power) = {max(power)}')\n",
    "print(f'min(vespa_dropna.IWV_THAAO) = {min(vespa_dropna.IWV_THAAO.values)}')\n",
    "print(f'max(vespa_dropna.IWV_THAAO) = {max(vespa_dropna.IWV_THAAO.values)}')\n",
    "print(f'min(vespa_dropna.timestamp) = {min(vespa_dropna.timestamp.values)}')\n",
    "print(f'max(vespa_dropna.timestamp) = {max(vespa_dropna.timestamp.values)}\\n')\n",
    "print(f'len(vespa_dropna) = {len(vespa_dropna)}')\n",
    "print(f'len(power) = {len(power)}')\n",
    "print(f'len(frequency) = {len(frequency)}\\n')\n",
    "\n",
    "#maximum_frequency = 2.0e-6\n",
    "frequency, power = LombScargle(vespa_dropna['timestamp'], \n",
    "                               vespa_dropna['IWV_THAAO'], \n",
    "                               dy#,\n",
    "                               #fit_mean=True#,\n",
    "                               # normalization='psd'\n",
    "                              ).autopower(minimum_frequency=minimum_frequency,\n",
    "                                          maximum_frequency=maximum_frequency,\n",
    "                                          samples_per_peak=10\n",
    "                                         ) \n",
    "\n",
    "#minimum_frequency = min(frequency)\n",
    "print(f'min(frequency) = {min(frequency)}, max(frequency) = {max(frequency)}')\n",
    "print(f'min(power) = {min(power)}, max(power) = {max(power)}')\n",
    "print(f'min(vespa_dropna.IWV_THAAO) = {min(vespa_dropna.IWV_THAAO.values)}')\n",
    "print(f'max(vespa_dropna.IWV_THAAO) = {max(vespa_dropna.IWV_THAAO.values)}')\n",
    "print(f'min(vespa_dropna.timestamp) = {min(vespa_dropna.timestamp.values)}')\n",
    "print(f'max(vespa_dropna.timestamp) = {max(vespa_dropna.timestamp.values)}\\n')\n",
    "print(f'len(vespa_dropna) = {len(vespa_dropna)}')\n",
    "print(f'len(power) = {len(power)}')\n",
    "print(f'len(frequency) = {len(frequency)}\\n')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "ba061322-4485-47b9-bfa2-c1b319e8ddff",
   "metadata": {
    "tags": []
   },
   "source": [
    "frequency, power = LombScargle(vespa_data['timestamp'], \n",
    "                               vespa_data['IWV_THAAO'], \n",
    "                               dy#, \n",
    "                               # normalization='psd'\n",
    "                              ).autopower(minimum_frequency=minimum_frequency,\n",
    "                                          maximum_frequency=1.e-6,\n",
    "                                          samples_per_peak=100)\n",
    "print(f'min(frequency) = {min(frequency)}, max(frequency) = {max(frequency)}')\n",
    "print(f'len(frequency) = {len(frequency)}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d8c55b3-2008-454e-910e-613f26a4e912",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PSD_LS = pd.DataFrame({'f_LS':frequency, 'f_LS_d':frequency*86400, 'PSD_LS':power}) \n",
    "PSD_LS.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0a7f57c-5101-4d8d-be43-22f1d7afa0a5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# print di controllo\n",
    "print(f'len(vespa_dropna) = {len(vespa_dropna)}')\n",
    "print(f'len(vespa_time_range) = {len(vespa_time_range)}')\n",
    "print(f'len(frequency) = {len(frequency)}, len(power) = {len(power)}')\n",
    "\n",
    "#print(f'len(PSD_LS) = {len(PSD_LS['PSD_LS'])}')#, len(f_LS) = {len(PSD_LSf_LS)}\\n')\n",
    "#print(f'dt = {dt} [s], fs = {fs} [s^-1]')\n",
    "#print(f'dt_o = {dt/3600} [hours]')\n",
    "#print(f'T   = {T} [s], df = {df} 1/(N*dt))[Hz]')\n",
    "#print(f'T_a = {T/(365*86400)} [years]\\n')\n",
    "#print(f'min(f_LS) = {min(f_LS)}, max(PSD_LS) = {max(PSD_LS)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c57adbd3-2496-4914-a4ed-939fe3a6e4d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "## drop first row to plot log(x)\n",
    "#LS_PSD_f = LS_PSD_f.iloc[1:].copy()\n",
    "print(len(PSD_LS['f_LS']))\n",
    "PSD_LS['t_LS_d'] = (1./PSD_LS.f_LS_d)  # unità di misura: xt_d [days] \n",
    "PSD_LS.set_index(['f_LS','f_LS_d'],drop=True,inplace=True)\n",
    "PSD_LS.head(2)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "c63c460c-7528-4768-9dd5-8ad97e7fef44",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "source": [
    "from scipy import signal\n",
    "fs = 1/dt \n",
    "T = len(vespa_time_range)*dt  # T=N*dt [s]\n",
    "df = 1/T  # df=1/T=fs/N , dim bin in freq. df = (1/(N*dt))[Hz]=fatt. di conv"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9efe0a86-6e32-4a0d-acf6-e9863e9bdb5e",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "source": [
    "f_ps, PS_f = signal.periodogram(vespa_time_range.PCHIP,fs, \n",
    "                                window='blackman', #'boxcar', #'hann', \n",
    "                                scaling='spectrum')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "3c057486-5860-4243-9204-25fb0038cade",
   "metadata": {
    "tags": []
   },
   "source": [
    "f_psd, PSD_f = signal.periodogram(vespa_time_range.PCHIP,fs, scaling='density')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "183972c5-2c9d-4a34-bcb5-9d241c2aec81",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# tooltips cambia il nome nel hovertool\n",
    "hover = HoverTool(tooltips=[('f_LS','@f_LS'),('t_LS_d','@t_LS_d'),('PSD_LS','@PSD_LS')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2db9a7d-3595-4d0a-bf1c-f8d74eba0ada",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "PSD_LS.hvplot.line(x='f_LS', y='PSD_LS', logx=True, logy=True, grid=True, \n",
    "                   xlabel='frequency [Hz]', ylabel='PSD LS [mm**2]', \n",
    "                   hover_cols=['f_LS','PSD_LS','t_LS_d']).opts(tools=[hover])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "0fbe7c9b-b9b2-4a71-9d12-79ca56d11d50",
   "metadata": {
    "tags": []
   },
   "source": [
    "print(len(frequency), frequency.min(), frequency.max())\n",
    "plt.semilogy(frequency, power)\n",
    "plt.loglog(frequency, power)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "2e4225c5-d845-44ed-96ee-bd3a8c48351e",
   "metadata": {
    "tags": []
   },
   "source": [
    "from astropy.timeseries import LombScargle\n",
    "import astropy.units as u\n",
    "\n",
    "dy = 0.05\n",
    "\n",
    "frequency, power = LombScargle(vespa_time_range['timestamp'], vespa_time_range['IWV_THAAO'], \n",
    "                               dy#, \n",
    "                               # normalization='psd'\n",
    "                              ).autopower()\n",
    "print(len(frequency), frequency.min(), frequency.max())\n",
    "plt.semilogy(frequency, power)\n",
    "plt.loglog(frequency, power)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "637a500b-8810-4e2e-9ba8-b0616eabd0be",
   "metadata": {},
   "source": [
    "# **Prove**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "",
   "language": "python",
   "name": ""
  },
  "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.8.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}