ARPA_functionalities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# ARPA Lombardia - Meteorological Stations Data Request and Processing Notebook\n",
    "\n",
    "Author: Emanuele Capizzi - Email: emanuele.capizzi@polimi.it\n",
    "\n",
    "This notebook allows to access to ARPA Lombardia meteorological stations using the datasets freely available on [Open Data Regione Lombardia Catalog](https://dati.lombardia.it/). This notebook allowed to test functionalities to be used in this QGIS plugin [ARPA Weather](https://github.com/gisgeolab/ARPA_Weather_plugin)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Datasets organization and description\n",
    "\n",
    "Firstly, it is necessary to understand the organization of the available datasets.<br>\n",
    "The dataset that will be used in this notebook are:\n",
    "1. Sensors information\n",
    "2. Time-series measured by each sensor\n",
    "\n",
    "By merging this datasets on the Sensor ID column it is possible to combine the information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ieKff42Du_uf",
    "tags": []
   },
   "source": [
    "### 1. Sensors Information\n",
    "\n",
    "The information about each sensor in ARPA's weather monitoring network (~ 1200/1300 sensors) are available at the following can be obtained directly from the API (Application Programming Interface), allowing a systematic way to request the data from the catalog. These information are provided separately from the sensors measured time-series.<br> \n",
    "<br> \n",
    "Link to the resource: https://www.dati.lombardia.it/Ambiente/Stazioni-Meteorologiche/nf78-nj6b\n",
    "\n",
    "\n",
    "The sensors information accessible through this portal include: <br>\n",
    "\n",
    "| Information description | Column Name |\n",
    "| :---: | :---: |\n",
    "| ID of the sensor | idsensore |\n",
    "| Type of the Sensor | tipologia |\n",
    "| Measure Unit | unit_dimisura |\n",
    "| Station ID (multiple sensors can be available at the same station) | idstazione |\n",
    "| Name of the Station | nomestazione |\n",
    "| Height of the station (m) | quota |\n",
    "| Start date of the time-series | datastart |\n",
    "| Whether the sensor is historical | storico |\n",
    "| WGS84 UTM 32N North Coordinate | cgb_nord |\n",
    "| WGS84 UTM 32N East Coordinate | cgb_est |\n",
    "| Longitude | lng |\n",
    "| Latitude | lat |\n",
    "| Dictionary containing latitude and longitude | location |\n",
    "| Eventual end date of the time-series, if the sensor is not measuring | datastop |\n",
    "\n",
    "Specifically, the following  are the types of sensors available in the ARPA Meteorological Network:<br>\n",
    "\n",
    "| Sensor Type | Sensor Name | Measure Unit |\n",
    "| :---: | :---: | :---: |\n",
    "| Rainfall | Precipitazione | mm |\n",
    "| Temperature | Temperatura | °C |\n",
    "| Relative Humidity | Umidità Relativa | % |\n",
    "| Global Solar Radiation | Radiazione Globale | W/m<sup>2</sup> |\n",
    "| Wind Direction | Direzione Vento | Degrees North (°) |\n",
    "| Wind Speed | Velocità Vento | m/s |\n",
    "| Hydrometric Level | Livello Idrometrico | cm |\n",
    "| Snow Depth | Altezza Neve | cm |\n",
    "\n",
    "### 2. Measured Time-series\n",
    "\n",
    "Subsequently, it is possible to request the time-series for each sensor. The temporal frequency of the measurement may vary depending on the sensor.\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">  \n",
    "<b>For accessing such data, it should be pointed out that:</b> \n",
    "    <br><li>The time-series obtained from the API are updated regularly and can be accessed in real-time, but only from the beginning of the current month.\n",
    "    <br><li>For historical data (past months and years), it's required to download as a zip folder containing CSV files containing yearly observations.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-success\"> This notebook allows to handles it automatically, requesting data both from API or downloading the CSV files.\n",
    "</div>\n",
    "\n",
    "Link to the resource (API) containing sensors data: https://www.dati.lombardia.it/Ambiente/Dati-sensori-meteo/i95f-5avh <br>\n",
    "<br>\n",
    "The following table describe the time-series data obtainable through the API or the CSV files:<br>\n",
    "\n",
    "| Information description | Column Name | Notes |\n",
    "| :---: | :---: | :---: |\n",
    "| ID of the sensor | IdSensore | Corresponding to the ID available in the sensors information table |\n",
    "| Date of observation | Data | - |\n",
    "| Observation value | Valore | LEGEND: -9999 = missing data / 888, 8888 = variable wind direction / 777, 7777 = quiet (for wind direction only) |\n",
    "| ID of the Operator | IdOperatore | LEGEND: 1: Mean value / 3: Maximum Value / 4: Valore cumulated value (for rain) |\n",
    "| Data status | Stato | LEGEND: VA, VV = Validated data / NA, NV, NC = Invalid data / NI = Uncertain data / ND = Data not available |\n",
    "\n",
    "\n",
    "### Notes about the notebook\n",
    "The functions used in this notebook are stored in a separated .py file that is imported at the beginning of the notebook, in order to keep the notebook as clean as possible.<br>\n",
    "In order to see the widgets it may be necessary to update ipywidgets library or to use the following command in the Command Line: ```jupyter nbextension enable --py widgetsnbextension```\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Access to ARPA Lombardia data through Socrata API\n",
    "### 1. Socrata API\n",
    "The [Socrata Open Data API](https://dev.socrata.com/) allows you to programmatically access a wealth of open data resources from governments, non-profits, and NGOs around the world. <br>\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">  \n",
    "<b>For accessing such data, it should be pointed out that:</b> \n",
    "While it is possible to perform simple unauthenticated queries against the Socrata Open Data API without making use of an application token, you’ll receive much higher throttling limits if you include an application token in your requests.\n",
    "    <b>You can try to request data without using a Token and check if it works, but in any case the procedure to obtain a token is described.</b>\n",
    "</div>\n",
    "\n",
    "### 2. Sodapy\n",
    "In this notebook [sodapy](https://github.com/xmunoz/sodapy) Python library is used. Sodapy is a python client for the Socrata Open Data API.\n",
    "\n",
    "### 3. How to obtain a token\n",
    " - In order to get a token it is necessary to open [Open Data Lombardia](https://dati.lombardia.it/) website.\n",
    " - Subscribe to the website and go to your profile settings by clicking on your name in the upper-right corner. Then click on  `Il mio profilo ` tab.\n",
    " - Once your are on your profile, click on the  `Pen symbol` near the  `Il tuo profilo`, as shown here:<br>\n",
    " <br>\n",
    " <img src=\"./img/ARPA_API.png\" style=\"display: block; margin: auto\";/><br>\n",
    " <br>\n",
    "- Modify your profile and open the  `Opzioni per lo sviluppatore` tab. Create a new  `App Token` to be be used by clicking  `Crea una nuova applicazione`, as shown: <br> \n",
    " <br>\n",
    " <img src=\"./img/developer_settings.png\"style=\"display: block; margin: auto\";/><br>\n",
    " <br>\n",
    " \n",
    "- Create your Token by inserting all the required information, as shown:<br> \n",
    " <br>\n",
    " <img src=\"./img/app_token_modify.png\" width=\"400px\"  style=\"display: block; margin: auto\";/><br>\n",
    " <br>\n",
    "\n",
    "- Finally the token will be available in the  `Token App` table, as shown:<br> \n",
    " <br>\n",
    " <img src=\"./img/final_token.png\" width=\"900px\"  style=\"display: block; margin: auto\";/><br>\n",
    " <br>\n",
    "\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">  \n",
    "<b><span>&#x2714;</span></b> Now you can copy and paste the token to be used for accessing the Socrata API.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries\n",
    "The most relevant libraries to be installed for this notebook are the following:\n",
    "- sodapy\n",
    "- pandas\n",
    "- dask\n",
    "- ipywidgets\n",
    "- matplotlib\n",
    "- seaborn\n",
    "- geopandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vIYOzSgg1Sji",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Filter warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# Libraries\n",
    "from sodapy import Socrata\n",
    "import pandas as pd\n",
    "import dask.dataframe as dd\n",
    "import numpy as np\n",
    "import ipywidgets as widgets\n",
    "import geopandas as gpd \n",
    "import requests\n",
    "from datetime import datetime, timedelta\n",
    "from io import BytesIO\n",
    "from zipfile import ZipFile\n",
    "import os\n",
    "import time\n",
    "\n",
    "# Plotting \n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,AutoMinorLocator)\n",
    "import matplotlib.dates as mdates\n",
    "import seaborn as sns\n",
    "\n",
    "# Import functions and set auto-reload\n",
    "import ARPA_functions as f\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Client Connection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function will create a client connection to ARPA API. It will use the token if it's provided by the user, otherwise it will proceed without using the token. <b> Meteo data can be requested withouth Token for now, withouth any evident limitation. This step can be skipped. </b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "token_text_w = widgets.Text(\n",
    "    placeholder='Enter token here',\n",
    "    description='Enter token (if available):',\n",
    "    disabled=False,\n",
    "    style= {'description_width': 'initial'},\n",
    "    layout = widgets.Layout(width='400px'))\n",
    "token_text_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 223
    },
    "executionInfo": {
     "elapsed": 266,
     "status": "error",
     "timestamp": 1673873146302,
     "user": {
      "displayName": "geolab como",
      "userId": "12419660887110986186"
     },
     "user_tz": -60
    },
    "id": "CWWE66sl02vk",
    "outputId": "51b5575c-1d89-4077-ae00-e0b277c0cb72",
    "tags": []
   },
   "outputs": [],
   "source": [
    "arpa_token = token_text_w.value  # Use empty string to use the API without token\n",
    "print(arpa_token)\n",
    "client = f.connect_ARPA_api(arpa_token)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "stationsId = \"nf78-nj6b\" # Select meteo stations dataset containing positions and information about sensors\n",
    "sensors_info = client.get_all(stationsId)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieve sensors information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As explained at the beginning of the notebook, it is necessary to retrieve all sensors information. The following function allows to obtain and print them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sensors_df = f.ARPA_sensors_info(client)\n",
    "sensors_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can select the sensor type from the following list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "unique_sensors_list = sensors_df['tipologia'].unique()\n",
    "sw = widgets.Dropdown(\n",
    "    options=unique_sensors_list,\n",
    "    value='Temperatura',\n",
    "    description='Sensor type:')\n",
    "sw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sensor_sel = sw.value # Store selected sensor value\n",
    "sensor_sel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a list of sensors that belong to the selected sensor type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sensors_list = (sensors_df.loc[sensors_df['tipologia'] == sensor_sel]).idsensore.tolist()  #& (sensors_df['storico'] == storic_data)\n",
    "print((\"Selected sensor: {sel}\").format(sel=sensor_sel))\n",
    "print((\"Number of selected sensor: {sens_len}\").format(sens_len=len(sensors_list)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Rt3NMBcH6osR"
   },
   "source": [
    "## Check API Time-Series Availability\n",
    "\n",
    "The following function is used to obtain minimum and maximum date of the corresponding time-series available in the ARPA API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "start_date_API, end_date_API = f.req_ARPA_start_end_date_API(client)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "api_start_limit = datetime(datetime.today().year, datetime.today().month, 1)\n",
    "api_start_limit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Date Range"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section is possible to select the time range for data request. \n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">  \n",
    "<li> Only dates within the same year can be selected.\n",
    "<li>  If dates from previous month are requested the CSV file of the corresponding file will be downloaded and used for processing.\n",
    "<li> The CSV are available from 2013 onwards.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create a date picker\n",
    "start_picker = widgets.DatetimePicker(description='Start Date:')\n",
    "end_picker = widgets.DatetimePicker(description='End Date:')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Display the date and time pickers\n",
    "widgets.VBox([start_picker, end_picker])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "start_date = start_picker.value.replace(tzinfo=None)\n",
    "end_date = end_picker.value.replace(tzinfo=None)\n",
    "\n",
    "# Check if dates are ok\n",
    "year, start_date, end_date = f.check_dates(start_date, end_date)\n",
    "print(\"Year:\", year,\"/ Start date:\", start_date,\"/ End date:\", end_date)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download and process observations from API or csv files\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">  \n",
    "If you see a limitation on data download you can launch the notebook using the following command: \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#!jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download and process observations from API or csv files\n",
    "\n",
    "The following code block will check if the chosen start date is before the start date available in the API.\n",
    "If if the date is before that date, the CSV file corresponding to the selected year will be downloaded and processed.\n",
    "For example, if the current  date is 20 February 2023 and you request data from 15 January 2023 to 15 February 2023:\n",
    "- the 2023 CSV data will be downloaded and processed up to 31 January 2023.\n",
    "- if you need the data from 01 February 2023 onwards, you need to set it as start date in order to request the data from the API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#If the chosen start date is before the start date available in the API -> use csv data\n",
    "if start_date < api_start_limit:\n",
    "    print(\"Requesting CSV\")\n",
    "    sensors_values = f.download_extract_csv_from_year(str(year)) #download the csv corresponding to the selected year\n",
    "    csv_file = str(year)+'.csv'\n",
    "    sensors_values = f.process_ARPA_csv(csv_file, start_date, end_date, sensors_list) #process csv file with dask\n",
    "    \n",
    "#If the chosen start date is equal or after the start date of API -> request data from API\n",
    "elif start_date >= api_start_limit:\n",
    "    print(\"Requesting from API\")\n",
    "    sensors_values = f.req_ARPA_data_API(client, start_date, end_date, sensors_list) #request data from ARPA API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sensors_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Filter out outliers\n",
    "Some filtering functions can be applied to remove invalid values and outlier, such as:\n",
    "- Interquantile range (IQR): Outliers are defined as observations that fall below Q1 − 1.5 IQR or above Q3 + 1.5 IQR.\n",
    "- Z-Score: calculate the Z-Score for the observations and remove those below a given threshold. (normally distributed data)\n",
    "\n",
    "Select the method to be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outlier_methods = ['iqr', 'zscore'] # List of methods for removing outliers\n",
    "checkboxes = [widgets.Checkbox(value=False, description=label) for label in outlier_methods]\n",
    "output = widgets.VBox(children=checkboxes)\n",
    "display(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "selected_data = []\n",
    "for i in range(0, len(checkboxes)):\n",
    "    if checkboxes[i].value == True:\n",
    "        selected_data = selected_data + [checkboxes[i].description]\n",
    "print(selected_data)\n",
    "if 'iqr' in selected_data:\n",
    "    sensors_values = sensors_values.groupby('idsensore').apply(f.outlier_filter_iqr)\n",
    "    print('Removed outliers using IQR')\n",
    "if 'zscore' in selected_data:\n",
    "    sensors_values = f.outlier_filter_zscore(sensors_values, sensors_list)\n",
    "    print('Removed outliers using zscore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data aggregation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to aggregate the data in the dataframe saved in the `sensors_values` dataframe.\n",
    "Must provide:\n",
    "- dataframe containing idsensore (int), data (datetime), valore (float).\n",
    "- Temporal aggregation: hour (H), day (D), week (W), month (M), year (Y).\n",
    "\n",
    "The statistics calculated are:\n",
    "- Mode, count $\\rightarrow$ for wind direction (since it's expressed in North Degrees).\n",
    "- Mean, min, max, std, count $\\rightarrow$ for all the remaining sensors.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the function to aggregate and group the data. You can pass the value obtained from the API or the values obtaines from the csv.\n",
    "\n",
    "Choose the temporal aggregation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tagg_w = widgets.ToggleButtons(\n",
    "    options=['H', 'D', 'W', 'M', 'Y'],\n",
    "    description='Period aggregation:',\n",
    "    disabled=False,\n",
    "    button_style='info', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltips=['Hour', 'Day', 'Week', 'Month', 'Year'])\n",
    "tagg_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tagg = tagg_w.value\n",
    "\n",
    "if sensor_sel != 'Direzione Vento':\n",
    "    sensor_test_agg = f.aggregate_group_data(sensors_values, tagg)\n",
    "if sensor_sel == 'Direzione Vento':\n",
    "    sensor_test_agg = f.aggregate_group_data_wind_dir(sensors_values, tagg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensor_test_agg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Join sensors information and time series \n",
    "Once the dataframes are created is possible to merge the information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_df = pd.merge(sensor_test_agg, sensors_df, on='idsensore')\n",
    "measure_unit = merged_df['unit_dimisura'].unique()[0]\n",
    "print(\"Measure unit: \" + measure_unit)\n",
    "merged_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if sensor_sel == 'Direzione Vento':\n",
    "    stats_list = ['mode', 'count']\n",
    "    default_value = 'mode'\n",
    "if sensor_sel != 'Direzione Vento':\n",
    "    stats_list = ['mean', 'max', 'min', 'std', 'count']\n",
    "    default_value = 'mean'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stat_w = widgets.Dropdown(\n",
    "    options=stats_list,\n",
    "    value= default_value,\n",
    "    description='Value:')\n",
    "stat_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "stat_sel = stat_w.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot all sensors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new figure\n",
    "fig, ax = plt.subplots(figsize=(15, 10))\n",
    "ax.set_title('ARPA Lombardia Sensors - ' + sensor_sel, fontdict={'fontsize': 10})\n",
    "ax.set_ylabel(measure_unit, fontdict={'fontsize': 10})\n",
    "ax.set_xlabel('Date', fontdict={'fontsize': 10})\n",
    "\n",
    "# Iterate over the sensor IDs\n",
    "for sensor_id in sensors_list:\n",
    "    # Get the data for the current sensor\n",
    "    sensor_data = merged_df[(merged_df['idsensore'] == sensor_id)]  # add & (merged_df['provincia'] == 'MI') in the square brackets if want to filter on Milan\n",
    "    # Plot the time series for the sensor\n",
    "    ax.plot(sensor_data['data'], sensor_data[stat_sel], label=sensor_id)\n",
    "\n",
    "\n",
    "# Show the plot\n",
    "plt.gca().xaxis.set_major_locator(mdates.DayLocator(interval=10))\n",
    "plt.gcf().autofmt_xdate()\n",
    "plt.xticks(rotation=90)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check highest and lowest values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select the sensor that provide the lowest value in the dataframe to check if there are possible errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_index = merged_df[stat_sel].idxmin()\n",
    "idsensor_min = merged_df.loc[min_index].idsensore\n",
    "sel_sensor_min = merged_df.loc[merged_df['idsensore']==idsensor_min]\n",
    "data_min = merged_df.loc[min_index]['data']\n",
    "merged_df.loc[min_index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select the sensor that provide the highest value in the dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_index = merged_df[stat_sel].idxmax()\n",
    "idsensor_max = merged_df.loc[max_index].idsensore\n",
    "sel_sensor_max = merged_df.loc[merged_df['idsensore']==idsensor_max]\n",
    "data_max = merged_df.loc[max_index]['data']\n",
    "merged_df.loc[max_index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot min and max selected sensors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,8))\n",
    "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))\n",
    "plt.gca().xaxis.set_major_locator(mdates.DayLocator())\n",
    "\n",
    "plt.gca().xaxis.set_major_locator(mdates.DayLocator(interval=10))\n",
    "plt.gcf().autofmt_xdate()\n",
    "plt.xticks(rotation=90)\n",
    "plt.plot(sel_sensor_min['data'],sel_sensor_min[stat_sel], '-',  color=\"tab:blue\", label=(\"Min time-series - {nome} {quota} m\").format(nome=sel_sensor_min[\"nomestazione\"].unique()[0], quota=sel_sensor_min[\"quota\"].unique()[0]) + \" - \" + stat_sel)\n",
    "plt.plot(sel_sensor_max['data'],sel_sensor_max[stat_sel], '-',  color=\"tab:red\", label=(\"Max time-series - {nome} {quota} m\").format(nome=sel_sensor_max[\"nomestazione\"].unique()[0], quota=sel_sensor_max[\"quota\"].unique()[0]) + \" - \" + stat_sel)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Box Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_pal = {0: \"tab:red\", 1: \"tab:blue\"}\n",
    "box = sns.boxplot(data=[sel_sensor_max[stat_sel], sel_sensor_min[stat_sel]], orient=\"v\", palette=my_pal)\n",
    "plt.xticks([0, 1],[\"Max time-series\",\"Min time-series\"])\n",
    "box.annotate((\"ID: {id}, Quota: {quota} m\\n {nome}\\n {stat}\").format(id=sel_sensor_max[\"idsensore\"].unique()[0],\n",
    "             quota=sel_sensor_max[\"quota\"].unique()[0], nome=sel_sensor_max[\"nomestazione\"].unique()[0], stat=stat_sel),\n",
    "             xy=(0, sel_sensor_max[stat_sel].max()+1),\n",
    "             horizontalalignment='center',\n",
    "             bbox=dict(boxstyle=\"round\", fc=\"0.8\"))\n",
    "box.annotate((\"ID: {id}, Quota: {quota} m\\n {nome}\\n {stat}\").format(id=sel_sensor_min[\"idsensore\"].unique()[0],\n",
    "             quota=sel_sensor_min[\"quota\"].unique()[0], nome=sel_sensor_min[\"nomestazione\"].unique()[0], stat=stat_sel),\n",
    "             xy=(1, sel_sensor_max[stat_sel].max()+1),\n",
    "             horizontalalignment='center',\n",
    "             bbox=dict(boxstyle=\"round\", fc=\"0.8\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import Lombardy and CMM vector files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lombardy_gdf = gpd.read_file('lombardia_boundary.gpkg')\n",
    "cmm_gdf = gpd.read_file('CMM.gpkg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert to Geodataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = gpd.GeoDataFrame(merged_df, geometry=gpd.points_from_xy(merged_df.cgb_est, merged_df.cgb_nord), crs=\"EPSG:32632\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter on the dates corresponding to the max and min values found before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_date = merged_df.loc[min_index].data\n",
    "max_date = merged_df.loc[max_index].data\n",
    "min_date, max_date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the gdf corresponding with the date where the maximum and minimum values are observed\n",
    "max_gdf = gdf.loc[merged_df['data'] == max_date]\n",
    "min_gdf = gdf.loc[merged_df['data'] == min_date]\n",
    "\n",
    "# Get the point (positions) where the minimum and maximum are observed\n",
    "min_point = min_gdf.loc[min_gdf['idsensore'] == idsensor_min]\n",
    "max_point = max_gdf.loc[max_gdf['idsensore'] == idsensor_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (9,9))\n",
    "min_gdf.plot(ax=ax, markersize=20, column=stat_sel, legend=True, cmap='coolwarm', legend_kwds={'label': measure_unit ,'orientation': \"vertical\"});\n",
    "lombardy_gdf.boundary.plot(ax=ax);\n",
    "cmm_gdf.boundary.plot(ax=ax, edgecolor='orange');\n",
    "min_point.plot(ax=ax, marker='o', markersize=400, color='black',facecolors='none')\n",
    "ax.set_title((\"ID: {id} - Quota: {quota} m - {nome} T°: {minT:0.2f}\\n Data: {data} - Min - ARPA Lombardia Sensors - {sensor_sel}\").format(id=sel_sensor_min[\"idsensore\"].unique()[0],\n",
    "             quota=sel_sensor_min[\"quota\"].unique()[0], nome=sel_sensor_min[\"nomestazione\"].unique()[0], minT=sel_sensor_min[stat_sel].min(),\n",
    "             data=data_min.strftime(\"%Y-%m-%d\"), sensor_sel=sensor_sel), fontdict = {'fontsize': 10})\n",
    "ax.set_ylabel('North [m]',fontdict = {'fontsize': 10})\n",
    "ax.set_xlabel('East [m]',fontdict = {'fontsize': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (9,9))\n",
    "max_gdf.plot(ax=ax, markersize=20, column=stat_sel, legend=True, cmap='coolwarm', legend_kwds={'label': measure_unit ,'orientation': \"vertical\"});\n",
    "lombardy_gdf.boundary.plot(ax=ax);\n",
    "cmm_gdf.boundary.plot(ax=ax, edgecolor='orange');\n",
    "max_point.plot(ax=ax, marker='o', markersize=400, color='black',facecolors='none')\n",
    "ax.set_title((\"ID: {id} - Quota: {quota} m -  {nome} T°: {maxT:0.2f}\\n Data: {data} - Max - ARPA Lombardia Sensors - {sensor_sel}\").format(id=sel_sensor_max[\"idsensore\"].unique()[0],\n",
    "             quota=sel_sensor_max[\"quota\"].unique()[0], nome=sel_sensor_max[\"nomestazione\"].unique()[0], maxT=sel_sensor_max[stat_sel].max(),\n",
    "             data=data_max.strftime(\"%Y-%m-%d\"), sensor_sel=sensor_sel), fontdict = {'fontsize': 10})\n",
    "\n",
    "ax.set_ylabel('North [m]',fontdict = {'fontsize': 10})\n",
    "ax.set_xlabel('East [m]',fontdict = {'fontsize': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (10,10))\n",
    "max_gdf.plot(ax=ax, markersize=20, column='quota', legend=True, cmap='coolwarm', legend_kwds={'label': 'm' ,'orientation': \"vertical\"});\n",
    "lombardy_gdf.boundary.plot(ax=ax);\n",
    "cmm_gdf.boundary.plot(ax=ax, edgecolor='orange');\n",
    "ax.set_title('ARPA Lombardia Sensors - Quota',fontdict = {'fontsize': 10})\n",
    "ax.set_ylabel('North [m]',fontdict = {'fontsize': 10})\n",
    "ax.set_xlabel('East [m]',fontdict = {'fontsize': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_test = gpd.GeoDataFrame(merged_df, geometry=gpd.points_from_xy(merged_df.lng, merged_df.lat), crs=\"EPSG:4326\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_test.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_test.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_gdf[\"lat\"] = [float(str(i).replace(\",\", \"\")) for i in min_gdf[\"lat\"]]\n",
    "min_gdf[\"lng\"] = [float(str(i).replace(\",\", \"\")) for i in min_gdf[\"lng\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_gdf[\"lat\"] = [float(str(i).replace(\",\", \"\")) for i in max_gdf[\"lat\"]]\n",
    "max_gdf[\"lng\"] = [float(str(i).replace(\",\", \"\")) for i in max_gdf[\"lng\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#using plotly for an animated choropleth map\n",
    "import plotly.express as px\n",
    "fig = px.scatter_mapbox(min_gdf,\n",
    "                        lat=\"lat\",\n",
    "                        lon=\"lng\",\n",
    "                        hover_name=\"idsensore\",\n",
    "                        hover_data=[\"nomestazione\",\"provincia\",stat_sel],\n",
    "                        color=stat_sel,\n",
    "                        color_continuous_scale=px.colors.cyclical.IceFire,\n",
    "                        zoom=7,\n",
    "                        height=600,\n",
    "                        size='lat',\n",
    "                        size_max=12,\n",
    "                        opacity=1,\n",
    "                        width=900)\n",
    "fig.update_layout(mapbox_style='carto-darkmatter')\n",
    "fig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter_mapbox(max_gdf,\n",
    "                        lat=\"lat\",\n",
    "                        lon=\"lng\",\n",
    "                        hover_name=\"idsensore\",\n",
    "                        hover_data=[\"nomestazione\",\"provincia\",stat_sel],\n",
    "                        color=stat_sel,\n",
    "                        color_continuous_scale=px.colors.cyclical.IceFire,\n",
    "                        zoom=7,\n",
    "                        height=600,\n",
    "                        size='lat',\n",
    "                        size_max=12,\n",
    "                        opacity=1,\n",
    "                        width=900)\n",
    "fig.update_layout(mapbox_style='carto-darkmatter')\n",
    "fig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "1a884f121133f7e48c0fb0d37d096822567e6dace2ca2327981e34c414aeb97b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}