Compute_MP.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "This code did data cleaning, computed KNDVI, and matrix profile.\n",
        "Users are advised to pay attention to all sections of the code.\n",
        "\n",
        "Last update Thursday Oct 31 2024\n",
        "@author: Pius N.Nwachukwu"
      ],
      "metadata": {
        "id": "P9Muf4XodA2K"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TGncVA4sUv_d"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install jupyter_bokeh\n"
      ],
      "metadata": {
        "id": "qIgnOYn-D94e"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install chord"
      ],
      "metadata": {
        "id": "hvtqFjfasYjC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "R7eKAHWJT27z"
      },
      "outputs": [],
      "source": [
        "pip install holoviews bokeh\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "l1zg5mBnXFUB"
      },
      "outputs": [],
      "source": [
        "!pip install stumpy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "w9sgu1hHjxIB"
      },
      "outputs": [],
      "source": [
        "import os\n",
        "import shutil\n",
        "import ee\n",
        "import geemap\n",
        "import pandas as pd\n",
        "from google.colab import drive\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "id": "Zwl9mjS3jq5g",
        "outputId": "5243250e-874f-46be-9e54-be4fab257c9b"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "            <style>\n",
              "                .geemap-dark {\n",
              "                    --jp-widgets-color: white;\n",
              "                    --jp-widgets-label-color: white;\n",
              "                    --jp-ui-font-color1: white;\n",
              "                    --jp-layout-color2: #454545;\n",
              "                    background-color: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-dark .jupyter-button {\n",
              "                    --jp-layout-color3: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-colab {\n",
              "                    background-color: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "\n",
              "                .geemap-colab .jupyter-button {\n",
              "                    --jp-layout-color3: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "            </style>\n",
              "            "
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ],
      "source": [
        "# Authenticate and initialize the Earth Engine module.\n",
        "ee.Authenticate()\n",
        "ee.Initialize(project='ee-')# Please provide a project name\n",
        "\n",
        "# Mount Google Drive\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 124
        },
        "id": "BQ-hwGjp4Euk",
        "outputId": "10d437b9-ea6c-4cc7-8a8a-5d99cd575fb9"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "            <style>\n",
              "                .geemap-dark {\n",
              "                    --jp-widgets-color: white;\n",
              "                    --jp-widgets-label-color: white;\n",
              "                    --jp-ui-font-color1: white;\n",
              "                    --jp-layout-color2: #454545;\n",
              "                    background-color: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-dark .jupyter-button {\n",
              "                    --jp-layout-color3: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-colab {\n",
              "                    background-color: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "\n",
              "                .geemap-colab .jupyter-button {\n",
              "                    --jp-layout-color3: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "            </style>\n",
              "            "
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "        ROI        date  percentage_valid  latitude  longitude  NDVI\n",
            "0  S_Americ  2000-02-01          29.70645 -6.609515 -58.598977  1355\n",
            "1  S_Americ  2000-02-01          29.70645 -6.609515 -58.508443  4857\n",
            "2  S_Americ  2000-02-01          29.70645 -6.654481 -58.694867  3963\n",
            "3  S_Americ  2000-02-01          29.70645 -6.654481 -58.649596  1773\n",
            "4  S_Americ  2000-02-01          29.70645 -6.654481 -58.604325  2665\n"
          ]
        }
      ],
      "source": [
        "# Calling up the Saved Time Series NDVI Values\n",
        "# csv_file_path = '/content/drive/My Drive/FoRes/gq_ndvi_timeseries_data.csv' # Get the NDVIs\n",
        "csv_file_path = '/content/drive/My Drive/FoRes/ndvi_pixelz_data.csv' # Get the NDVIs in 5000 SR Pixels\n",
        "\n",
        "\n",
        "# Read the CSV file into a pandas DataFrame\n",
        "ndvi_df = pd.read_csv(csv_file_path)\n",
        "print(ndvi_df.head())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 17
        },
        "id": "tAL0ZvF9I3VJ",
        "outputId": "e2980a7d-35bc-4c94-ba6d-2392f7f7d473"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "            <style>\n",
              "                .geemap-dark {\n",
              "                    --jp-widgets-color: white;\n",
              "                    --jp-widgets-label-color: white;\n",
              "                    --jp-ui-font-color1: white;\n",
              "                    --jp-layout-color2: #454545;\n",
              "                    background-color: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-dark .jupyter-button {\n",
              "                    --jp-layout-color3: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-colab {\n",
              "                    background-color: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "\n",
              "                .geemap-colab .jupyter-button {\n",
              "                    --jp-layout-color3: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "            </style>\n",
              "            "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "ndvi_df.set_index('date', inplace=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 17
        },
        "id": "rsA9II1Lt8iC",
        "outputId": "5c8d9fc6-e0e3-4322-b9aa-6eb66b4cb58e"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "            <style>\n",
              "                .geemap-dark {\n",
              "                    --jp-widgets-color: white;\n",
              "                    --jp-widgets-label-color: white;\n",
              "                    --jp-ui-font-color1: white;\n",
              "                    --jp-layout-color2: #454545;\n",
              "                    background-color: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-dark .jupyter-button {\n",
              "                    --jp-layout-color3: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-colab {\n",
              "                    background-color: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "\n",
              "                .geemap-colab .jupyter-button {\n",
              "                    --jp-layout-color3: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "            </style>\n",
              "            "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# Converting the NDVI to Reflectance\n",
        "ndvi_df[\"NDVI_2\"] = ndvi_df[\"NDVI\"] / 10000"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 17
        },
        "id": "MQ7vL2-gqxdG",
        "outputId": "021aec8b-a891-4e6e-c19f-629ef2bb646c"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "            <style>\n",
              "                .geemap-dark {\n",
              "                    --jp-widgets-color: white;\n",
              "                    --jp-widgets-label-color: white;\n",
              "                    --jp-ui-font-color1: white;\n",
              "                    --jp-layout-color2: #454545;\n",
              "                    background-color: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-dark .jupyter-button {\n",
              "                    --jp-layout-color3: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-colab {\n",
              "                    background-color: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "\n",
              "                .geemap-colab .jupyter-button {\n",
              "                    --jp-layout-color3: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "            </style>\n",
              "            "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# # Calculate KNDVI\n",
        "ndvi_df['KNDVI'] = (ndvi_df['NDVI_2'] + 1) / 2\n",
        "\n",
        "# print(ndvi_df)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "id": "aHDciD-vO1Uk",
        "outputId": "a8fab9ee-7247-438b-8f6d-f7849e025d12"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "\n",
              "            <style>\n",
              "                .geemap-dark {\n",
              "                    --jp-widgets-color: white;\n",
              "                    --jp-widgets-label-color: white;\n",
              "                    --jp-ui-font-color1: white;\n",
              "                    --jp-layout-color2: #454545;\n",
              "                    background-color: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-dark .jupyter-button {\n",
              "                    --jp-layout-color3: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-colab {\n",
              "                    background-color: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "\n",
              "                .geemap-colab .jupyter-button {\n",
              "                    --jp-layout-color3: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "            </style>\n",
              "            "
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "['S_Americ' 'N_Americ' 'Afric' 'Europ' 'Asia' 'Aust_Ocean']\n"
          ]
        }
      ],
      "source": [
        "# ndvi_df['ROI']\n",
        "\n",
        "# Print all the names of the ROIs\n",
        "print(ndvi_df['ROI'].unique())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 17
        },
        "id": "OHe3BRWQ-Y9S",
        "outputId": "11713940-fa8e-49ec-e1c2-d999eaf003e4"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "            <style>\n",
              "                .geemap-dark {\n",
              "                    --jp-widgets-color: white;\n",
              "                    --jp-widgets-label-color: white;\n",
              "                    --jp-ui-font-color1: white;\n",
              "                    --jp-layout-color2: #454545;\n",
              "                    background-color: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-dark .jupyter-button {\n",
              "                    --jp-layout-color3: #383838;\n",
              "                }\n",
              "\n",
              "                .geemap-colab {\n",
              "                    background-color: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "\n",
              "                .geemap-colab .jupyter-button {\n",
              "                    --jp-layout-color3: var(--colab-primary-surface-color, white);\n",
              "                }\n",
              "            </style>\n",
              "            "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# Filter the DataFrame to get rows where the 'Region' column is correspnd to a continent eg 'Africa'\n",
        "af = ndvi_df[ndvi_df[\"ROI\"] == \"Afric\"]\n",
        "asia = ndvi_df[ndvi_df[\"ROI\"] == \"Asia\"]\n",
        "eu = ndvi_df[ndvi_df[\"ROI\"] == \"Europ\"]\n",
        "na = ndvi_df[ndvi_df[\"ROI\"] == \"N_Americ\"]\n",
        "sa = ndvi_df[ndvi_df[\"ROI\"] == \"S_Americ\"]\n",
        "oc = ndvi_df[ndvi_df[\"ROI\"] == \"Aust_Ocean\"]\n",
        "# print(eu)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DqCIR1B-koG6"
      },
      "outputs": [],
      "source": [
        "# Define the ROI name you want to filter\n",
        "selected_roi = 'Europ'  # Replace current continent/ROI\n",
        "\n",
        "# Filter the DataFrame for pixels within the selected ROI\n",
        "df_filtered = ndvi_df[ndvi_df['ROI'] == selected_roi]\n",
        "\n",
        "# Display the first few rows of the filtered DataFrame\n",
        "print(df_filtered.head())\n",
        "\n",
        "# # Optionally, save the filtered DataFrame to a CSV file\n",
        "# df_filtered.to_csv(f'{selected_roi}_ndvi_pixel_data.csv', index=False)\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Rounding the coordinates to group nearby coordinates to effectively clusters locations within a certain decimal precision\n",
        "#\n",
        "df_filtered['latitude'] = df_filtered_copy['latitude'].round(2)\n",
        "df_filtered['longitude'] = df_filtered_copy['longitude'].round(2)\n",
        "df_filtered\n",
        "\n"
      ],
      "metadata": {
        "id": "89CELKW8Wz-O"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "# Group the data by latitude and longitude to count the unique pixels\n",
        "unique_pixels = df_filtered.groupby(['latitude', 'longitude']).size()\n",
        "\n",
        "# Print the number of unique pixels\n",
        "print(f\"Number of unique pixels: {len(unique_pixels)}\")\n",
        "\n",
        "# If you want to see a sample of the pixel data\n",
        "print(unique_pixels.head())"
      ],
      "metadata": {
        "id": "4KXysn9TXioN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Confirm Pixel Characteristics\n",
        "\n",
        "# Define the range for valid KNDVI values\n",
        "valid_ndvi_min = -1\n",
        "valid_ndvi_max = 1\n",
        "\n",
        "# Filter the dataset for valid NDVI values\n",
        "valid_ndvi_data = df_filtered[(df_filtered['NDVI_2'] >= valid_ndvi_min) & (df_filtered['NDVI_2'] <= valid_ndvi_max)]\n",
        "\n",
        "# Group by latitude and longitude to identify unique pixels with valid data\n",
        "valid_pixels = valid_ndvi_data.groupby(['latitude', 'longitude']).size()\n",
        "\n",
        "# Calculate the percentage of valid pixels\n",
        "percentage_valid_pixels = (len(valid_pixels) / 1277) * 100  # Assuming 5460 is the total number of unique pixels\n",
        "\n",
        "print(f\"Percentage of pixels with valid NDVI data: {percentage_valid_pixels:.2f}%\")\n"
      ],
      "metadata": {
        "id": "H6JzTL4aXuJ5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Reshapes the data so that each row represents a unique pixel (defined by latitude and longitude)\n",
        "\n",
        "# Pivot the data to create a time series for each pixel\n",
        "ndvi_pivot = df_filtered.pivot_table(index=['latitude', 'longitude'], columns='date', values='KNDVI')\n",
        "\n",
        "# Display the reshaped DataFrame\n",
        "print(ndvi_pivot.head())\n",
        "\n",
        "# Now, check the length of the time series for one pixel\n",
        "print(f\"Number of time points: {ndvi_pivot.shape[1]}\")"
      ],
      "metadata": {
        "id": "bLfsKwa3YS6W"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# df_filtered.index"
      ],
      "metadata": {
        "id": "qaZPPOJaZjpx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Compute MP with Stumpy (With zero Exclusion)\n",
        "\n",
        "import pandas as pd\n",
        "import stumpy\n",
        "\n",
        "# Fill missing values using interpolation or fill with mean (as needed)\n",
        "ndvi_pivot_filled = ndvi_pivot.interpolate(method='linear', axis=1).fillna(ndvi_pivot.mean(axis=1))\n",
        "\n",
        "# Initialize a list to store results\n",
        "results = []\n",
        "\n",
        "# Loop through each pixel (each row in the pivot table)\n",
        "for index, row in ndvi_pivot_filled.iterrows():\n",
        "    # Extract the NDVI time series for this pixel\n",
        "    ndvi_series = row.values\n",
        "\n",
        "    # Define a reasonable window size\n",
        "    window_size = 24\n",
        "\n",
        "    # Compute the matrix profile for the time series\n",
        "    mp = stumpy.stump(ndvi_series, m=window_size)\n",
        "\n",
        "    # Identify the motif (minimum matrix profile value) and its nearest neighbor\n",
        "    motif_idx = mp[:, 0].argmin()\n",
        "    distance = mp[motif_idx, 0]  # Matrix profile distance\n",
        "    neighbor_idx = int(mp[motif_idx, 1])\n",
        "\n",
        "    # Only proceed if the distance is greater than zero (exclude exact matches)\n",
        "    if distance > 0:\n",
        "        # Identify the discord (maximum matrix profile value)\n",
        "        discord_idx = mp[:, 0].argmax()\n",
        "\n",
        "        # Extract the dates corresponding to the motif, neighbor, and discord\n",
        "        motif_date = ndvi_pivot_filled.columns[motif_idx]\n",
        "        neighbor_date = ndvi_pivot_filled.columns[neighbor_idx]\n",
        "        discord_date = ndvi_pivot_filled.columns[discord_idx]\n",
        "\n",
        "        # Append results for this pixel to the list\n",
        "        results.append({\n",
        "            'latitude': index[0],\n",
        "            'longitude': index[1],\n",
        "            'motif_date': motif_date,\n",
        "            'neighbor_date': neighbor_date,\n",
        "            'discord_date': discord_date,\n",
        "            'motif_value': ndvi_series[motif_idx],\n",
        "            'neighbor_value': ndvi_series[neighbor_idx],\n",
        "            'discord_value': ndvi_series[discord_idx],\n",
        "            'distance': distance  # Include the distance in the results\n",
        "        })\n",
        "\n",
        "# Convert the results list to a DataFrame\n",
        "results_df = pd.DataFrame(results)\n",
        "\n",
        "# # Save the results to a CSV file\n",
        "# results_df.to_csv('ndvi_motifs_discords_excluding_zero_distance.csv', index=False)\n",
        "\n",
        "# # Display the first few rows of the results\n",
        "# print(results_df.head())\n"
      ],
      "metadata": {
        "id": "vDc5bqOmZxCT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Plotting  data point (month).\n",
        "\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.dates as mdates\n",
        "import pandas as pd\n",
        "\n",
        "# Ensure that the date columns are in datetime format\n",
        "ndvi_pivot_filled.columns = pd.to_datetime(ndvi_pivot_filled.columns)\n",
        "\n",
        "# Select a pixel (e.g., the first pixel in the pivot table)\n",
        "selected_pixel_series = ndvi_pivot_filled.iloc[3].values  # NDVI time series for the selected pixel\n",
        "selected_pixel_index = ndvi_pivot_filled.index[3]  # Get the corresponding latitude and longitude\n",
        "\n",
        "# Define a reasonable window size\n",
        "window_size = 24\n",
        "\n",
        "# Compute the matrix profile for the selected pixel\n",
        "mp = stumpy.stump(selected_pixel_series, m=window_size)\n",
        "\n",
        "# Identify the motif, neighbor, and discord indices\n",
        "motif_idx = mp[:, 0].argmin()\n",
        "neighbor_idx = int(mp[motif_idx, 1])\n",
        "discord_idx = mp[:, 0].argmax()\n",
        "\n",
        "# Plot the NDVI time series\n",
        "plt.figure(figsize=(7, 4))\n",
        "\n",
        "# NDVI Time Series\n",
        "plt.subplot(2, 1, 1)\n",
        "plt.plot(ndvi_pivot_filled.columns, selected_pixel_series, label='NDVI Time Series')\n",
        "\n",
        "# Highlight motif, neighbor, and discord\n",
        "plt.plot(ndvi_pivot_filled.columns[motif_idx:motif_idx+window_size], selected_pixel_series[motif_idx:motif_idx+window_size], color='green', label='Motif')\n",
        "plt.plot(ndvi_pivot_filled.columns[neighbor_idx:neighbor_idx+window_size], selected_pixel_series[neighbor_idx:neighbor_idx+window_size], color='orange', label='Neighbor')\n",
        "plt.plot(ndvi_pivot_filled.columns[discord_idx:discord_idx+window_size], selected_pixel_series[discord_idx:discord_idx+window_size], color='red', label='Discord')\n",
        "\n",
        "plt.title(f'KNDVI Time Series for Nth America at Latitude {selected_pixel_index[0]}, Longitude {selected_pixel_index[1]}')\n",
        "plt.xlabel('Date')\n",
        "plt.ylabel('KNDVI')\n",
        "plt.xticks(rotation=45)\n",
        "\n",
        "# Set major ticks to the first month of each year\n",
        "plt.gca().xaxis.set_major_locator(mdates.YearLocator())\n",
        "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))\n",
        "\n",
        "# Matrix Profile\n",
        "plt.subplot(2, 1, 2)\n",
        "plt.plot(ndvi_pivot_filled.columns[window_size-1:], mp[:, 0], label='Matrix Profile', color='orange')\n",
        "\n",
        "# Highlight motif, neighbor, and discord in the matrix profile\n",
        "plt.scatter(ndvi_pivot_filled.columns[motif_idx + window_size - 1], mp[motif_idx, 0], color='green', label='Motif')\n",
        "plt.scatter(ndvi_pivot_filled.columns[neighbor_idx + window_size - 1], mp[neighbor_idx, 0], color='orange', label='Neighbor')\n",
        "plt.scatter(ndvi_pivot_filled.columns[discord_idx + window_size - 1], mp[discord_idx, 0], color='black', label='Discord')\n",
        "\n",
        "plt.title('Matrix Profile')\n",
        "plt.xlabel('Date')\n",
        "plt.ylabel('Matrix Profile Value')\n",
        "plt.xticks(rotation=45)\n",
        "\n",
        "# Set major ticks to the first month of each year\n",
        "plt.gca().xaxis.set_major_locator(mdates.YearLocator())\n",
        "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))\n",
        "\n",
        "# Adjust the layout and place the legend below the plots\n",
        "plt.tight_layout(rect=[0, 0, 1, 1.2])\n",
        "plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.70), ncol=4)\n",
        "# 0.5, 0.1\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "hHEGj75XaD5W"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Compute MP and Reurn the Year Occurence of Motif and Discord\n",
        "\n",
        "import pandas as pd\n",
        "from collections import defaultdict\n",
        "import stumpy\n",
        "\n",
        "# Create dictionaries to store motif and discord counts by year\n",
        "motif_counts_by_year = defaultdict(int)\n",
        "discord_counts_by_year = defaultdict(int)\n",
        "\n",
        "# Loop through the motif-neighbor pairs and count occurrences by year\n",
        "for index, row in ndvi_pivot_filled.iterrows():\n",
        "    ndvi_series = row.values\n",
        "    dates = row.index  # Assuming the index of the row contains the dates\n",
        "\n",
        "    window_size = 24\n",
        "    mp = stumpy.stump(ndvi_series, m=window_size)\n",
        "\n",
        "    # Identify the motif (minimum matrix profile value) and its nearest neighbor\n",
        "    motif_idx = mp[:, 0].argmin()\n",
        "    motif_distance = mp[motif_idx, 0]  # Matrix profile distance for motif\n",
        "\n",
        "    # Only count motifs where the distance is greater than zero\n",
        "    if motif_distance > 0:\n",
        "        motif_year = pd.to_datetime(dates[motif_idx]).year  # Extract year from the date\n",
        "        motif_counts_by_year[motif_year] += 1  # Increase count for this year\n",
        "\n",
        "    # Identify the discord (maximum matrix profile value)\n",
        "    discord_idx = mp[:, 0].argmax()\n",
        "    discord_distance = mp[discord_idx, 0]  # Matrix profile distance for discord\n",
        "\n",
        "    # Only count discords where the distance is greater than zero\n",
        "    if discord_distance > 0:\n",
        "        discord_year = pd.to_datetime(dates[discord_idx]).year  # Extract year from the date\n",
        "        discord_counts_by_year[discord_year] += 1  # Increase discord count for this year\n",
        "\n",
        "# Convert motif counts to a DataFrame for easier manipulation\n",
        "motif_counts_by_year_df = pd.DataFrame(list(motif_counts_by_year.items()), columns=['Year', 'Motif Count'])\n",
        "\n",
        "# Calculate the total motifs\n",
        "total_motifs = motif_counts_by_year_df['Motif Count'].sum()\n",
        "\n",
        "# Calculate the percentage of motifs for each year\n",
        "motif_counts_by_year_df['Percentage'] = (motif_counts_by_year_df['Motif Count'] / total_motifs) * 100\n",
        "\n",
        "# Sort by year\n",
        "motif_counts_by_year_df = motif_counts_by_year_df.sort_values('Year')\n",
        "\n",
        "# Convert discord counts to a DataFrame\n",
        "discord_counts_by_year_df = pd.DataFrame(list(discord_counts_by_year.items()), columns=['Year', 'Discord Count'])\n",
        "\n",
        "# Calculate the total discords\n",
        "total_discords = discord_counts_by_year_df['Discord Count'].sum()\n",
        "\n",
        "# Calculate the percentage of discords for each year\n",
        "discord_counts_by_year_df['Percentage'] = (discord_counts_by_year_df['Discord Count'] / total_discords) * 100\n",
        "\n",
        "# Sort by year\n",
        "discord_counts_by_year_df = discord_counts_by_year_df.sort_values('Year')\n",
        "\n",
        "# Display the motif and discord results\n",
        "print(\"Motif Counts by Year:\")\n",
        "print(motif_counts_by_year_df)\n",
        "\n",
        "print(\"\\nDiscord Counts by Year:\")\n",
        "print(discord_counts_by_year_df)\n"
      ],
      "metadata": {
        "id": "reQGGPcCaztO"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# V1 Compute MP and return Only Motif\n",
        "import pandas as pd\n",
        "from collections import defaultdict\n",
        "import stumpy\n",
        "\n",
        "# Assuming your data has a 'date' column to group by year or period\n",
        "# Create a dictionary to store motif occurrences by time, month, and pixel\n",
        "motif_occurrences_by_month = defaultdict(list)\n",
        "\n",
        "# Loop through the motif-neighbor pairs and count occurrences by month\n",
        "for pixel_index, row in ndvi_pivot_filled.iterrows():\n",
        "    ndvi_series = row.values\n",
        "    dates = row.index  # Assuming the index of the row contains the dates\n",
        "\n",
        "    window_size = 24\n",
        "    mp = stumpy.stump(ndvi_series, m=window_size)\n",
        "\n",
        "    motif_idx = mp[:, 0].argmin()\n",
        "    distance = mp[motif_idx, 0]  # Matrix profile distance\n",
        "\n",
        "    # Exclude motifs where the distance is exactly zero\n",
        "    if distance > 0:\n",
        "        motif_sequence = tuple(ndvi_series[motif_idx:motif_idx + window_size])\n",
        "        motif_time = dates[motif_idx]  # Get the time of the motif\n",
        "\n",
        "        # Convert motif_time to datetime object if it's a string\n",
        "        if isinstance(motif_time, str):\n",
        "            motif_time = pd.to_datetime(motif_time)\n",
        "\n",
        "        motif_month = motif_time.month  # Extract the month\n",
        "        motif_year = motif_time.year  # Extract the year\n",
        "\n",
        "        # Record the motif, the pixel index where it was found, the month, and the year\n",
        "        motif_occurrences_by_month[(motif_year, motif_month, motif_time)].append((pixel_index, motif_sequence))\n",
        "\n",
        "# Filter motifs that occur in the same month and year across multiple pixels\n",
        "common_motifs_by_month = {\n",
        "    (year, month, time): pixels for (year, month, time), pixels in motif_occurrences_by_month.items() if len(pixels) > 1\n",
        "}\n",
        "\n",
        "# Convert the common motifs to a DataFrame for easier manipulation\n",
        "common_motifs_list = []\n",
        "for (year, month, time), pixels in common_motifs_by_month.items():\n",
        "    for pixel_index, motif_sequence in pixels:\n",
        "        common_motifs_list.append({\n",
        "            'Year': year,\n",
        "            'Month': month,\n",
        "            'Time': time,\n",
        "            'Pixel Index': pixel_index,\n",
        "            'Motif Sequence': motif_sequence\n",
        "        })\n",
        "\n",
        "common_motifs_df = pd.DataFrame(common_motifs_list)\n",
        "\n",
        "# Sort the DataFrame by year and month\n",
        "common_motifs_df = common_motifs_df.sort_values(by=['Year', 'Month', 'Time'])\n",
        "\n",
        "# Display the results\n",
        "print(\"Common motifs observed at the same time across multiple pixels, sorted by year and month:\")\n",
        "print(common_motifs_df)\n"
      ],
      "metadata": {
        "id": "4Cr_C3C0bbiQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# V2 Compute MP and return Motif & 1 Nearest Neighbours\n",
        "import pandas as pd\n",
        "from collections import defaultdict\n",
        "import stumpy\n",
        "\n",
        "# Dictionary to store motif occurrences with nearest neighbors by time, month, and pixel\n",
        "motif_occurrences_by_month = defaultdict(list)\n",
        "\n",
        "# Loop through the motif-neighbor pairs and count occurrences by month\n",
        "for pixel_index, row in ndvi_pivot_filled.iterrows():\n",
        "    ndvi_series = row.values\n",
        "    dates = row.index  # Assuming the index of the row contains the dates\n",
        "\n",
        "    window_size = 24\n",
        "    mp = stumpy.stump(ndvi_series, m=window_size)\n",
        "\n",
        "    # Identify the motif and its nearest neighbors\n",
        "    motif_idx = mp[:, 0].argmin()  # Index of motif with the smallest distance\n",
        "    distance = mp[motif_idx, 0]  # Matrix profile distance of the motif\n",
        "\n",
        "    # Exclude motifs where the distance is exactly zero\n",
        "    if distance > 0:\n",
        "        motif_sequence = tuple(ndvi_series[motif_idx:motif_idx + window_size])\n",
        "        motif_time = dates[motif_idx]  # Get the time of the motif\n",
        "\n",
        "        # Convert motif_time to datetime if it's a string\n",
        "        if isinstance(motif_time, str):\n",
        "            motif_time = pd.to_datetime(motif_time)\n",
        "\n",
        "        motif_month = motif_time.month  # Extract the month\n",
        "        motif_year = motif_time.year  # Extract the year\n",
        "\n",
        "        # Capture the indices and sequences for the nearest neighbors\n",
        "        neighbors = []\n",
        "        for i in range(1, 2):  # Assuming you want the two nearest neighbors\n",
        "            neighbor_idx = mp[motif_idx, i]\n",
        "            neighbor_sequence = tuple(ndvi_series[int(neighbor_idx):int(neighbor_idx) + window_size])\n",
        "            neighbor_time = dates[int(neighbor_idx)]\n",
        "            neighbors.append({\n",
        "                'Neighbor Sequence': neighbor_sequence,\n",
        "                'Neighbor Time': neighbor_time\n",
        "            })\n",
        "\n",
        "        # Record the motif, pixel index, month, year, and nearest neighbors\n",
        "        motif_occurrences_by_month[(motif_year, motif_month, motif_time)].append({\n",
        "            'Pixel Index': pixel_index,\n",
        "            'Motif Sequence': motif_sequence,\n",
        "            'Neighbors': neighbors\n",
        "        })\n",
        "\n",
        "# Filter motifs with nearest neighbors that occur in the same month and year across multiple pixels\n",
        "common_motifs_by_month = {\n",
        "    (year, month, time): pixels for (year, month, time), pixels in motif_occurrences_by_month.items() if len(pixels) > 1\n",
        "}\n",
        "\n",
        "# Convert the common motifs with neighbors to a DataFrame\n",
        "common_motifs_list = []\n",
        "for (year, month, time), pixels in common_motifs_by_month.items():\n",
        "    for data in pixels:\n",
        "        pixel_index = data['Pixel Index']\n",
        "        motif_sequence = data['Motif Sequence']\n",
        "        neighbors = data['Neighbors']\n",
        "\n",
        "        # Add each neighbor to the DataFrame\n",
        "        for neighbor in neighbors:\n",
        "            common_motifs_list.append({\n",
        "                'Year': year,\n",
        "                'Month': month,\n",
        "                'Time': time,\n",
        "                'Pixel Index': pixel_index,\n",
        "                'Motif Sequence': motif_sequence,\n",
        "                'Neighbor Sequence': neighbor['Neighbor Sequence'],\n",
        "                'Neighbor Time': neighbor['Neighbor Time']\n",
        "            })\n",
        "\n",
        "# Create DataFrame and sort by year and month\n",
        "common_motifs_df = pd.DataFrame(common_motifs_list)\n",
        "common_motifs_df = common_motifs_df.sort_values(by=['Year', 'Month', 'Time'])\n",
        "\n",
        "# Display the results\n",
        "print(\"Common motifs with nearest neighbors observed at the same time across multiple pixels, sorted by year and month:\")\n",
        "print(common_motifs_df)\n"
      ],
      "metadata": {
        "id": "trg7JlhvbleT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# V3 Compute MP and return Motif & 2 Nearest Neighbours\n",
        "import pandas as pd\n",
        "from collections import defaultdict\n",
        "import stumpy\n",
        "\n",
        "# Dictionary to store motif occurrences with nearest neighbors by time, month, and pixel\n",
        "motif_occurrences_by_month = defaultdict(list)\n",
        "\n",
        "# Loop through the motif-neighbor pairs and count occurrences by month\n",
        "for pixel_index, row in ndvi_pivot_filled.iterrows():\n",
        "    ndvi_series = row.values\n",
        "    dates = row.index  # Assuming the index of the row contains the dates\n",
        "\n",
        "    window_size = 24\n",
        "    mp = stumpy.stump(ndvi_series, m=window_size)\n",
        "\n",
        "    # Identify the motif and its nearest neighbors\n",
        "    motif_idx = mp[:, 0].argmin()  # Index of motif with the smallest distance\n",
        "    distance = mp[motif_idx, 0]  # Matrix profile distance of the motif\n",
        "\n",
        "    # Exclude motifs where the distance is exactly zero\n",
        "    if distance > 0:\n",
        "        motif_sequence = tuple(ndvi_series[motif_idx:motif_idx + window_size])\n",
        "        motif_time = dates[motif_idx]  # Get the time of the motif\n",
        "\n",
        "        # Convert motif_time to datetime if it's a string\n",
        "        if isinstance(motif_time, str):\n",
        "            motif_time = pd.to_datetime(motif_time)\n",
        "\n",
        "        motif_month = motif_time.month  # Extract the month\n",
        "        motif_year = motif_time.year  # Extract the year\n",
        "\n",
        "        # Capture the indices and sequences for the nearest neighbors\n",
        "        neighbors = []\n",
        "        for i in range(1, 3):  # Assuming you want the two nearest neighbors\n",
        "            neighbor_idx = mp[motif_idx, i]\n",
        "            neighbor_sequence = tuple(ndvi_series[int(neighbor_idx):int(neighbor_idx) + window_size])\n",
        "            neighbor_time = dates[int(neighbor_idx)]\n",
        "            neighbors.append({\n",
        "                'Neighbor Sequence': neighbor_sequence,\n",
        "                'Neighbor Time': neighbor_time\n",
        "            })\n",
        "\n",
        "        # Record the motif, pixel index, month, year, and nearest neighbors\n",
        "        motif_occurrences_by_month[(motif_year, motif_month, motif_time)].append({\n",
        "            'Pixel Index': pixel_index,\n",
        "            'Motif Sequence': motif_sequence,\n",
        "            'Neighbors': neighbors\n",
        "        })\n",
        "\n",
        "# Filter motifs with nearest neighbors that occur in the same month and year across multiple pixels\n",
        "common_motifs_by_month = {\n",
        "    (year, month, time): pixels for (year, month, time), pixels in motif_occurrences_by_month.items() if len(pixels) > 1\n",
        "}\n",
        "\n",
        "# Convert the common motifs with neighbors to a DataFrame\n",
        "common_motifs_list = []\n",
        "for (year, month, time), pixels in common_motifs_by_month.items():\n",
        "    for data in pixels:\n",
        "        pixel_index = data['Pixel Index']\n",
        "        motif_sequence = data['Motif Sequence']\n",
        "        neighbors = data['Neighbors']\n",
        "\n",
        "        # Add each neighbor to the DataFrame\n",
        "        for neighbor in neighbors:\n",
        "            common_motifs_list.append({\n",
        "                'Year': year,\n",
        "                'Month': month,\n",
        "                'Time': time,\n",
        "                'Pixel Index': pixel_index,\n",
        "                'Motif Sequence': motif_sequence,\n",
        "                'Neighbor Sequence': neighbor['Neighbor Sequence'],\n",
        "                'Neighbor Time': neighbor['Neighbor Time']\n",
        "            })\n",
        "\n",
        "# Create DataFrame and sort by year and month\n",
        "common_motifs_df = pd.DataFrame(common_motifs_list)\n",
        "common_motifs_df = common_motifs_df.sort_values(by=['Year', 'Month', 'Time'])\n",
        "\n",
        "# Display the results\n",
        "print(\"Common motifs with nearest neighbors observed at the same time across multiple pixels, sorted by year and month:\")\n",
        "print(common_motifs_df)\n"
      ],
      "metadata": {
        "id": "cQk-_3ombGm4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Plottings"
      ],
      "metadata": {
        "id": "GESg6D1TcYUN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Heatmap\n",
        "\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Filter for valid motifs and neighbor sequences\n",
        "filtered_df = common_motifs_df[(common_motifs_df['Motif Sequence'].notna()) & (common_motifs_df['Neighbor Sequence'] != '()')]\n",
        "\n",
        "# Count the occurrences of motifs and nearest neighbors by Year and Month\n",
        "monthly_counts = filtered_df.groupby(['Year', 'Month']).size().reset_index(name='Total_Count')\n",
        "\n",
        "# Pivot the data for heatmap compatibility\n",
        "heatmap_df = monthly_counts.pivot(index='Month', columns='Year', values='Total_Count')\n",
        "\n",
        "# Plot the heatmap\n",
        "plt.figure(figsize=(8, 3))\n",
        "sns.heatmap(heatmap_df, annot=True, cmap=\"YlGnBu\", fmt=\".0f\", linewidths=.5)\n",
        "# plt.title('Monthly Count of Motifs and Nearest Neighbors in Europe')\n",
        "plt.xlabel('Year')\n",
        "plt.xticks(rotation=45, fontsize=10)\n",
        "plt.ylabel('Month')\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "1qIAav50cGVO"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}