VRI_anomaly_detection_bites.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "48ea4439-32c5-4d32-8d51-e4b90e49466c",
   "metadata": {},
   "source": [
    "# VRI anomaly detection bites"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e858118-e213-4f67-93be-0ba957abfa24",
   "metadata": {},
   "source": [
    "This notebook has the aim to study how to detect anomalies in the VRI computed by our models displayed here: https://labs.mosquitoalert.com/MosquitoAlertES/\n",
    "\n",
    "Data gathered from models bites."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "238d783b-8961-46da-965d-bc259227c90f",
   "metadata": {},
   "source": [
    "## Requirements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "be769fb1-d114-425b-b875-ffea20848d2e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Importing plotly failed. Interactive plots will not work.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from prophet import Prophet\n",
    "from prophet.plot import seasonality_plot_df\n",
    "import os\n",
    "import pandas as pd\n",
    "import json\n",
    "from datetime import datetime\n",
    "import re\n",
    "import geopandas as gpd\n",
    "from tqdm import tqdm\n",
    "from scipy.signal import savgol_filter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6dc1cd3-51ea-435d-bbbc-d2cb58f5da07",
   "metadata": {},
   "source": [
    "## Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c62acf6c-ade0-4674-8c79-079192bd2eb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "BITING_DATA_PATH = './biting/'\n",
    "\n",
    "# Collect all CSV file paths\n",
    "files = [os.path.join(root, file)\n",
    "        for root, _, files in os.walk(BITING_DATA_PATH)\n",
    "        for file in files if file == 'gadm4_monthly.csv']\n",
    "\n",
    "# Initialize an empty list to hold the data\n",
    "dfs = []\n",
    "# Loop through the files\n",
    "for file in files:\n",
    "    try:\n",
    "        dfs.append(pd.read_csv(file))\n",
    "    except Exception as e:\n",
    "        print(f\"Error processing file {file}: {e}\")\n",
    "\n",
    "# Create a DataFrame from the list of data\n",
    "df = pd.concat(dfs)\n",
    "del dfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "de4f8b00-2a9f-437a-8645-ccec1a58c1f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>year</th>\n",
       "      <th>month</th>\n",
       "      <th>gid_4</th>\n",
       "      <th>est</th>\n",
       "      <th>se</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2014</td>\n",
       "      <td>5</td>\n",
       "      <td>ALA.2.4_2</td>\n",
       "      <td>0.155015</td>\n",
       "      <td>0.719975</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2014</td>\n",
       "      <td>5</td>\n",
       "      <td>ALA.2.5_2</td>\n",
       "      <td>0.256332</td>\n",
       "      <td>0.869462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2014</td>\n",
       "      <td>5</td>\n",
       "      <td>ALA.3.1_2</td>\n",
       "      <td>0.252280</td>\n",
       "      <td>0.864868</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2014</td>\n",
       "      <td>5</td>\n",
       "      <td>ALB.1.1.1_1</td>\n",
       "      <td>0.541033</td>\n",
       "      <td>0.997685</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2014</td>\n",
       "      <td>5</td>\n",
       "      <td>ALB.1.1.10_1</td>\n",
       "      <td>0.526849</td>\n",
       "      <td>0.999200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70826</th>\n",
       "      <td>2021</td>\n",
       "      <td>9</td>\n",
       "      <td>ZNC.1_1</td>\n",
       "      <td>0.693000</td>\n",
       "      <td>0.922499</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70827</th>\n",
       "      <td>2021</td>\n",
       "      <td>9</td>\n",
       "      <td>ZNC.2_1</td>\n",
       "      <td>0.747000</td>\n",
       "      <td>0.869462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70828</th>\n",
       "      <td>2021</td>\n",
       "      <td>9</td>\n",
       "      <td>ZNC.3_1</td>\n",
       "      <td>0.675000</td>\n",
       "      <td>0.936750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70829</th>\n",
       "      <td>2021</td>\n",
       "      <td>9</td>\n",
       "      <td>ZNC.4_1</td>\n",
       "      <td>0.713000</td>\n",
       "      <td>0.904723</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70830</th>\n",
       "      <td>2021</td>\n",
       "      <td>9</td>\n",
       "      <td>ZNC.5_1</td>\n",
       "      <td>0.506000</td>\n",
       "      <td>0.999928</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>9061729 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       year  month         gid_4       est        se\n",
       "0      2014      5     ALA.2.4_2  0.155015  0.719975\n",
       "1      2014      5     ALA.2.5_2  0.256332  0.869462\n",
       "2      2014      5     ALA.3.1_2  0.252280  0.864868\n",
       "3      2014      5   ALB.1.1.1_1  0.541033  0.997685\n",
       "4      2014      5  ALB.1.1.10_1  0.526849  0.999200\n",
       "...     ...    ...           ...       ...       ...\n",
       "70826  2021      9       ZNC.1_1  0.693000  0.922499\n",
       "70827  2021      9       ZNC.2_1  0.747000  0.869462\n",
       "70828  2021      9       ZNC.3_1  0.675000  0.936750\n",
       "70829  2021      9       ZNC.4_1  0.713000  0.904723\n",
       "70830  2021      9       ZNC.5_1  0.506000  0.999928\n",
       "\n",
       "[9061729 rows x 5 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "35db900a-6918-45cf-99be-0909cd5f0a9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ds</th>\n",
       "      <th>gid_4</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2990468</th>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.393110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2990469</th>\n",
       "      <td>2014-02-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.402229</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2990470</th>\n",
       "      <td>2014-03-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.315096</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2990471</th>\n",
       "      <td>2014-04-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.482270</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2990472</th>\n",
       "      <td>2014-05-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.426545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4028456</th>\n",
       "      <td>2024-03-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.219000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4028457</th>\n",
       "      <td>2024-04-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.208000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4028458</th>\n",
       "      <td>2024-05-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.354000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4028459</th>\n",
       "      <td>2024-06-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.442000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4028460</th>\n",
       "      <td>2024-07-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.513000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1037993 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                ds           gid_4         y\n",
       "2990468 2014-01-01  ESP.1.1.1.10_1  0.393110\n",
       "2990469 2014-02-01  ESP.1.1.1.10_1  0.402229\n",
       "2990470 2014-03-01  ESP.1.1.1.10_1  0.315096\n",
       "2990471 2014-04-01  ESP.1.1.1.10_1  0.482270\n",
       "2990472 2014-05-01  ESP.1.1.1.10_1  0.426545\n",
       "...            ...             ...       ...\n",
       "4028456 2024-03-01   ESP.9.1.9.9_1  0.219000\n",
       "4028457 2024-04-01   ESP.9.1.9.9_1  0.208000\n",
       "4028458 2024-05-01   ESP.9.1.9.9_1  0.354000\n",
       "4028459 2024-06-01   ESP.9.1.9.9_1  0.442000\n",
       "4028460 2024-07-01   ESP.9.1.9.9_1  0.513000\n",
       "\n",
       "[1037993 rows x 3 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Rename columns for Prophet\n",
    "df['ds'] = pd.to_datetime(df[['year', 'month']].assign(day=1))\n",
    "df.rename(columns={\"est\": \"y\"}, inplace=True)\n",
    "\n",
    "df.sort_values(by=['gid_4', 'ds'], inplace=True, ignore_index=True)\n",
    "\n",
    "# Keep only values for gid_4, ds, y\n",
    "df = df[['ds', 'gid_4', 'y']]\n",
    "\n",
    "# Keep only values for Spanish municipalities\n",
    "df = df[df['gid_4'].str.startswith('ESP')]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b9200e32-2596-4ee6-8243-d47d601a9a15",
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "logger = logging.getLogger('cmdstanpy')\n",
    "logger.addHandler(logging.NullHandler())\n",
    "logger.propagate = False\n",
    "logger.setLevel(logging.CRITICAL)\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=pd.errors.SettingWithCopyWarning)\n",
    "\n",
    "# Function to train a model and detect anomalies for each city\n",
    "def detect_anomalies_for_city(city_data):\n",
    "    group_name, city_df = city_data\n",
    "    if (city_df['y'] == 0).all():  # Skip if all original items are 0\n",
    "        return None, None\n",
    "\n",
    "    first_non_zero = city_df[city_df[\"y\"] != 0].iloc[0]\n",
    "    holidays_df = city_df[(city_df['y']==0) & (city_df['ds'] < first_non_zero['ds'])]['ds'].reset_index()\n",
    "    holidays_df['holiday'] = 'no-prediction-yet'\n",
    "    \n",
    "    # Step 3: Initialize Prophet with logistic growth\n",
    "    model = Prophet(growth='logistic', yearly_seasonality=True, weekly_seasonality=False, daily_seasonality=False, holidays=holidays_df[['ds','holiday']])\n",
    "    city_df.loc[:,'cap'] = 1\n",
    "    city_df.loc[:,'floor'] = 0\n",
    "    model.fit(city_df)\n",
    "    \n",
    "    # Make predictions for historical data (no future periods)\n",
    "    future = model.make_future_dataframe(periods=0)\n",
    "    future['cap'] = 1  # Ensure the future data has the cap\n",
    "    future['floor'] = 0  # Ensure the future data has the floor\n",
    "    forecast = model.predict(future)\n",
    "\n",
    "    forecast['fact'] = city_df['y'].reset_index(drop = True)\n",
    "\n",
    "    forecast['anomaly'] = 0\n",
    "    forecast.loc[forecast['fact'] > forecast['yhat_upper'], 'anomaly'] = 1\n",
    "    forecast.loc[forecast['fact'] < forecast['yhat_lower'], 'anomaly'] = -1\n",
    "\n",
    "     #anomaly importances\n",
    "    forecast['importance'] = 0.0\n",
    "    forecast.loc[forecast['anomaly'] ==1, 'importance'] = \\\n",
    "        (forecast['fact'] - forecast['yhat_upper'])/forecast['fact']\n",
    "    forecast.loc[forecast['anomaly'] ==-1, 'importance'] = \\\n",
    "        (forecast['yhat_lower'] - forecast['fact'])/forecast['fact']\n",
    "\n",
    "    # Merge forecast with the original data\n",
    "    city_df_forecast = forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper', 'trend', 'anomaly', 'importance']]\n",
    "    result_df = city_df[['gid_4', 'ds']].merge(city_df_forecast, on='ds', how='left')\n",
    "\n",
    "    # Seasonality component\n",
    "    df_w = seasonality_plot_df(m=model, ds=pd.date_range(start='2017-01-01', periods=365))\n",
    "    seas_df = model.predict_seasonal_components(df_w)\n",
    "    yearly_df = seas_df['yearly'].reset_index()\n",
    "    yearly_df.loc[:,'gid_4'] = city_df.iloc[0]['gid_4']\n",
    "\n",
    "    return result_df, yearly_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a4e0c809-24e5-4594-bb3d-0c67b9834e3a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████| 8205/8205 [02:32<00:00, 53.67it/s]\n"
     ]
    }
   ],
   "source": [
    "from concurrent.futures import ProcessPoolExecutor\n",
    "import os\n",
    "import math\n",
    "\n",
    "# Apply the anomaly detection for each city in parallel\n",
    "with ProcessPoolExecutor(max_workers=math.floor(max(os.cpu_count() * 0.8, 1))) as executor:\n",
    "    results = list(\n",
    "        tqdm(\n",
    "            executor.map(\n",
    "                detect_anomalies_for_city,\n",
    "                df.groupby('gid_4')\n",
    "            ), \n",
    "            total=len(\n",
    "                df['gid_4'].unique()\n",
    "            )\n",
    "        )\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c06ffedb-b3f1-4e25-b247-d99ddb5d12a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine the results for all cities\n",
    "result_df = df.merge(\n",
    "    pd.concat([arr[0] for arr in results if arr is not None]),\n",
    "    on=['gid_4', 'ds'], \n",
    "    how='left'\n",
    ")\n",
    "# Setting a 0 for the prediction value that hasn't been predicted because was all 0.\n",
    "result_df[['yhat', 'yhat_lower', 'yhat_upper', 'trend', 'anomaly', 'importance']] = result_df[['yhat', 'yhat_lower', 'yhat_upper', 'trend', 'anomaly', 'importance']].fillna(0)\n",
    "yearly_seasonality_df = pd.concat([arr[1] for arr in results if arr is not None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f46a6d27-768a-4ca3-bd60-e26dd3cd9148",
   "metadata": {},
   "outputs": [],
   "source": [
    "result_df.to_csv('spain_activty_anomaly_bites.csv', index=False)\n",
    "yearly_seasonality_df.to_csv('spain_seasonality_bites.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "460c2fa5-fe57-44b2-98ab-5774d12a8b79",
   "metadata": {},
   "source": [
    "## Part 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "28a60289-0f3d-41a1-a830-86a215aebc10",
   "metadata": {},
   "outputs": [],
   "source": [
    "result_df = pd.read_csv('./spain_activty_anomaly_bites.csv')\n",
    "yearly_seasonality_df = pd.read_csv('./spain_seasonality_bites.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7131f23e-d54e-41b6-a566-e5b9ff8f3f68",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ds</th>\n",
       "      <th>gid_4</th>\n",
       "      <th>y</th>\n",
       "      <th>yhat</th>\n",
       "      <th>yhat_lower</th>\n",
       "      <th>yhat_upper</th>\n",
       "      <th>trend</th>\n",
       "      <th>anomaly</th>\n",
       "      <th>importance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.393110</td>\n",
       "      <td>0.322870</td>\n",
       "      <td>0.250874</td>\n",
       "      <td>0.392607</td>\n",
       "      <td>0.528668</td>\n",
       "      <td>1</td>\n",
       "      <td>0.001281</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2014-02-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.402229</td>\n",
       "      <td>0.326708</td>\n",
       "      <td>0.257482</td>\n",
       "      <td>0.394605</td>\n",
       "      <td>0.527618</td>\n",
       "      <td>1</td>\n",
       "      <td>0.018955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2014-03-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.315096</td>\n",
       "      <td>0.365371</td>\n",
       "      <td>0.298259</td>\n",
       "      <td>0.437156</td>\n",
       "      <td>0.526669</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2014-04-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.482270</td>\n",
       "      <td>0.420077</td>\n",
       "      <td>0.357624</td>\n",
       "      <td>0.490649</td>\n",
       "      <td>0.525618</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2014-05-01</td>\n",
       "      <td>ESP.1.1.1.10_1</td>\n",
       "      <td>0.426545</td>\n",
       "      <td>0.488816</td>\n",
       "      <td>0.422413</td>\n",
       "      <td>0.554239</td>\n",
       "      <td>0.524601</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1037988</th>\n",
       "      <td>2024-03-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.219000</td>\n",
       "      <td>0.200083</td>\n",
       "      <td>0.142992</td>\n",
       "      <td>0.260890</td>\n",
       "      <td>0.337588</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1037989</th>\n",
       "      <td>2024-04-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.208000</td>\n",
       "      <td>0.235622</td>\n",
       "      <td>0.178164</td>\n",
       "      <td>0.296882</td>\n",
       "      <td>0.336815</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1037990</th>\n",
       "      <td>2024-05-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.354000</td>\n",
       "      <td>0.345233</td>\n",
       "      <td>0.288684</td>\n",
       "      <td>0.400505</td>\n",
       "      <td>0.336067</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1037991</th>\n",
       "      <td>2024-06-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.442000</td>\n",
       "      <td>0.473369</td>\n",
       "      <td>0.418315</td>\n",
       "      <td>0.531617</td>\n",
       "      <td>0.335295</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1037992</th>\n",
       "      <td>2024-07-01</td>\n",
       "      <td>ESP.9.1.9.9_1</td>\n",
       "      <td>0.513000</td>\n",
       "      <td>0.579231</td>\n",
       "      <td>0.519388</td>\n",
       "      <td>0.639722</td>\n",
       "      <td>0.334549</td>\n",
       "      <td>-1</td>\n",
       "      <td>0.012451</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1037993 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                ds           gid_4         y      yhat  yhat_lower  \\\n",
       "0       2014-01-01  ESP.1.1.1.10_1  0.393110  0.322870    0.250874   \n",
       "1       2014-02-01  ESP.1.1.1.10_1  0.402229  0.326708    0.257482   \n",
       "2       2014-03-01  ESP.1.1.1.10_1  0.315096  0.365371    0.298259   \n",
       "3       2014-04-01  ESP.1.1.1.10_1  0.482270  0.420077    0.357624   \n",
       "4       2014-05-01  ESP.1.1.1.10_1  0.426545  0.488816    0.422413   \n",
       "...            ...             ...       ...       ...         ...   \n",
       "1037988 2024-03-01   ESP.9.1.9.9_1  0.219000  0.200083    0.142992   \n",
       "1037989 2024-04-01   ESP.9.1.9.9_1  0.208000  0.235622    0.178164   \n",
       "1037990 2024-05-01   ESP.9.1.9.9_1  0.354000  0.345233    0.288684   \n",
       "1037991 2024-06-01   ESP.9.1.9.9_1  0.442000  0.473369    0.418315   \n",
       "1037992 2024-07-01   ESP.9.1.9.9_1  0.513000  0.579231    0.519388   \n",
       "\n",
       "         yhat_upper     trend  anomaly  importance  \n",
       "0          0.392607  0.528668        1    0.001281  \n",
       "1          0.394605  0.527618        1    0.018955  \n",
       "2          0.437156  0.526669        0    0.000000  \n",
       "3          0.490649  0.525618        0    0.000000  \n",
       "4          0.554239  0.524601        0    0.000000  \n",
       "...             ...       ...      ...         ...  \n",
       "1037988    0.260890  0.337588        0    0.000000  \n",
       "1037989    0.296882  0.336815        0    0.000000  \n",
       "1037990    0.400505  0.336067        0    0.000000  \n",
       "1037991    0.531617  0.335295        0    0.000000  \n",
       "1037992    0.639722  0.334549       -1    0.012451  \n",
       "\n",
       "[1037993 rows x 9 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2b26195e-a513-46bf-9ea6-4640c4c52bb3",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_462799/1002646056.py:3: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.\n",
      "  ).groupby('gid_4').apply(lambda x: x.iloc[-1])[['y', 'yhat', 'yhat_lower', 'yhat_upper', 'trend', 'anomaly', 'importance', 'ds']]\n"
     ]
    }
   ],
   "source": [
    "current_status_df = result_df.sort_values(\n",
    "    by=['gid_4', 'ds']\n",
    ").groupby('gid_4').apply(lambda x: x.iloc[-1])[['y', 'yhat', 'yhat_lower', 'yhat_upper', 'trend', 'anomaly', 'importance', 'ds']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "08c5a5ab-1d1e-4df3-8da6-b85c97e1f0d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "current_status_df.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e26d1601-5689-42b2-84d5-2549c0efba3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "current_status_df.rename(columns={'ds': 'last_update'}, inplace=True)\n",
    "current_status_df['gid_4'] = current_status_df['gid_4'].astype(str)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fb5e491-3c67-42d4-a78c-a6e4edb8a4ab",
   "metadata": {},
   "source": [
    "### Load shapefiles & save geopackage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d5f01d23-889a-46c2-b5c7-00dbd0a824ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "peninsula_gdf = gpd.read_file('lineas_limite/SHP_ETRS89/recintos_municipales_inspire_peninbal_etrs89')\n",
    "peninsula_gdf = peninsula_gdf.to_crs(epsg=4326)\n",
    "\n",
    "canarias_gdf = gpd.read_file('lineas_limite/SHP_REGCAN95/recintos_municipales_inspire_canarias_regcan95')\n",
    "canarias_gdf = canarias_gdf.to_crs(epsg=4326)\n",
    "\n",
    "spain_gdf = gpd.GeoDataFrame(pd.concat([peninsula_gdf, canarias_gdf], ignore_index=True))\n",
    "spain_gdf['NAMEUNIT'] = spain_gdf['NAMEUNIT'].str.split('/').str[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8c79c72b-b1ee-47ab-bc18-310b2a957683",
   "metadata": {},
   "outputs": [],
   "source": [
    "peninsula_ccaa_gdf = gpd.read_file('lineas_limite/SHP_ETRS89/recintos_autonomicas_inspire_peninbal_etrs89')\n",
    "peninsula_ccaa_gdf = peninsula_ccaa_gdf.to_crs(epsg=4326)\n",
    "\n",
    "canarias_ccaa_gdf = gpd.read_file('lineas_limite/SHP_REGCAN95/recintos_autonomicas_inspire_canarias_regcan95')\n",
    "canarias_ccaa_gdf = canarias_ccaa_gdf.to_crs(epsg=4326)\n",
    "\n",
    "spain_ccaa_gdf = gpd.GeoDataFrame(pd.concat([peninsula_ccaa_gdf, canarias_ccaa_gdf], ignore_index=True))\n",
    "spain_ccaa_gdf['NAMEUNIT'] = spain_ccaa_gdf['NAMEUNIT'].str.split('/').str[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3771a53a-9981-487c-9f14-738b27d8e296",
   "metadata": {},
   "outputs": [],
   "source": [
    "gadm4_gdf = gpd.read_file('gadm41_ESP.gpkg', layer='ADM_ADM_4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1ba0d570-f58a-44f8-85a9-a749093615ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "gadm4_gdf['geometry'] = gadm4_gdf.representative_point()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "04ebdf2f-935f-4553-8d86-87ab1b9bba6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "municipalities_gdf = gpd.sjoin(spain_gdf, gadm4_gdf, how=\"left\")[[\n",
    "    'GID_4', 'NATCODE', 'NAMEUNIT', 'CODNUT2', 'geometry'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b0ce2776-6aeb-4132-99f2-51c1c7a5cd30",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = municipalities_gdf[['GID_4', 'NATCODE', 'NAMEUNIT', 'CODNUT2', 'geometry']].merge(\n",
    "    spain_ccaa_gdf[['NAMEUNIT', 'CODNUT2']].rename(columns={'NAMEUNIT': 'NAMEUNIT_NUT2'}),\n",
    "    on='CODNUT2',\n",
    "    how='inner'\n",
    ")\n",
    "gdf['NATCODE'] = gdf['NATCODE'].astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "fab4cebf-eaf0-4ddc-b4f0-150a1d6f2910",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>GID_4</th>\n",
       "      <th>NATCODE</th>\n",
       "      <th>NAMEUNIT</th>\n",
       "      <th>CODNUT2</th>\n",
       "      <th>geometry</th>\n",
       "      <th>NAMEUNIT_NUT2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ESP.17.1.5.3_1</td>\n",
       "      <td>34033333022</td>\n",
       "      <td>Degaña</td>\n",
       "      <td>ES12</td>\n",
       "      <td>MULTIPOLYGON (((-6.6574 42.96745, -6.64737 42....</td>\n",
       "      <td>Principado de Asturias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ESP.17.1.3.4_1</td>\n",
       "      <td>34033333023</td>\n",
       "      <td>El Franco</td>\n",
       "      <td>ES12</td>\n",
       "      <td>MULTIPOLYGON (((-6.87709 43.56358, -6.87705 43...</td>\n",
       "      <td>Principado de Asturias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ESP.17.1.4.2_1</td>\n",
       "      <td>34033333024</td>\n",
       "      <td>Gijón</td>\n",
       "      <td>ES12</td>\n",
       "      <td>MULTIPOLYGON (((-5.81929 43.50727, -5.8184 43....</td>\n",
       "      <td>Principado de Asturias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ESP.17.1.1.6_1</td>\n",
       "      <td>34033333025</td>\n",
       "      <td>Gozón</td>\n",
       "      <td>ES12</td>\n",
       "      <td>MULTIPOLYGON (((-5.91545 43.60853, -5.91537 43...</td>\n",
       "      <td>Principado de Asturias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ESP.17.1.8.4_1</td>\n",
       "      <td>34033333026</td>\n",
       "      <td>Grado</td>\n",
       "      <td>ES12</td>\n",
       "      <td>POLYGON ((-6.20021 43.18357, -6.20121 43.1856,...</td>\n",
       "      <td>Principado de Asturias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8335</th>\n",
       "      <td>ESP.14.2.1.3_1</td>\n",
       "      <td>34053838003</td>\n",
       "      <td>Alajeró</td>\n",
       "      <td>ES70</td>\n",
       "      <td>MULTIPOLYGON (((-17.22374 28.0254, -17.22373 2...</td>\n",
       "      <td>Canarias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8336</th>\n",
       "      <td>ESP.14.2.1.4_1</td>\n",
       "      <td>34053838004</td>\n",
       "      <td>Arafo</td>\n",
       "      <td>ES70</td>\n",
       "      <td>POLYGON ((-16.48414 28.33504, -16.48377 28.336...</td>\n",
       "      <td>Canarias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8337</th>\n",
       "      <td>ESP.14.2.1.5_1</td>\n",
       "      <td>34053838005</td>\n",
       "      <td>Arico</td>\n",
       "      <td>ES70</td>\n",
       "      <td>MULTIPOLYGON (((-16.47283 28.10376, -16.47287 ...</td>\n",
       "      <td>Canarias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8338</th>\n",
       "      <td>ESP.14.2.1.6_1</td>\n",
       "      <td>34053838006</td>\n",
       "      <td>Arona</td>\n",
       "      <td>ES70</td>\n",
       "      <td>MULTIPOLYGON (((-16.6964 28.00131, -16.69639 2...</td>\n",
       "      <td>Canarias</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8339</th>\n",
       "      <td>ESP.14.2.1.7_1</td>\n",
       "      <td>34053838007</td>\n",
       "      <td>Barlovento</td>\n",
       "      <td>ES70</td>\n",
       "      <td>MULTIPOLYGON (((-17.8601 28.76452, -17.85895 2...</td>\n",
       "      <td>Canarias</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>8340 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "               GID_4      NATCODE    NAMEUNIT CODNUT2  \\\n",
       "0     ESP.17.1.5.3_1  34033333022      Degaña    ES12   \n",
       "1     ESP.17.1.3.4_1  34033333023   El Franco    ES12   \n",
       "2     ESP.17.1.4.2_1  34033333024       Gijón    ES12   \n",
       "3     ESP.17.1.1.6_1  34033333025       Gozón    ES12   \n",
       "4     ESP.17.1.8.4_1  34033333026       Grado    ES12   \n",
       "...              ...          ...         ...     ...   \n",
       "8335  ESP.14.2.1.3_1  34053838003     Alajeró    ES70   \n",
       "8336  ESP.14.2.1.4_1  34053838004       Arafo    ES70   \n",
       "8337  ESP.14.2.1.5_1  34053838005       Arico    ES70   \n",
       "8338  ESP.14.2.1.6_1  34053838006       Arona    ES70   \n",
       "8339  ESP.14.2.1.7_1  34053838007  Barlovento    ES70   \n",
       "\n",
       "                                               geometry  \\\n",
       "0     MULTIPOLYGON (((-6.6574 42.96745, -6.64737 42....   \n",
       "1     MULTIPOLYGON (((-6.87709 43.56358, -6.87705 43...   \n",
       "2     MULTIPOLYGON (((-5.81929 43.50727, -5.8184 43....   \n",
       "3     MULTIPOLYGON (((-5.91545 43.60853, -5.91537 43...   \n",
       "4     POLYGON ((-6.20021 43.18357, -6.20121 43.1856,...   \n",
       "...                                                 ...   \n",
       "8335  MULTIPOLYGON (((-17.22374 28.0254, -17.22373 2...   \n",
       "8336  POLYGON ((-16.48414 28.33504, -16.48377 28.336...   \n",
       "8337  MULTIPOLYGON (((-16.47283 28.10376, -16.47287 ...   \n",
       "8338  MULTIPOLYGON (((-16.6964 28.00131, -16.69639 2...   \n",
       "8339  MULTIPOLYGON (((-17.8601 28.76452, -17.85895 2...   \n",
       "\n",
       "               NAMEUNIT_NUT2  \n",
       "0     Principado de Asturias  \n",
       "1     Principado de Asturias  \n",
       "2     Principado de Asturias  \n",
       "3     Principado de Asturias  \n",
       "4     Principado de Asturias  \n",
       "...                      ...  \n",
       "8335                Canarias  \n",
       "8336                Canarias  \n",
       "8337                Canarias  \n",
       "8338                Canarias  \n",
       "8339                Canarias  \n",
       "\n",
       "[8340 rows x 6 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "9dd2efac-9b82-411e-8b88-071d3e97e54d",
   "metadata": {},
   "outputs": [],
   "source": [
    "current_gdf = gdf.merge(current_status_df.rename(columns={'gid_4': 'GID_4'}), on='GID_4')\n",
    "current_gdf.set_index('NATCODE', inplace=True)\n",
    "current_gdf.drop(columns=['GID_4'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "7e26897f-3782-4d08-aae8-82a949e3b376",
   "metadata": {},
   "outputs": [],
   "source": [
    "historic_gdf = gpd.GeoDataFrame(\n",
    "    result_df.merge(\n",
    "        gdf[['NATCODE', 'GID_4']].rename(columns={'GID_4': 'gid_4'}),\n",
    "        on='gid_4',\n",
    "        how='inner'\n",
    "    ).drop(columns=['gid_4']),\n",
    "    geometry=None\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "f5f9c416-4981-4741-81b7-5234fa5bbfc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpk_path = 'output_bites.gpkg'\n",
    "# Save the GeoPandas DataFrame (geometries)\n",
    "current_gdf.to_file(gpk_path, layer='geometries', driver=\"GPKG\")\n",
    "historic_gdf.to_file(gpk_path, layer='histories', driver=\"GPKG\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f9a0e5b1-564e-4822-a7b3-aaa4751a50df",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpd.GeoDataFrame(yearly_seasonality_df.merge(\n",
    "        gdf[['NATCODE', 'GID_4']].rename(columns={'GID_4': 'gid_4'}),\n",
    "        on='gid_4',\n",
    "        how='inner'\n",
    "    ).drop(columns=['gid_4']), geometry=None).to_file(gpk_path, layer='seasonality', driver=\"GPKG\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa74897a-0ddf-467c-aad3-799c9e8c30ac",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}