Thule-mix-LAADS_v06.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "tested-function",
   "metadata": {},
   "source": [
    "# __MODIS vs THAAO comparison of Integrated Water Vapor (IWV)__\n",
    "[MODIS](https://modis.gsfc.nasa.gov/about/specifications.php): Moderate Resolution Imaging Spectroradiometer  <br>\n",
    "[THAAO](https://www.thuleatmos-it.it/index.php): Thule High Arctic Atmospheric Observatory, Greenland  (coordinates  __76.5145, -68.7432__, altitude a.s.l. __220 m__)  <br>\n",
    "\n",
    "__MXD05_L2 IWV product__ <br>\n",
    "_References_ <br>\n",
    "[MYD05_L2](https://atmosphere-imager.gsfc.nasa.gov/products/water-vapor) <br>\n",
    "\n",
    "_Algorithm Theoretical Basis Documents (ATBDs)_ <br>\n",
    "[ATBD](https://atmosphere-imager.gsfc.nasa.gov/documentation/atbds-plans-guides) <br>\n",
    "[The MODIS Near-IR Water Vapor Algorithm](https://atmosphere-imager.gsfc.nasa.gov/sites/default/files/ModAtmo/atbd_mod03_0.pdf) <br>\n",
    "[MODIS Atmospheric Profile Retrieval\n",
    "ATBD MOD07 Collection 6](https://modis-images.gsfc.nasa.gov/_docs/MOD07_atbd_v7_April2011.pdf) <br>\n",
    "[Atmospheric Correction Algorithm: Spectral\n",
    "Reflectances MOD08](https://modis.gsfc.nasa.gov/data/atbd/atbd_mod08.pdf) <br>\n",
    "[MODIS Atmosphere L3 Gridded Product\n",
    "(ATBD) & Users Guide](https://atmosphere-imager.gsfc.nasa.gov/sites/default/files/ModAtmo/L3_ATBD_C6_C61_2019_02_20.pdf) <br>\n",
    "\n",
    "__THAAO IWV product__ <br>\n",
    "_References_ <br>\n",
    "[The water Vapor Emission Spectrometer for Polar Atmospheres at 22 GHz](https://www.thuleatmos-it.it/data/vespa22/index.php)<br> \n",
    "[Integrated Water Vapor](https://www.thuleatmos-it.it/dataaccess/vapor/index.php)<br>\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "contrary-spread",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.pandas\n",
    "import holoviews as hv\n",
    "import numpy as np         \n",
    "import pandas as pd\n",
    "import panel as pn\n",
    "import datetime\n",
    "from holoviews.plotting.links import RangeToolLink\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "pd.options.plotting.backend = 'holoviews'\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fuzzy-automation",
   "metadata": {},
   "source": [
    "#### __MXD05_L2 and THAAO IWV products__ \n",
    "\n",
    "MODIS data were extracted as timeseries of IWV Near-IR data on a Point geometry centered on the THAAO Observatory from the MXD05_L2 products uploaded in the ADAM repository through the ADAM GUI. Within the framework of the RELIANCE project the repository has been populated of MODIS MXD05_L2 data covering the Greenland region upon our request. Other MODIS products have been requested and will be uploaded on the ADAM repository to explore other satellite products, including MODIS LST and AOD, and possibly validate them by close comparison with similar ground-based measurements routinely carried out at THAAO Observatory. The choice of the ADAM framework is motivated by possibility to upload several years of selected satelllite products on the repository, to easily extract temporal and regional subsets in the form of timeseries both through the ADAM GUI and the ADAM API developed within the RELIANCE project. The latter possibility in particular allows more flexibility to fully exploit the repository content and process the data by means of Jupyter notebooks and Python language, while the GUI is designed for an easier usage and allows straightforwrd dataset exploration and simpler data subset extraction. <br>\n",
    "\n",
    "Thule (Greenland) IWV data timeseries have been downloaded from the [Integrated Water Vapor](https://www.thuleatmos-it.it/dataaccess/vapor/index.php) THAAO data access facility were the data are freely accessible. Likewise MODIS satellite IWV product, these data contain vertical column integrated amount  water vapor or precipitable water retrieved by means of the water Vapor Emission Spectrometer for Polar Atmospheres at 22 GHz (VESPA-22).  <br> \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "flexible-assessment",
   "metadata": {},
   "source": [
    "#### __Loading MODIS and Thule timeseries__  \n",
    "IWV timeseries of MODIS in .csv and THAAO in .txt formats are read into pandas dataframes. <br>\n",
    "MODIS data have an additional 13% error column derived from the maximum error claimed in the NearIR IWV product ATBD (see sect. 3.1.4, page 16). <br>\n",
    "THAAO IWV data are given in $kg/m^2$ units and thus divided for factor of 10 to have homogeoeous units with MODIS IWV data that are given in $cm$. <br>  "
   ]
  },
  {
   "cell_type": "raw",
   "id": "16d9fb4e-50fd-4b1c-9baf-00288ef70a69",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "source": [
    "modis=pd.read_csv('mxd05_l2_wvni_thule.csv',header=1,\n",
    "                  names=[\"time\", \"IWV_MODIS\"],sep=';', #index_col=0,\n",
    "                  parse_dates=['time'],\n",
    "                  date_parser=lambda col: pd.to_datetime(col, utc=True))\n",
    "modis['time']=modis['time'].dt.tz_localize(None)\n",
    "modis['Err_IWV_MODIS']=modis[\"IWV_MODIS\"]*0.13\n",
    "modis=modis.set_index('time',drop=False)\n",
    "modis.index=modis.index.rename('Time')\n",
    "#modis.info()\n",
    "modis_adam = modis\n",
    "modis_adam.tail()\n",
    "#modis.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "answering-universal",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def vespa_read(fname):\n",
    "    thaao = pd.read_csv(fname, sep=\"\\t\", header=5, \n",
    "                      names=[\"time\",\"IWV_THAAO\"],\n",
    "                      parse_dates=['time'],\n",
    "                      date_parser=lambda col: pd.to_datetime(col))\n",
    "    thaao['IWV_THAAO']=thaao['IWV_THAAO']/10.\n",
    "    thaao['err_IWV']=thaao['IWV_THAAO']*0.05\n",
    "    thaao=thaao.set_index('time', drop=False)\n",
    "    thaao.index=thaao.index.rename('Time')\n",
    "    thaao['time2']=thaao['time']\n",
    "    thaao['key']=thaao.index\n",
    "    return thaao"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7febf48d-5249-41df-8872-5e22485b4116",
   "metadata": {
    "tags": []
   },
   "source": [
    "# MODIS LAADS IWV NIR\n",
    "\n",
    "path_laads = 'MYD05_5km.csv'\n",
    "modis=pd.read_csv(path_laads,header=0,sep='\\t',# index_col=0,\n",
    "                  #names=[\"time\", \"IWV_MODIS\"]\n",
    "                  parse_dates=['time'],\n",
    "                  date_parser=lambda col: pd.to_datetime(col, utc=True))\n",
    "modis['time']=modis['time'].dt.tz_localize(None)\n",
    "modis['Err_IWV_MODIS']=modis[\"IWV_MODIS\"]*0.13\n",
    "modis=modis.set_index('time',drop=False)\n",
    "modis.index=modis.index.rename('Time')\n",
    "#modis.info()\n",
    "modis.head()\n",
    "#modis.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f371385-18a7-4d88-8857-13544b36330e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def modis_read(fname):\n",
    "    path_laads = fname\n",
    "    modis = pd.read_csv(path_laads,header=0,sep='\\t',\n",
    "                      parse_dates=['time'],\n",
    "                      date_parser=lambda col: pd.to_datetime(col, utc=True))\n",
    "    modis['time'] = modis['time'].dt.tz_localize(None)\n",
    "    modis['Err_IWV_MODIS'] = modis[\"IWV_MODIS\"]*0.13\n",
    "    modis = modis.set_index('time',drop=False)\n",
    "    modis.index=modis.index.rename('Time')\n",
    "    modis.head\n",
    "    return modis"
   ]
  },
  {
   "cell_type": "raw",
   "id": "04a8c118-15fd-450c-b007-4db4e04942b8",
   "metadata": {},
   "source": [
    "wdir_hatpro = 'HATPRO/'\n",
    "DATA_dir = '/home/luca/WORK/DATA/Groundbased/THULE/'\n",
    "ELAB_dir = '/home/luca/WORK/ELAB/Groundbased/THULE/'\n",
    "\n",
    "\n",
    "def read_hatpro_year(fname):\n",
    "    year = fname[-15:-11]\n",
    "    return year\n",
    "import glob\n",
    "from pathlib import *\n",
    "\n",
    "p = Path('.')\n",
    "butta = list(p.glob('**/*.ipynb'))\n",
    "\n",
    "p_data = PurePath(p.home(),'WORK','DATA','Groundbased','THULE')\n",
    "p_elab = PurePath(p.home(),'WORK','ELAB','Groundbased','THULE')\n",
    "#year = [read_hatpro_year(y) for y in fnlist]\n",
    "butta\n",
    "#p_elab.mkdir(mode=0o777, parents=False)\n",
    "dir = p_elab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8bfdbe6-b40d-4b70-948c-d24daac59078",
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_IO_dirs(path_segments):\n",
    "    \"\"\" Final path_segments as a list, to be added to\n",
    "    the baseline directory fo a complete absolute path\n",
    "    to DATA and ELAB dirs.\n",
    "    \"\"\"\n",
    "    from pathlib import Path\n",
    "    \n",
    "    # writing the path segments\n",
    "    homedir = Path.home()          # similar to ~/\n",
    "    data_dir_segments = [homedir,'WORK','DATA','Groundbased','THULE']\n",
    "    [data_dir_segments.append(s) for s in path_segments]\n",
    "    #elab_dir_segments = data_dir_segments\n",
    "    #elab_dir_segments = data_dir_segments[:]\n",
    "    import copy\n",
    "    elab_dir_segments = copy.deepcopy(data_dir_segments)\n",
    "    \n",
    "    elab_dir_segments[2] = 'ELAB'\n",
    "    return data_dir_segments, elab_dir_segments\n",
    "\n",
    "def check_dirs_existence_mkdir(p_data, p_elab):\n",
    "    \"\"\" Check out DATA and ELAB I/O dirs existence and \n",
    "    eventually creates them.\n",
    "    \"\"\"\n",
    "    p_data.mkdir(parents=True, exist_ok=True)\n",
    "    p_elab.mkdir(parents=True, exist_ok=True)\n",
    "    return\n",
    "\n",
    "def field_campaign_name():\n",
    "    pass\n",
    "    return fc_name\n",
    "\n",
    "def sensor_name():\n",
    "    pass\n",
    "    return sensor_name\n",
    "\n",
    "\n",
    "\n",
    "from pathlib import *\n",
    "\n",
    "#transforming path segments into concrete paths\n",
    "path_segments  = ['HATPRO']\n",
    "path_data, path_elab = set_IO_dirs(path_segments)\n",
    "\n",
    "print(f'data: {path_data}    id: {id(path_data)}')\n",
    "print(f'elab: {path_elab}    id: {id(path_elab)}')\n",
    "\n",
    "#p_data, p_elab = Path(*final_path_segments)\n",
    "p_data = Path(*path_data)\n",
    "p_elab = Path(*path_elab)\n",
    "\n",
    "print(f'data: {p_data}    id: {id(p_data)}')\n",
    "print(f'elab: {p_elab}    id: {id(p_elab)}')\n",
    "\n",
    "check_dirs_existence_mkdir(p_data, p_elab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea11493-a0c4-4de6-891c-de9a7230f831",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hatpro_read(fname, year='2016'):\n",
    "    #year = read_hatpro_year(fname)\n",
    "    hatpro = pd.read_csv(fname, delim_whitespace=True, header=0)\n",
    "    hatpro['Time'] = hatpro.JD_2016.to_frame().applymap(lambda x: \n",
    "                                                        datetime.datetime(\n",
    "                                                            2016,1,1,0,0,0)+\n",
    "                                                        datetime.timedelta(days=x-1.))\n",
    "    hatpro.Time = hatpro.Time.dt.tz_localize(None)\n",
    "    hatpro.set_index('Time', drop=True,inplace=True)\n",
    "    hatpro.drop(columns=['JD_rif', 'JD_2016', 'mcorIWV', 'stdcoIWV', 'biasIWV',\n",
    "                         'n_med','RF_med'], inplace=True)\n",
    "    return hatpro\n",
    "\n",
    "\n",
    "HATPRO = hatpro_read('/home/luca/WORK/DATA/Groundbased/THULE/HATPRO/IWV_2016_HATPRO.txt')\n",
    "\n",
    "HATPRO"
   ]
  },
  {
   "cell_type": "raw",
   "id": "c2f6a69f-7522-48cc-9ff0-d09061b6d904",
   "metadata": {},
   "source": [
    "def gnss_read(filename, year):\n",
    "    gnss_data = pd.read_csv(filename,header=[0,1],delim_whitespace=True,\n",
    "                             date_parser=lambda col: pd.to_datetime(col,\n",
    "                                                                    utc=True))\n",
    "    gnss_data['Time'] = gnss_data.DOY.applymap(lambda x: datetime.datetime(\n",
    "        int(year),1,1,0,0,0)+datetime.timedelta(days=x-1.))\n",
    "    gnss_data.Time = gnss_data.Time.dt.tz_localize(None)\n",
    "    gnss_data = gnss_data.set_index('Time', drop=True)\n",
    "    return gnss_data\n",
    "\n",
    "\n",
    "year=[2016,2017,2018,2019]\n",
    "gnss_data = pd.DataFrame()\n",
    "for year in year:\n",
    "    gnss_data = gnss_data.append(gnss_read(f'GNSS/THU3_{year}global.txt',year))\n",
    "gnss_data.head() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18460672-3480-4471-aec8-7de0e14c4017",
   "metadata": {},
   "outputs": [],
   "source": [
    "HATPRO = hatpro_read('/home/luca/WORK/DATA/Groundbased/THULE/HATPRO/IWV_2017_HATPRO.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c354b5c0-9923-4a94-bfa2-03e054d5c6be",
   "metadata": {},
   "outputs": [],
   "source": [
    "HATPRO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa0b8217-f509-4681-a6f2-1788eb5baa7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "VESPA = vespa_read('/home/luca/WORK/DATA/Groundbased/THULE/VESPA22/Vapor_20160716_20220311.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fd945c3-fa35-4cff-ac75-4d2f3ba0e8ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "MOD_IR = modis_read('/home/luca/WORK/DATA/MODIS/THULE/Thule_LAADS_IWV_IR/MOD05_5km_totale.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d76fb7a-949b-4a6d-8726-fbb0ab5a9555",
   "metadata": {},
   "outputs": [],
   "source": [
    "MYD_IR = modis_read('/home/luca/WORK/DATA/MODIS/THULE/Thule_LAADS_IWV_IR/MYD05_5km_totale.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61144d91-8d54-489f-b813-4e5eb5ae21e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "VESPA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8841b37a-f180-4811-b51d-1edb7626d8c8",
   "metadata": {},
   "source": [
    "### Drop rows with values >5 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd061e5a-cade-405f-a2f6-c2a1b6c9e17e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_drop(df, condition):\n",
    "    df = df.drop(df[condition].index)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bea5f082-0fbf-41db-a7c5-3ffc96a02e91",
   "metadata": {},
   "outputs": [],
   "source": [
    "MOD_IR = df_drop(MOD_IR, MOD_IR.IWV_MODIS > 5. )\n",
    "MOD_IR = df_drop(MOD_IR, MOD_IR.longitude < -65000. )\n",
    "\n",
    "MYD_IR = df_drop(MYD_IR, MYD_IR.IWV_MODIS > 5. )\n",
    "MYD_IR = df_drop(MYD_IR, MYD_IR.longitude < -65000. )\n",
    "number_MXD_data_points = len(MOD_IR.IWV_MODIS) + len(MYD_IR.IWV_MODIS)\n",
    "MYD_IR.head(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a46c1640-3b8b-47a7-be48-32c30bcb38ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_MXD_data_points = len(MOD_IR.IWV_MODIS) + len(MYD_IR.IWV_MODIS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "attempted-florence",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### **MODIS and THAAO Integrated Water Vapor (IWV) timeseries**\n",
    "##### Simple linear plots of the whole MODIS and THAAO IVW timeseries."
   ]
  },
  {
   "cell_type": "raw",
   "id": "fd8f52cb-45a3-4c73-946a-5a6934bfd763",
   "metadata": {
    "tags": []
   },
   "source": [
    "(MOD_IR.hvplot(yaxis=None, height=350, fields={'date': 'Date'}, color=[\"#718dbf\",'red'],title ='MOD') \n",
    "+ MYD_IR.hvplot(yaxis=None, height=350, fields={'date': 'Date'}, color=[\"#718dbf\",'red'],title='MYD'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fc7aa22-4d0e-4df4-8ee9-48475095ac35",
   "metadata": {},
   "outputs": [],
   "source": [
    "MOD_IR['Err_IWV_MODIS']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "impressed-aquatic",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "width = 1100\n",
    "height = 300 \n",
    "size = 3\n",
    "scale = 0.3\n",
    "color_MOD = 'cyan'\n",
    "color_MYD = '#F2583E'\n",
    "color_THAAO = 'green'\n",
    "\n",
    "(MOD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', ylabel='IWV  (cm)', \n",
    "                                 height=height, width=width, size=size, \n",
    "                                 color=color_MOD, title ='IWV MOD')\n",
    "* MOD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "+ MYD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', ylabel='IWV  (cm)', \n",
    "                                  height=height, width=width, size=size, \n",
    "                                  color=color_MYD,title ='IWV MYD')\n",
    "* MYD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "+ (MOD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', ylabel='IWV  (cm)', \n",
    "                                   height=height, width=width, size=size, \n",
    "                                   color=color_MOD,title ='IWV MOD')\n",
    "* MOD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "* MYD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', ylabel='IWV  (cm)', \n",
    "                                  height=height, width=width, size=size, \n",
    "                                  color=color_MYD,title ='IWV MYD')\n",
    "* MYD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "  ).opts(title='MOD+MYD')).cols(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17bb0e93-6873-499d-98fa-9098605c375b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#MOD_IR #['Err_IWV_MODIS']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "light-turkish",
   "metadata": {},
   "outputs": [],
   "source": [
    "thaao_timeseries = VESPA.hvplot.scatter(yaxis=True, \n",
    "                                        height=height, width=width, \n",
    "                                        size=size, \n",
    "                                        color=color_THAAO, \n",
    "                                        fields={'date': 'Date'},\n",
    "                                        title='IWV THAAO - VESPA', label='VESPA',\n",
    "                                        legend='top_right')\n",
    "thaao_timeseries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "geological-benefit",
   "metadata": {},
   "source": [
    "Overlay of the IWV plots on the same time reference <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "danish-civilian",
   "metadata": {},
   "outputs": [],
   "source": [
    "MOD_timeseries = MOD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', \n",
    "                                                 ylabel='IWV (cm)',label='MOD',\n",
    "                                                 height=height, width=width, \n",
    "                                                 size=size, color=color_MOD, \n",
    "                                                 title='IWV MOD - VESPA',\n",
    "                                                 legend='top_right')\n",
    "MYD_timeseries = MYD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', ylabel='IWV (cm)',label='MYD',\n",
    "                                                 height=height, width=width, \n",
    "                                                 size=size, color=color_MYD, \n",
    "                                                 title='IWV MYD - VESPA',\n",
    "                                                 legend='top_right')\n",
    "\n",
    "\n",
    "(MOD_timeseries *thaao_timeseries\n",
    "+ MYD_timeseries * thaao_timeseries\n",
    "+ MOD_timeseries*MYD_timeseries * thaao_timeseries\n",
    "* MOD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "* MYD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)).cols(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e56fba-9733-481d-8ac6-c58e8fcf293c",
   "metadata": {},
   "outputs": [],
   "source": [
    "title = f'IWV MXD-VESPA                                         # data points={number_MXD_data_points}'\n",
    "(MOD_timeseries*MYD_timeseries * thaao_timeseries).opts(title=title, show_grid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "506f14f5-1133-471a-aa26-46fccc74625f",
   "metadata": {},
   "source": [
    "### __TODO filtrare i dati MODIS con il QA__ !!!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "treated-disease",
   "metadata": {},
   "source": [
    "#### __Merging pandas dataframes to filter MODIS and THAAO data points closest in time__\n",
    "\n",
    "Due to time lenght of the timeseries, even if the notebook representation allows to zoom in into the datasets, it is still not easy to see how close in time two data points are. Since it makes little sense to compare two MODIS and THAAO IWV data points if they are too far apart, and they can easily be hours or days far apart, therefore easily mistaken if one relies solely on the visual inspection of the plots, a tool to select and compare only the data points which are closer in time is required. The following code allows to filter the data on a timespan basis by suitably merging modis and thaao pandas dataframes while defining the Timedelta parameter (expressed in s,m,d, i.e. seconds, minutes or days respectively) that allows to operate the time-based filter. Different choices of Timedelta allow to widen or narrow the time span around each MODIS data point and select for the comparison only the THAAO data points that fall within the chosen interval. This selection is based on MODIS data points because THAAO measures are far more dense. <br>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8325912-e6f5-486d-ba48-6dcacc5d8948",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def closest_in_time(df1,df2,delta_time='30m'):\n",
    "    \"\"\"delta_time : string\n",
    "    add RSME to new_df\n",
    "    \"\"\"\n",
    "    newdf=pd.merge_asof(df1,df2,on='time',tolerance=pd.Timedelta(delta_time),direction='nearest')\n",
    "    newdf['diff']=(newdf.time-newdf.time2).abs()\n",
    "    newdf['diff_int']=newdf['diff'].astype(int)/1000000000\n",
    "    newdf['diff_IWV']=(newdf.IWV_MODIS-newdf.IWV_THAAO)\n",
    "    newdf['diff_perc_IWV']=(newdf.IWV_MODIS/newdf.IWV_THAAO)*100. - 100.\n",
    "    newdf['diff_distance']=np.sqrt((newdf.latitude-76.5145)**2 + (newdf.longitude-68.7432)**2)\n",
    "    newdf=newdf.sort_values('diff').drop_duplicates('key')\n",
    "    newdf=newdf.set_index(['time','time2'], drop=False).sort_index()\n",
    "    newdf.index=newdf.index.rename(['Time','Time2'])\n",
    "    newdf=newdf[1:]\n",
    "    newdf['RMSE_IWV']=((newdf.IWV_MODIS-newdf.IWV_THAAO)**2).mean()**.5\n",
    "    newdf.info()\n",
    "    newdf.head()\n",
    "    return newdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df5edb01-92b1-489c-9cfb-d7602d0dca2a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "deltatime = '10m'\n",
    "newdf_MOD = closest_in_time(MOD_IR,VESPA, delta_time=deltatime)\n",
    "newdf_MYD = closest_in_time(MYD_IR,VESPA, delta_time=deltatime)\n",
    "number_MXD_data_points_deltatime = len(newdf_MOD.IWV_MODIS) + len(newdf_MYD.IWV_MODIS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd1ab70e-19d7-4293-9527-53990a2d0db2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "newdf_MOD.head()\n",
    "newdf_MOD.vza.max()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "cb6b712d-234e-4977-9033-990fbd090ee7",
   "metadata": {
    "tags": []
   },
   "source": [
    "newdf_MYD.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37ae7e04-348e-484e-a2b9-02d46569f51e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "newdf_TOT = pd.concat([newdf_MOD,newdf_MYD], sort=True)\n",
    "#newdf_TOT.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ac7de8b-2374-4181-96be-ed9b3643d8cc",
   "metadata": {},
   "source": [
    "### Numero punti per mese "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "464ea274-d9d5-4863-9c51-26381832eedd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "newdf_MOD.time.groupby([newdf_MOD['time'].dt.month.rename('month')]).agg({'count'}).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14f8f89e-db35-471a-a9df-0d3ac7eca2c0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "newdf_MYD.time.groupby([newdf_MYD['time'].dt.month.rename('month')]).agg({'count'}).transpose()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "49051c07-ee51-428f-a0e2-2a1173fbe570",
   "metadata": {
    "tags": []
   },
   "source": [
    "newdf['RMSE_IWV']=((newdf.IWV_MODIS-newdf.IWV_THAAO)**2).mean()**.5\n",
    "newdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d05c4c22-349c-4d26-bfdd-7d033f9063d5",
   "metadata": {},
   "source": [
    "### Timeseries su diff temporali "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mysterious-keyboard",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "newdf_MOD_timeseries = newdf_MOD.IWV_MODIS.hvplot.scatter(x='Time', \n",
    "                                                          ylabel='IWV (cm)',\n",
    "                                                          label='MOD',\n",
    "                                                          height=height, \n",
    "                                                          width=width, \n",
    "                                                          size=size, \n",
    "                                                          color=color_MOD, \n",
    "                                                          title= 'IWV MOD-VESPA',\n",
    "                                                          legend='top_right')\n",
    "\n",
    "newdf_MYD_timeseries = newdf_MYD.IWV_MODIS.hvplot.scatter(x='Time', ylabel='IWV (cm)',\n",
    "                                                          label='IWV_MYD',\n",
    "                                                          height=height, \n",
    "                                                          width=width, \n",
    "                                                          size=size, \n",
    "                                                          color=color_MYD, \n",
    "                                                          title= 'IWV MYD-VESPA',\n",
    "                                                          legend='top_right')\n",
    "\n",
    "(newdf_MOD_timeseries*newdf_MOD.IWV_THAAO.hvplot.scatter(x='Time2',\n",
    "                                                         height=height, \n",
    "                                                         width=width, \n",
    "                                                         size=size, \n",
    "                                                         color=color_THAAO)\n",
    "+ newdf_MYD_timeseries*newdf_MYD.IWV_THAAO.hvplot.scatter(x='Time2',\n",
    "                                                          height=height, \n",
    "                                                          width=width, \n",
    "                                                          size=size, \n",
    "                                                          color=color_THAAO)\n",
    "+ (newdf_MOD_timeseries *newdf_MYD_timeseries\n",
    "* MOD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "* MYD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)   \n",
    "* newdf_MOD.IWV_THAAO.hvplot.scatter(x='Time2', height=height, \n",
    "                                     width=width, size=size, \n",
    "                                     color=color_THAAO)\n",
    "* newdf_MYD.IWV_THAAO.hvplot.scatter(x='Time2',height=height, \n",
    "                                     width=width, size=size, \n",
    "                                     color=color_THAAO)).opts(title='IWV MXD-VESPA')).cols(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecff34d4-bec6-45b0-ac28-9bcb7b8e911f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ok\n",
    "title = f'IWV MXD-VESPA              deltatime={deltatime}      # data points={number_MXD_data_points_deltatime}'\n",
    "(newdf_MOD_timeseries *newdf_MYD_timeseries \n",
    "* newdf_MOD.IWV_THAAO.hvplot.scatter(x='Time2', height=height, \n",
    "                                     width=width, size=size, \n",
    "                                     color=color_THAAO)\n",
    "* newdf_MYD.IWV_THAAO.hvplot.scatter(x='Time2',height=height, \n",
    "                                     width=width, size=size, \n",
    "                                     color=color_THAAO)).opts(title=title, show_grid=True)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b2ee3d72-cb20-4768-8ad1-c56c73398343",
   "metadata": {
    "tags": []
   },
   "source": [
    "newdf_MOD.hvplot.table(['Time','Time2','Err_IWV_MODIS','IWV_MODIS','IWV_THAAO','diff','diff_IWV', 'vza', 'vaa'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92a4599a-8dc2-4180-b391-2289fd844a0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "(MOD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', ylabel='IWV  (cm)', \n",
    "                                   height=height, width=width, size=size, \n",
    "                                   color=color_MOD,title ='IWV MOD')\n",
    "* MOD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "* MYD_IR.IWV_MODIS.hvplot.scatter(xlabel='Time', ylabel='IWV  (cm)', \n",
    "                                  height=height, width=width, size=size, \n",
    "                                  color=color_MYD,title ='IWV MYD')\n",
    "* MYD_IR.hvplot.errorbars(y='IWV_MODIS', yerr1='Err_IWV_MODIS', hover=False)\n",
    "  ).opts(title='MOD+MYD') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "659980cf-baac-44d8-b7cf-f994faf20ba1",
   "metadata": {},
   "source": [
    "### bisettrice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf036732-2c84-489f-8be9-4de5b374f36d",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = hv.Slope(1, 0).opts(line_color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb2e3cdd-2d86-4cde-a3de-aa6f9dfb5049",
   "metadata": {},
   "source": [
    "### Regressione con Holoviews"
   ]
  },
  {
   "cell_type": "raw",
   "id": "04731b36-91b5-4882-8259-1db7a5511418",
   "metadata": {},
   "source": [
    "def regression_holoviews(x,y,color='orange', overlay_values=True):\n",
    "    \"\"\"\n",
    "    Input: x values and y values\n",
    "    overlay values: overlay of coefficients on the plot\n",
    "\n",
    "    Return a list having:\n",
    "    \n",
    "    1- all regress element\n",
    "    2- slope value,\n",
    "    3- intercept value \"\"\"\n",
    "\n",
    "    coords = {'x':x,'y':y}\n",
    "    xxxx= hv.Scatter(coords)\n",
    "    regress=hv.Slope.from_scatter(xxxx).opts(color=color)\n",
    "    \n",
    "    if overlay_values == True:      \n",
    "        return [regress*hv.Text(max(x)-0.3,min(y)+0.1, \n",
    "                                f'y = {regress.slope} x + {regress.y_intercept}'),\n",
    "                regress.slope,\n",
    "                regress.y_intercept]\n",
    "    else:\n",
    "        return [regress,\n",
    "                regress.slope,\n",
    "                regress.y_intercept]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19f9d312-02b8-41b7-a64f-6118c683842a",
   "metadata": {},
   "source": [
    "### **Regressione con sklearn**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "022d2c21-f23e-4141-b086-603a52ca11ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def regression_sklearn(x,y, color='orange', overlay_values=True):\n",
    "    \"\"\"\n",
    "    Input: Pd.dataframe column\n",
    "    overlay values: overlay of coefficients on the plot\n",
    "       \n",
    "    Return a list having:\n",
    "    1- all regress element\n",
    "    2- R coefficient\n",
    "    3- slope value,\n",
    "    4- intercept value \"\"\"\n",
    "    \n",
    "    model = LinearRegression()\n",
    "    x = x.values.reshape(-1,1)\n",
    "    y = y.values.reshape(-1,1)\n",
    "    reg = LinearRegression().fit(x,y)\n",
    "    regress = hv.Slope(reg.coef_[0][0],reg.intercept_[0]).opts(color=color)\n",
    "    \n",
    "    if overlay_values == True:\n",
    "        return [regress,#*hv.Text(max(x)-0.3,min(y)+0.1, \n",
    "                        #        f'y = {reg.coef_[0][0]} x + {reg.intercept_[0]}'),\n",
    "                reg.score(x, y),\n",
    "                reg.coef_[0][0],\n",
    "                reg.intercept_[0]]\n",
    "    else:\n",
    "        return [regress,reg.score(x, y),\n",
    "                reg.coef_[0][0],reg.intercept_[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50959868-7497-4b94-9df6-4a7759e30928",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def regression_eq_text(df_x, df_y, regress, x_offset=-0.3,y_offset=+0.1):\n",
    "    return hv.Text(max(df_x)+x_offset, min(df_y)+y_offset, \n",
    "                                f'y = {regress[2]:.3f} x + {regress[3]:.3f} \\n R = {regress[1]:.3f}')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "021a1863-4767-4fc0-892c-47b64ffca2f9",
   "metadata": {
    "tags": []
   },
   "source": [
    "regress_MOD = regression_holoviews(newdf_MOD.IWV_MODIS,newdf_MOD.IWV_THAAO)\n",
    "regress_MYD = regression_holoviews(newdf_MYD.IWV_MODIS,newdf_MYD.IWV_THAAO,color='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "865cc697-613f-4fe4-96c2-28f85f7d72bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "regress_MOD_sk = regression_sklearn(newdf_MOD.IWV_MODIS,newdf_MOD.IWV_THAAO)\n",
    "regress_MOD_sk_text = regression_eq_text(newdf_MOD.IWV_MODIS,\n",
    "                                         newdf_MOD.IWV_THAAO, regress_MOD_sk)\n",
    "regress_MYD_sk = regression_sklearn(newdf_MYD.IWV_MODIS,\n",
    "                                    newdf_MYD.IWV_THAAO, color='black')\n",
    "regress_MYD_sk_text = regression_eq_text(newdf_MYD.IWV_MODIS,\n",
    "                                         newdf_MYD.IWV_THAAO, regress_MYD_sk)\n",
    "regress_TOT_sk = regression_sklearn(newdf_TOT.IWV_MODIS,newdf_TOT.IWV_THAAO, \n",
    "                                    color='cyan')\n",
    "regress_TOT_sk_text = regression_eq_text(newdf_TOT.IWV_MODIS,\n",
    "                                         newdf_TOT.IWV_THAAO, regress_TOT_sk)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cognitive-continuity",
   "metadata": {},
   "source": [
    "### **MODIS vs THAAO Integrated Water Vapor (IWV) timeseries scatter plots**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "845e2d68-5d97-4740-b514-d4d2898fca24",
   "metadata": {},
   "source": [
    "### Scatter su differenza secondi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aaceb79-6d5e-4a70-82ea-3dd706c47e95",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ThaaoVsMOD = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                      height=500, width=600, c='vza',\n",
    "                                      s='diff_int', scale=0.5,\n",
    "                                      xlabel='IWV MOD (cm)', \n",
    "                                      ylabel='IWV THAAO (cm)'\n",
    "                                     ).hist(num_bins=21,dimension=['IWV_MODIS','IWV_THAAO'])\n",
    "\n",
    "\n",
    "ThaaoVsMYD = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO', \n",
    "                                      height=500, width=600, c='vza',\n",
    "                                      s='diff_int', scale=0.5, \n",
    "                                      xlabel='IWV MYD (cm)', \n",
    "                                      ylabel='IWV THAAO (cm)'\n",
    "                                     ).hist(num_bins=21,dimension=['IWV_MODIS','IWV_THAAO'])\n",
    "\n",
    "(ThaaoVsMOD*slope*regress_MOD_sk[0]*regress_MOD_sk_text\n",
    "+ ThaaoVsMYD*slope*regress_MYD_sk[0]\n",
    " *regress_MYD_sk_text\n",
    "+ ThaaoVsMOD*ThaaoVsMYD*slope*regress_TOT_sk_text\n",
    "* regress_MOD_sk[0]*regress_MYD_sk[0]*regress_TOT_sk[0]).cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "837f76d8-2936-4fe9-b2bc-744706f2d4c3",
   "metadata": {},
   "source": [
    "### Scatter su differenza distanza geografica"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc06bdbb-ee4e-4aa1-961c-23c5443ebae6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ThaaoVsMOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=500, width=600, c='vza',\n",
    "                                              s='diff_distance', scale=0.3,\n",
    "                                              xlabel='IWV MOD (cm)', \n",
    "                                              ylabel='IWV THAAO (cm)'\n",
    "                                             ).hist(num_bins=21,dimension=['IWV_MODIS','IWV_THAAO'])\n",
    "\n",
    "ThaaoVsMYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=500, width=600, c='vza',\n",
    "                                              s='diff_distance', scale=0.3,\n",
    "                                              xlabel='IWV MYD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)'\n",
    "                                             ).hist(num_bins=21,dimension=['IWV_MODIS','IWV_THAAO'])\n",
    "\n",
    "(ThaaoVsMOD_scatter*slope*regress_MOD_sk[0]*regress_MOD_sk_text\n",
    "+ ThaaoVsMYD_scatter*slope*regress_MYD_sk[0]*regress_MYD_sk_text\n",
    "+ ThaaoVsMOD_scatter*ThaaoVsMYD_scatter*slope*regress_TOT_sk_text\n",
    " *regress_MOD_sk[0]*regress_MYD_sk[0]*regress_TOT_sk[0]).opts(width=400, height=300).cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51993457-e623-43f8-8c15-3a7ee9658cd4",
   "metadata": {},
   "source": [
    "### Scatter totale per piattaforma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1bd5997-2b32-424e-8f06-b82cecb8980f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ThaaoVsMOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width, \n",
    "                                              c='vza',\n",
    "                                              s='diff_distance', scale=scale,\n",
    "                                              xlabel='IWV MOD  (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)'\n",
    "                                             ).hist(num_bins=21,dimension=['IWV_MODIS','IWV_THAAO'])\n",
    "\n",
    "\n",
    "ThaaoVsMYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width, \n",
    "                                              c='vza',\n",
    "                                              s='diff_distance', scale=scale, \n",
    "                                              xlabel='IWV MYD (cm)', \n",
    "                                              ylabel='IWV THAAO (cm)'\n",
    "                                             ).hist(num_bins=21,dimension=['IWV_MODIS','IWV_THAAO'])\n",
    "\n",
    "(ThaaoVsMOD_scatter*slope*regress_MOD_sk[0]*regress_MOD_sk_text\n",
    "+ ThaaoVsMYD_scatter*slope*regress_MYD_sk[0]*regress_MYD_sk_text\n",
    "+ ThaaoVsMOD_scatter*ThaaoVsMYD_scatter*slope*regress_TOT_sk_text\n",
    " *regress_MOD_sk[0]*regress_MYD_sk[0]*regress_TOT_sk[0]).cols(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c96de5c3-4475-4833-899d-9210d3c6482f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ok\n",
    "height=500\n",
    "width=700\n",
    "title = f'IWV MXD-VESPA              deltatime={deltatime}      # data points={number_MXD_data_points_deltatime}'\n",
    "\n",
    "ThaaoVsMOD_scatter_ok = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width, \n",
    "                                              c='vza',\n",
    "                                              s='diff_distance', scale=scale,\n",
    "                                              xlabel='IWV MXD  (cm)',\n",
    "                                              ylabel='IWV VESPA (cm)')\n",
    "\n",
    "ThaaoVsMYD_scatter_ok = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width, \n",
    "                                              c='vza',\n",
    "                                              s='diff_distance', scale=scale, \n",
    "                                              #xlabel='IWV MYD (cm)', \n",
    "                                              ylabel='IWV VESPA (cm)')\n",
    "\n",
    "(ThaaoVsMOD_scatter_ok*ThaaoVsMYD_scatter_ok\n",
    "*slope\n",
    "*regress_TOT_sk_text\n",
    "#*regress_MOD_sk[0]*regress_MYD_sk[0]\n",
    "*regress_TOT_sk[0]).opts(title=title, show_grid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98c19afa-a1f6-4af6-9cd2-141971171e2b",
   "metadata": {},
   "source": [
    "### Stagionalità"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf55dea3-1006-476d-9d39-1dad8ff9fecd",
   "metadata": {},
   "outputs": [],
   "source": [
    "xlim = [0,2]\n",
    "ylim = [0,2]\n",
    "regress_MOD_s_text = regression_eq_text(newdf_MOD.IWV_MODIS,\n",
    "                                        newdf_MOD.IWV_THAAO, regress_MOD_sk,\n",
    "                                        x_offset=xlim[0]-0.5,y_offset=ylim[0]+0.5)\n",
    "regress_MYD_s_text = regression_eq_text(newdf_MYD.IWV_MODIS,\n",
    "                                         newdf_MYD.IWV_THAAO, regress_MYD_sk,\n",
    "                                        x_offset=xlim[0]-0.5,y_offset=ylim[0]+0.5)\n",
    "regress_TOT_s_text = regression_eq_text(newdf_TOT.IWV_MODIS,\n",
    "                                         newdf_TOT.IWV_THAAO, regress_TOT_sk,\n",
    "                                        x_offset=xlim[0]-0.5,y_offset=ylim[0]+0.5)\n",
    "\n",
    "ThaaoVsMOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO', \n",
    "                                              height=height, width=width, \n",
    "                                              c='vza', \n",
    "                                              xlim=xlim, ylim=ylim, \n",
    "                                              s='diff_distance', scale=0.3, \n",
    "                                              xlabel='IWV MOD (cm)', \n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              groupby=['time.month'])\n",
    "\n",
    "ThaaoVsMYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO', \n",
    "                                              height=height, width=width, c='vza',\n",
    "                                              xlim=xlim, ylim=ylim, \n",
    "                                              s='diff_distance', scale=0.3, \n",
    "                                              xlabel='IWV MYD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              groupby=['time.month'])\n",
    "\n",
    "(ThaaoVsMOD_scatter*slope*regress_MOD_sk[0]*regress_MOD_s_text\n",
    "+ ThaaoVsMYD_scatter*slope*regress_MYD_sk[0]*regress_MYD_s_text\n",
    "+ (ThaaoVsMOD_scatter*ThaaoVsMYD_scatter*slope*regress_TOT_s_text\n",
    "   *regress_MOD_sk[0]*regress_MYD_sk[0]*regress_TOT_sk[0]\n",
    "  ).opts(title='DDDDD', xlabel= 'IWV MXD (cm)')).cols(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "115d8cde-2512-4d26-87f4-0b6d3f327ff3",
   "metadata": {},
   "source": [
    "### GroupBy month"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31585dc2-87f7-4a1e-9923-48b128992cdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "ThaaoVsMYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=300,xlim=[0,2],ylim=[0,2],\n",
    "                                              scale=0.3, xlabel='IWV MYD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              by=['time.month'], groupby=['time.month'])\n",
    "\n",
    "ThaaoVsMOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=300, xlim=[0,2],ylim=[0,2],\n",
    "                                              scale=0.3, xlabel='IWV MOD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              by=['time.month'], groupby=['time.month'])\n",
    "\n",
    "(ThaaoVsMOD_scatter*slope*regress_MOD_sk[0]*regress_MOD_s_text\n",
    "+ ThaaoVsMYD_scatter*slope*regress_MYD_sk[0]*regress_MYD_s_text\n",
    "+ (ThaaoVsMOD_scatter*ThaaoVsMYD_scatter*slope*regress_MOD_sk[0]*regress_MYD_sk[0]\n",
    "*regress_TOT_sk[0]*regress_TOT_s_text\n",
    "  ).opts(title='DDDDD', xlabel = 'IWV MXD (cm)')).cols(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3af6b39c-2140-4fa4-b1f6-17ff24197e26",
   "metadata": {},
   "outputs": [],
   "source": [
    "ThaaoVsMYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=300,xlim=[0,2],ylim=[0,2],\n",
    "                                              scale=0.3, xlabel='IWV MYD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              by=['time.month'], groupby=['time.month'])\n",
    "\n",
    "ThaaoVsMOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=300, xlim=[0,2],ylim=[0,2],\n",
    "                                              scale=0.3, xlabel='IWV MOD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              by=['time.month'], groupby=['time.month'])\n",
    "\n",
    "pn.panel((ThaaoVsMOD_scatter*slope*regress_MOD_sk[0]*regress_MOD_s_text\n",
    "+ ThaaoVsMYD_scatter*slope*regress_MYD_sk[0]*regress_MYD_s_text\n",
    "+ (ThaaoVsMOD_scatter*ThaaoVsMYD_scatter*slope*regress_MOD_sk[0]*regress_MYD_sk[0]\n",
    "*regress_TOT_sk[0]*regress_TOT_s_text\n",
    "  ).opts(title='DDDDD', xlabel = 'IWV MXD (cm)')).cols(2),\n",
    "         widget_type='scrubber', widget_location='top')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6f4edce-0221-4458-b63f-70e1a5d18279",
   "metadata": {},
   "source": [
    "### By month"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ebb92c4-6994-44fc-a3de-7524a2f7c5bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#cmap = 'hsv'\n",
    "ThaaoVsMOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width,\n",
    "#                                              cmap=cmap,\n",
    "                                              scale=0.3, xlabel='IWV MOD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)', by=['time.month'])\n",
    "\n",
    "ThaaoVsMYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width,\n",
    "#                                              cmap=cmap,\n",
    "                                              scale=0.3, xlabel='IWV MYD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              by=['time.month'])\n",
    "\n",
    "(ThaaoVsMOD_scatter*slope*regress_MOD_sk[0]*regress_MOD_sk_text\n",
    "+ ThaaoVsMYD_scatter*slope*regress_MYD_sk[0]*regress_MYD_sk_text\n",
    "+ (ThaaoVsMOD_scatter*ThaaoVsMYD_scatter*slope*regress_MOD_sk[0]*regress_MYD_sk[0]\n",
    "   *regress_TOT_sk[0]*regress_TOT_sk_text \n",
    "  ).opts(title='DDDDD', xlabel = 'IWV MXD (cm)')).cols(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7eac32c-9671-4fc1-baf0-d87e84efdf54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ok stagionalità\n",
    "#cmap = 'hsv'\n",
    "ThaaoVsMOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width,\n",
    "#                                              cmap=cmap,\n",
    "                                              scale=0.3, xlabel='IWV MOD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)', by=['time.month'])\n",
    "\n",
    "ThaaoVsMYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_MODIS', y='IWV_THAAO',\n",
    "                                              height=height, width=width,\n",
    "#                                              cmap=cmap,\n",
    "                                              scale=0.3, xlabel='IWV MYD (cm)',\n",
    "                                              ylabel='IWV THAAO (cm)',\n",
    "                                              by=['time.month'])\n",
    "\n",
    "\n",
    "(ThaaoVsMOD_scatter*ThaaoVsMYD_scatter\n",
    "# *slope\n",
    "#*regress_MOD_sk[0]*regress_MYD_sk[0]\n",
    "*regress_TOT_sk[0]*regress_TOT_sk_text\n",
    ").opts(title='IWV seasonality', xlabel = 'IWV MXD (cm)', show_grid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0e8a856-fb05-49f6-b4f2-bd52d775c9f4",
   "metadata": {},
   "source": [
    "### Diff IWV MXD-THAAO vs IWV THAAO "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bb712a8-a303-4d63-bfaa-09bbf80b223c",
   "metadata": {},
   "outputs": [],
   "source": [
    "height=380\n",
    "width=1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ade3e0ba-9607-4c09-a3aa-ba6b00ddc76d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Diff_iwv_Vs_MOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_THAAO',\n",
    "                                                   y='diff_IWV', by=['time.month'],\n",
    "                                                   xlabel='IWV VESPA (cm)',\n",
    "                                                   ylabel='Diff MOD-VESPA (cm)',\n",
    "                                                   height=height, width=width, \n",
    "                                                   scale=scale,grid=True)\n",
    "\n",
    "Diff_iwv_Vs_MYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_THAAO',\n",
    "                                                   y='diff_IWV', by=['time.month'],\n",
    "                                                   xlabel='IWV VESPA (cm)',\n",
    "                                                   ylabel='Diff MYD-VESPA (cm)',\n",
    "                                                   height=height, width=width, \n",
    "                                                   scale=scale,grid=True)\n",
    "\n",
    "(Diff_iwv_Vs_MOD_scatter\n",
    "+ Diff_iwv_Vs_MYD_scatter\n",
    "+ (Diff_iwv_Vs_MOD_scatter*Diff_iwv_Vs_MYD_scatter).opts(ylabel='Diff MXD-VESPA (cm)')).cols(1)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d0af74c0-9632-4eef-accb-0829fb704a89",
   "metadata": {},
   "source": [
    "pippo = (Diff_iwv_Vs_MOD_scatter*Diff_iwv_Vs_MYD_scatter).opts(ylabel='Diff MXD-THAAO (cm)')\n",
    "pippo = pn.panel(pippo)\n",
    "pippo.save('pippo.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40600c0c-f34e-4f98-9e25-7f5f8b8bea68",
   "metadata": {},
   "outputs": [],
   "source": [
    "title = f'IWV MXD-VESPA              deltatime={deltatime}      # data points={number_MXD_data_points_deltatime}'\n",
    "((Diff_iwv_Vs_MOD_scatter*Diff_iwv_Vs_MYD_scatter).\n",
    " opts(title=title, ylabel='Diff MXD-VESPA (cm)'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfaa6b1c-1016-4d61-a5ac-5bafe2086881",
   "metadata": {},
   "source": [
    "### Diff % IWV MXD-THAAO vs IWV THAAO "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a02ce4d-1aa5-4107-87bc-2ff7190642b0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Diffperc_iwv_Vs_MOD_scatter = newdf_MOD.hvplot.scatter(x='IWV_THAAO',\n",
    "                                                   y='diff_perc_IWV', by=['time.month'],\n",
    "                                                   xlabel='IWV VESPA (cm)',\n",
    "                                                   ylabel='100*(MOD-VESPA)/VESPA',\n",
    "                                                   height=height, width=width, \n",
    "                                                   scale=scale,grid=True)\n",
    "\n",
    "Diffperc_iwv_Vs_MYD_scatter = newdf_MYD.hvplot.scatter(x='IWV_THAAO',\n",
    "                                                   y='diff_perc_IWV', by=['time.month'],\n",
    "                                                   xlabel='IWV VESPA (cm)',\n",
    "                                                   ylabel='100*(MYD-VESPA)/VESPA',\n",
    "                                                   height=height, width=width, \n",
    "                                                   scale=scale,grid=True)\n",
    "\n",
    "(Diffperc_iwv_Vs_MOD_scatter\n",
    "+ Diffperc_iwv_Vs_MYD_scatter\n",
    "+ (Diffperc_iwv_Vs_MOD_scatter*Diffperc_iwv_Vs_MYD_scatter)\n",
    " .opts(ylabel='100*(MXD-VESPA)/VESPA')).cols(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28c7def1-1751-441e-8327-0e93a92547f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "title = f'IWV 100*(MXD-VESPA)/VESPA              deltatime={deltatime}      # data points={number_MXD_data_points_deltatime}'\n",
    "((Diffperc_iwv_Vs_MOD_scatter*Diffperc_iwv_Vs_MYD_scatter)\n",
    " .opts(title=title, ylabel='100*(MXD-VESPA)/VESPA'))"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a50faca6-1d49-446c-b471-4d697d728834",
   "metadata": {},
   "source": [
    "type(read_hatpro_year('/home/luca/WORK/DATA/Groundbased/THULE/HATPRO/IWV_2017_HATPRO.txt'))"
   ]
  }
 ],
 "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
}