python/tutorials/Visualizing_Methane_Plume_Timeseries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing Methane Plume Timeseries\n",
    "\n",
    "**Summary**  \n",
    "\n",
    "In this notebook, we'll conduct a search and visualize the available methane plume complex observations detected by the Earth Mineral Dust Source Investigation (EMIT) instrument, then select a region of interest (ROI) that shows several plumes that appear to be from the same source. We will then visualize a time series of plume data for the ROI, calculate the integrated methane enhancement, or mass of methane observed in each plume. Lastly, to add more context to the plume timeseries, we will look at all EMIT observations over the target regions and determine if there were factors such as clouds limiting plume detection, or simply no methane emissions during additional overpasses.\n",
    "\n",
    "> Throughout this notebook, several complex functions and workflows are used to process and visualize data. These can be found in the `emit_tools.py` and `tutorial_utils.py` modules.\n",
    "\n",
    "**Background**\n",
    "\n",
    "Methane is major contributor to atmospheric radiative forcing because its more efficient at trapping radiation than other greenhouse gases, like carbon dioxide. Because the lifetime of methane in the atmosphere is only about 10 years, reducing methane emissions offers an effective way to curb anthropogenic contributions to atmospheric radiative forcing.\n",
    "\n",
    "The EMIT instrument is an imaging spectrometer that measures light in visible and infrared wavelengths. These measurements display unique spectral signatures that correspond to chemical composition. Although the primary mission focuses on mapping the mineral composition of Earth's surface, the EMIT instrument has also been used to successfully map methane point source emissions within its target mask using the unique spectral fingerprint of methane. More details about EMIT and its associated products can be found in the **README.md** and on the [EMIT website](https://earth.jpl.nasa.gov/emit/).\n",
    "\n",
    "The L2B Estimated Methane Plume Complexes ([EMITL2BCH4PLM](https://lpdaac.usgs.gov/products/emitl2bch4plmv001/)) product provides delineated methane plume complexes in parts per million meter (ppm m) along with associated uncertainties. This product is a plume specific subset of the EMIT L2B Methane Enhancement Data ([EMITL2BCH4ENH](https://lpdaac.usgs.gov/products/emitl2bch4enhv001/)) product. Each EMITL2BCH4PLM granule is sized to a specific plume complex but may cross multiple EMITL2BCH4ENH granules. A list of source EMITL2BCH4ENH granules is included in the GeoTIFF file metadata as well as the GeoJSON file. Each EMITL2BCH4PLM granule contains two files: one Cloud Optimized GeoTIFF (COG) file at a spatial resolution of 60 meters (m) and one GeoJSON file. The EMITL2BCH4PLM COG file contains a raster image of a methane plume complex extracted from EMITL2BCH4ENH v001 data. The EMITL2BCH4PLM GeoJSON file contains a vector outline of the plume complex, a list of source scenes, coordinates of the maximum enhancement values with associated uncertainties. \n",
    "\n",
    "The EMITL2BCH4ENH product only includes granules where methane plume complexes have been identified. To reduce the risk of false positives, all EMITL2BCH4ENH data undergo a manual review (or identification and confirmation) process before being designated as a plume complex. For more information on the manual review process, see Section 4.2.2 of the [EMIT GHG Algorithm Theoretical Basis Document (ATBD)](https://lpdaac.usgs.gov/documents/1696/EMIT_GHG_ATBD_V1.pdf). \n",
    "\n",
    "\n",
    "**References**\n",
    "\n",
    "Andrew K. Thorpe et al., Attribution of individual methane and carbon dioxide emission sources using EMIT observations from space. Sci. Adv.9, eadh2391 (2023). DOI:[10.1126/sciadv.adh2391](https://www.science.org/doi/10.1126/sciadv.adh2391)\n",
    "\n",
    "**Requirements** \n",
    " - Set up Python Environment - See **setup_instructions.md** in the `/setup/` folder\n",
    " - NASA Earthdata Login Account. [Sign Up](https://urs.earthdata.nasa.gov/users/new)\n",
    "\n",
    "**Data Used** \n",
    " - [EMIT L2B Estimated Methane Plume Complexes (EMITL2BCH4PLM)](https://lpdaac.usgs.gov/products/emitl2bch4plmv001/)\n",
    " - [EMIT L2A Estimated Surface Reflectance and Uncertainty and Masks (EMITL2ARFL)](https://lpdaac.usgs.gov/products/emitl2arflv001/)\n",
    "\n",
    "**Learning Objectives** \n",
    " - Search for EMIT L2B Estimated Methane Plume Complexes\n",
    " - Visualize search results on a map\n",
    " - Retrieve and visualize the EMIT L2B Estimated Methane Plume Complexes Metadata\n",
    " - Select a region of interest and build a timeseries of plume data\n",
    " - Further investigate plume detection by looking at browse images and quality information\n",
    "\n",
    "**Tutorial Outline**  \n",
    "\n",
    "1. [**Setup**](#setup)\n",
    "2. [**Search for EMIT L2B Estimated Methane Plume Complexes**](#search)\n",
    "3. [**Creating a Timeseries from Plume Data**](#plume-timeseries)\n",
    "4. [**Further investigation into plume detection**](#plume-detection)\n",
    "5. [**Calculating the Integrated Mass Enhancement for Plumes**](#ime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Set up<a id='setup'></a>\n",
    "\n",
    "Import the necessary Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required libraries\n",
    "import sys\n",
    "import os\n",
    "import glob\n",
    "import requests\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "from osgeo import gdal\n",
    "import geopandas as gpd\n",
    "\n",
    "from datetime import datetime\n",
    "import folium\n",
    "import earthaccess\n",
    "import folium.plugins\n",
    "import rioxarray as rxr\n",
    "\n",
    "import hvplot.xarray\n",
    "import hvplot.pandas\n",
    "\n",
    "from branca.element import Figure\n",
    "from IPython.display import display\n",
    "from shapely.geometry.polygon import orient\n",
    "from shapely.geometry import Point\n",
    "\n",
    "sys.path.append('../modules/')\n",
    "from emit_tools import emit_xarray, ortho_xr, ortho_browse\n",
    "from tutorial_utils import list_metadata_fields, results_to_geopandas, convert_bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Log into `earthaccess` and create a `.netrc` file if necessary. This file will store your NASA Earthdata Login credentials and use them to authenticate when necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earthaccess.login(persist=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the data we use or save will go into the `methane_tutorial` directory, so we can go ahead and define that filepath now, relative to this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "methane_dir = '../../data/methane_tutorial/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Search for EMIT L2B Estimated Methane Plume Complexes<a id='search'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `earthaccess` to find all EMIT L2B Estimated Methane Plume Complexes (EMITL2BCH4PLM) data available from 2023. Define the date range, and concept-ids (unique product identifier) for the EMIT products that we want to search for, but leave the spatial arguments like `polygon` and `bbox` empty so we can preview detected methane plumes globally. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data Collections for our search, using a dictionary\n",
    "concept_ids = {'plumes':'C2748088093-LPCLOUD', 'reflectance':'C2408750690-LPCLOUD'}\n",
    "# Define Date Range\n",
    "date_range = ('2023-01-01','2023-12-31')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = earthaccess.search_data(\n",
    "    concept_id=concept_ids['plumes'],\n",
    "    temporal=date_range,\n",
    "    count=2000\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the results to a `geopandas.GeoDataFrame` using a function from our `tutorial_utils` module. This gives a nice way to organize and visualize the search results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = results_to_geopandas(results)\n",
    "gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default this function includes some fields, but you can add fields with a `fields` argument. To see all of the metadata available use the `list_metadata_fields` function imported from the `tutorial_utils.py` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "list_metadata_fields(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add an index column to the dataframe to include it in the tooltips for our visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify index so we can reference it with gdf.explore()\n",
    "gdf['index']=gdf.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up Figure and Basemap tiles\n",
    "fig = Figure(width=\"1080px\",height=\"540\")\n",
    "map1 = folium.Map(tiles=None)\n",
    "folium.TileLayer(tiles='https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}',name='Google Satellite', attr='Google', overlay=True).add_to(map1)\n",
    "folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png',\n",
    "                name='ESRI World Imagery',\n",
    "                attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',\n",
    "                overlay='True').add_to(map1)\n",
    "fig.add_child(map1)\n",
    "# Add Search Results gdf\n",
    "gdf.explore(\"_single_date_time\",\n",
    "            categorical=True,\n",
    "            style_kwds={\"fillOpacity\":0.1,\"width\":2},\n",
    "            name=\"EMIT L2B CH4PLM\",\n",
    "            tooltip=[\n",
    "                \"index\",\n",
    "                \"native-id\",\n",
    "                \"_single_date_time\",\n",
    "            ],\n",
    "            m=map1,\n",
    "            legend=False\n",
    ")\n",
    "\n",
    "# Zoom to Data\n",
    "map1.fit_bounds(bounds=convert_bounds(gdf.unary_union.bounds))\n",
    "# Add Layer controls\n",
    "map1.add_child(folium.LayerControl(collapsed=False))\n",
    "display(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we've chosen a region that looks like it has a several plumes being emitted from the same source, a landfill in Jordan. To create a simple bounding box around our target region, we can use the plumes that extend furthest in the cardinal directions to generate a bounding box around the region that we can use in our upcoming analysis. \n",
    "\n",
    "Instead of this we will import an existing `geojson` as a `GeoDataFrame` with a bounding box around this region because the plume indices may change, which results in inconsistent outputs. A commented out cell below as included as an example of how the `GeoDataFrame` was created before being written to `geojson`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Select a list of plumes to create geometry we can use for a spatial subset\n",
    "# plumes = [146,198,243]\n",
    "# bbox = gdf.loc[plumes].geometry.unary_union.envelope\n",
    "# bbox = orient(bbox, sign=1)\n",
    "# plume_bbox = gpd.GeoDataFrame({\"name\":['plume_bbox'], \"geometry\":[bbox]},crs=gdf.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the predefined `geojson` of our plume bounding box as a `GeoDataFrame` then visualize on our `folium` figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the geojson with our plume bbox\n",
    "plume_bbox = gpd.read_file(f'{methane_dir}/plume_bbox.geojson')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up Figure and Basemap tiles\n",
    "fig = Figure(width=\"1080px\",height=\"540\")\n",
    "map1 = folium.Map(tiles=None)\n",
    "folium.TileLayer(tiles='https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}',name='Google Satellite', attr='Google', overlay=True).add_to(map1)\n",
    "folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png',\n",
    "                name='ESRI World Imagery',\n",
    "                attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',\n",
    "                overlay='True').add_to(map1)\n",
    "fig.add_child(map1)\n",
    "# Add Search Results gdf\n",
    "plume_bbox.explore(\"name\",\n",
    "                   name='Plume BBox',\n",
    "                   style_kwds={\"fillOpacity\":0,\"width\":2},\n",
    "                   m=map1,\n",
    "                   legend=False)\n",
    "\n",
    "gdf.explore(\"_single_date_time\",\n",
    "            categorical=True,\n",
    "            style_kwds={\"fillOpacity\":0.1,\"width\":2},\n",
    "            name=\"EMIT L2B CH4PLM\",\n",
    "            tooltip=[\n",
    "                \"index\",\n",
    "                \"native-id\",\n",
    "                \"_single_date_time\",\n",
    "            ],\n",
    "            m=map1,\n",
    "            legend=False\n",
    ")\n",
    "# Zoom to Data\n",
    "map1.fit_bounds(bounds=convert_bounds(plume_bbox.unary_union.bounds))\n",
    "# Add Layer controls\n",
    "map1.add_child(folium.LayerControl(collapsed=False))\n",
    "display(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subset our geodataframe of plumes to only those that intersect our bounding box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_gdf = gdf[gdf.geometry.intersects(plume_bbox.geometry[0])]\n",
    "plm_gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we look at an example of the `_related_urls` column in our geodataframe, we can see it contains the various links to assets associated with that plume."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_gdf['_related_urls'].iloc[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write a function to return the asset `URL` for a given asset and row in our dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_asset_url(row,asset, key='Type',value='GET DATA'):\n",
    "    \"\"\"\n",
    "    Retrieve a url from the list of dictionaries for a row in the _related_urls column.\n",
    "    Asset examples: CH4PLM, CH4PLMMETA, RFL, MASK, RFLUNCERT \n",
    "    \"\"\"\n",
    "    # Add _ to asset so string matching works\n",
    "    asset = f\"_{asset}_\"\n",
    "    # Retrieve URL matching parameters\n",
    "    for _dict in row['_related_urls']:\n",
    "        if _dict.get(key) == value and asset in _dict['URL'].split('/')[-1]:\n",
    "            return _dict['URL']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write another function to make a request to retrieve the json metadata for a given plume, using our `get_asset_url` function to select the correct URL from the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to fetch CH4 Plume Metadata\n",
    "# Speed could be improved here by using asyncio/aiohttp\n",
    "def fetch_ch4_metadata(row):\n",
    "    response = requests.get(get_asset_url(row, 'CH4PLMMETA'))\n",
    "    json = response.json()\n",
    "    return json['features'][0]['properties']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fetch_ch4_metadata(plm_gdf.iloc[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrieve additional plume metadata contained in the EMIT L2B Estimated Methane Plume Complexes (EMITL2BCH4PLM) data product, which contains the maximum enhancement value, the uncertainty of the plume complex, and the list of source scenes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the function to each row and convert the result to a DataFrame\n",
    "plm_meta = plm_gdf.apply(fetch_ch4_metadata, axis=1).apply(pd.Series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_meta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add the points with highest methane concentration to our visualization.\n",
    "\n",
    "Create an index column, as we did for the plumes, then convert the latitude and longitude of max concentration to a shapely `Point` object and add it to our `GeoDataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify index so we can reference it with gdf.explore()\n",
    "plm_meta['index'] = plm_meta.index\n",
    "# Add Geometry and convert to geodataframe\n",
    "plm_meta['geometry'] = plm_meta.apply(lambda row: Point(row['Longitude of max concentration'], row['Latitude of max concentration']), axis=1)\n",
    "plm_meta = gpd.GeoDataFrame(plm_meta, geometry='geometry', crs='EPSG:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_meta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now add this to our visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up Figure and Basemap tiles\n",
    "fig = Figure(width=\"1080px\",height=\"540\")\n",
    "map1 = folium.Map(tiles=None)\n",
    "folium.TileLayer(tiles='https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}',name='Google Satellite', attr='Google', overlay=True).add_to(map1)\n",
    "folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png',\n",
    "                name='ESRI World Imagery',\n",
    "                attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',\n",
    "                overlay='True').add_to(map1)\n",
    "fig.add_child(map1)\n",
    "# Add Search Results gdf\n",
    "plume_bbox.explore(\"name\",\n",
    "                   name='Plume BBox',\n",
    "                   style_kwds={\"fillOpacity\":0,\"width\":2},\n",
    "                   m=map1,\n",
    "                   legend=False)\n",
    "\n",
    "plm_gdf.explore(\"index\",\n",
    "            categorical=True,\n",
    "            style_kwds={\"fillOpacity\":0.1,\"width\":2},\n",
    "            name=\"EMIT L2B CH4PLM\",\n",
    "            tooltip=[\n",
    "                \"index\",\n",
    "                \"native-id\",\n",
    "                \"_single_date_time\",\n",
    "            ],\n",
    "            m=map1,\n",
    "            legend=False\n",
    ")\n",
    "\n",
    "plm_meta.explore(\"index\",\n",
    "            categorical=True,\n",
    "            style_kwds={\"fillOpacity\":0.1,\"width\":2},\n",
    "            name=\"Location of Max Concentration (ppm m)\",\n",
    "            tooltip=[\n",
    "                \"DAAC Scene Names\",\n",
    "                \"UTC Time Observed\",\n",
    "                \"Max Plume Concentration (ppm m)\",\n",
    "                \"Concentration Uncertainty (ppm m)\",\n",
    "                \"Orbit\"\n",
    "            ],\n",
    "            m=map1,\n",
    "            legend=False\n",
    ")\n",
    "# Zoom to Data\n",
    "map1.fit_bounds(bounds=convert_bounds(plume_bbox.unary_union.bounds))\n",
    "# Add Layer controls\n",
    "map1.add_child(folium.LayerControl(collapsed=False))\n",
    "display(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Creating a Timeseries from Plume Data<a id='plume-timeseries'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize a timeseries of these plumes that appear to be from the same source. To do this we'll generate a list of the COG urls for the plumes, then use rioxarray to stream the data and build a timeseries dataset based on the timestamps in the filenames."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `get_asset_url` function to retrieve the CH4PLM asset URLs by applying it to our dataframe and converting the output to a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Iterate over rows in the plm_gdf and get the CH4PLM urls and store them in a list\n",
    "plm_urls = plm_gdf.apply(lambda row: get_asset_url(row, asset='CH4PLM'), axis=1).tolist()\n",
    "plm_urls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a list of COG urls we can set our `gdal` configuration options to pass our NASA Earthdata login credentials when we access each COG."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GDAL configurations used to successfully access LP DAAC Cloud Assets via vsicurl \n",
    "gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt')\n",
    "gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt')\n",
    "gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR')\n",
    "gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF')\n",
    "gdal.SetConfigOption('GDAL_HTTP_UNSAFESSL', 'YES')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To build our timeseries we will start by opening all the necessary data. Loop over our list of urls, open each plume, merge plumes acquired at the same time, and store them in a dictionary where keys correspond to the acquisition time and values are the plume data in an `xarray.DataArray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_ts_dict = {}\n",
    "# Set max retries for vsicurl errors\n",
    "max_retries=5\n",
    "# Iterate over plm urls\n",
    "for url in plm_urls:\n",
    "    # retrieve acquisition time from url\n",
    "    acquisition_time = url.split('/')[-1].split('.')[-2].split('_')[-2]\n",
    "    # list plumes identified in same scene if there are any\n",
    "    same_scene = [url for url in plm_urls if acquisition_time in url.split('/')[-1].split('.')[-2].split('_')[-2]]\n",
    "    to_merge = []\n",
    "    # prevent duplicate processing of plumes from the same scene\n",
    "    if acquisition_time not in list(plm_ts_dict.keys()):\n",
    "        # Open and merge plumes identified from each scene\n",
    "        for _plm in same_scene:\n",
    "            print(f\"Opening {_plm.split('/')[-1]}\")\n",
    "            # Try loop for vsicurl/unrecongnized format error\n",
    "            for retry in range(max_retries):\n",
    "                try:\n",
    "                    # Open COG and squeeze band dimension\n",
    "                    plm = rxr.open_rasterio(_plm).squeeze('band', drop=True)\n",
    "                    break\n",
    "                except Exception as e:\n",
    "                    print(f'{e} Retrying...')\n",
    "                else:\n",
    "                    print(f\"Failed to process {url} after {max_retries} retries. Please check to see you're authenticated with earthaccess.\")\n",
    "            # Add to list of plumes to merge\n",
    "            to_merge.append(plm)\n",
    "            # Merge plumes and add to timeseries\n",
    "            plm_ts_dict[acquisition_time] = rxr.merge.merge_arrays(to_merge)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a plume object for each date in our timeseries, we need to put them all on a common grid so they are spatially aligned before we can stack them along the time dimension.\n",
    "\n",
    "To do this, we will find the minimum and maximum bounds of each dataarray in our dictionary, then create a common grid to reproject to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import List\n",
    "def create_common_grid(data_arrays: List[xr.DataArray]) -> xr.DataArray:\n",
    "    \"\"\"\n",
    "    Create a common grid for a list of xarray DataArrays, matching the resolution of the first data array.\n",
    "    \"\"\"\n",
    "    # Initial Bounds for common grid\n",
    "    minx = miny = float('inf')\n",
    "    maxx = maxy = float('-inf')\n",
    "\n",
    "    for array in data_arrays:\n",
    "        left, bottom, right, top = array.rio.bounds()\n",
    "        minx = min(minx, left)\n",
    "        miny = min(miny, bottom)\n",
    "        maxx = max(maxx, right)\n",
    "        maxy = max(maxy, top)\n",
    "\n",
    "    bounds = (minx, miny, maxx, maxy)\n",
    "\n",
    "    res = data_arrays[0].rio.resolution()\n",
    "    crs = data_arrays[0].rio.crs\n",
    "    nodata = data_arrays[0].rio.nodata\n",
    "\n",
    "    # # Calculate new raster shape using the new extent, maintaining the original resolution\n",
    "    height = int(np.ceil((bounds[3] - bounds[1]) / abs(res[1])))\n",
    "    width = int(np.ceil((bounds[2] - bounds[0]) / abs(res[0])))\n",
    "    data = np.full((height,width),nodata)\n",
    "    coords = {'y':(['y'],np.arange(bounds[1], bounds[3], abs(res[1]))),\n",
    "              'x':(['x'],np.arange(bounds[0], bounds[2], abs(res[0])))}\n",
    "    common_grid = xr.DataArray(data, coords=coords)\n",
    "    common_grid.rio.write_crs(crs, inplace=True)\n",
    "    return(common_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "common_grid = create_common_grid(list(plm_ts_dict.values()))\n",
    "common_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reproject each of the plume dataarrays to the common grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_ts_dict = {key: value.rio.reproject_match(common_grid) for key, value in plm_ts_dict.items()}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have all of our plumes on a standard grid, we can concatenate them along a time dimension to create a timeseries. Create an xarray variable called 'time' from our dictionary keys, then use `xarray.concat` to concatenate all of our plumes along the time dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_time = xr.Variable('time', [datetime.strptime(t,'%Y%m%dT%H%M%S') for t in list(plm_ts_dict.keys())])\n",
    "plm_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_ts_ds = xr.concat(list(plm_ts_dict.values()), dim=plm_time)\n",
    "plm_ts_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set our no_data values to `np.nan` to make sure they are transparent for improved visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_ts_ds.data[plm_ts_ds.data == -9999] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the plume time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_ts_plot = plm_ts_ds.hvplot.image(x='x',y='y',geo=True, tiles='ESRI', crs='EPGS:4326', cmap='inferno',clim=(0,np.nanmax(plm_ts_ds.data)),clabel=f'Methane Concentation ({plm_ts_ds.Units})', frame_width=600, frame_height=600, rasterize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_ts_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Calculating the IME for each plume<a id='ime'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integrated mass enhancement (IME), a sum of the mass of methane present in each plume is calculated by summing the mass of methane present in all plume pixels. The IME in kilograms (kg) can be estimated by using the equation below which includes the mixing ratio length per pixel (ppmm), the area of the pixel (m^2), where 22.4 is the volume of 1 mole of gas at standard temperature and pressure (STP) in liters, and 0.01604 is the molar mass of methane in kg/mol.\n",
    "\n",
    "The IME can be combined with a plume length and windspeed to estimate emissions in units of kg per hour, but in this tutorial, we will not calculate emissions, rather keep things simple and calculate the IME for each plume and observe how these values change over time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\ kg\\ \\ (per \\ \\ pixel) = \\frac{pixel \\ \\ value \\ \\ ppm \\cdot m}{1} \\frac{1}{1 \\cdot 10^6 \\ \\ ppm} \\frac {60 \\ \\ m \\cdot 60 \\ \\ m} {1} \\frac {1000 \\ \\ L} {m^3} \\frac {1 \\ \\ mol} {22.4 \\ \\ L} \\frac {0.01604 \\ \\ kg} {1 \\ \\ mol}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write this as a function for each pixel, then apply it to the entire timeseries to calculate the IME for each plume."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_ime(plume_da):\n",
    "    molar_volume = 22.4 # L/mol at STP\n",
    "    molar_mass_ch4 = 0.01604 #kg/mol\n",
    "\n",
    "    kg = plume_da * (1/1e6) * (60*60) * (1000) * (1/molar_volume) * molar_mass_ch4\n",
    "    ime = np.nansum(kg)\n",
    "    return ime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the function along the 'x' and 'y' dimensions\n",
    "ime_ts = xr.apply_ufunc(calc_ime, plm_ts_ds, input_core_dims=[['y', 'x']], vectorize=True)\n",
    "ime_ts.name = 'value'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ime_plot = ime_ts.hvplot.scatter(x='time',y='value', title='Observed Methane IME over 2023', color='black', xticks=list(ime_ts.time.data), rot=90, grid=True, xlabel='Date Observed', ylabel='kg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ime_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Further investigation into plume detection<a id='plume-detection'></a>\n",
    "\n",
    "The timeseries shown above doesn't necessarily give us a full a full picture methane emissions at the landfill. In addition to cases where no emissions were observed, gaps in data can result from absence of observations due the variable revisit period of the ISS, or clouds. To add more context to this plume timeseries, we will look at all EMIT acquisitions over the target regions and try to determine if there were factors affecting plume detection, or simply no methane emissions during those for any other overpass.\n",
    "\n",
    "For this search we'll want to restrict our search to a smaller ROI than our bounding box. We want to pick something smaller more centered around the source of the emission, so we avoid retrieving irrelevant data, for example an overpass barely crossing the corner of our large bounding box. To do this we will use the maximum concentration points from our plumes to create a smaller ROI that is likely to include a methane plume if one is present.\n",
    "\n",
    "Create a new polygon using the maximum concentration points from our plumes as vertices, then orienting the points in counter-clockwise order so they are compatible with an `earthaccess` search. This is just a simple example approach, there are several ways we could define a smaller ROI using the plume data or ancillary information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_conc_poly = plm_meta.geometry.unary_union.envelope\n",
    "max_conc_poly = orient(max_conc_poly, sign=1.0)\n",
    "max_conc_gdf = gpd.GeoDataFrame({\"name\":['max_points'], \"geometry\":[max_conc_poly]},crs=plm_meta.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roi = list(max_conc_gdf.geometry[0].exterior.coords)\n",
    "roi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now conduct a search for reflectance data over our ROI and convert the results to a `geopandas.GeoDataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rfl_results = earthaccess.search_data(\n",
    "    concept_id=concept_ids['reflectance'],\n",
    "    polygon=roi,\n",
    "    temporal= date_range, #('2023-03-01','2023-03-31')\n",
    "    count=2000\n",
    ")\n",
    "rfl_gdf = results_to_geopandas(rfl_results)\n",
    "rfl_gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the results we can see we have 20 scenes that intersect our ROI. We can visualize the footprints of these scenes to gain some insight into coverage over our ROI. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up Figure and Basemap tiles\n",
    "fig = Figure(width=\"1080px\",height=\"540\")\n",
    "map1 = folium.Map(tiles=None)\n",
    "folium.TileLayer(tiles='https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}',name='Google Satellite', attr='Google', overlay=True).add_to(map1)\n",
    "folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png',\n",
    "                name='ESRI World Imagery',\n",
    "                attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',\n",
    "                overlay='True').add_to(map1)\n",
    "fig.add_child(map1)\n",
    "\n",
    "# Add Search Reflectance Scenes with no CH4\n",
    "rfl_gdf.explore(color='red',\n",
    "               style_kwds={\"fillOpacity\":0,\"width\":2},\n",
    "               name=\"Scenes with no CH4 Plumes\",\n",
    "               tooltip=[\n",
    "                \"native-id\",\n",
    "                \"_beginning_date_time\",\n",
    "                ],\n",
    "                m=map1,\n",
    "                legend=False)\n",
    "\n",
    "# Add Plume BBox to Map\n",
    "max_conc_gdf.explore(m=map1,\n",
    "                   name='Plumes Bounding Box',\n",
    "                   legend=False)\n",
    "\n",
    "# Zoom to Data\n",
    "map1.fit_bounds(bounds=convert_bounds(rfl_gdf.unary_union.bounds))\n",
    "# Add Layer controls\n",
    "map1.add_child(folium.LayerControl(collapsed=False))\n",
    "display(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this we can see that there are several scenes that intersect with our ROI, but likely have no relevant information since they only cover a small portion or corner of the ROI."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use a similar process as with the methane product to construct a time series to better understand the data gathered on each overpass. To do this, we will use the browse imagery and the masks included in the EMITL2ARFL product. The by default the mask files and browse images are not orthorectified, so we must do that as part of our workflow.\n",
    "\n",
    "First, get the urls for the browse images and masks for each scene in our `rfl_gdf` search results using the `get_asset_urls` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_urls = rfl_gdf.apply(lambda row: get_asset_url(row, asset='RFL', value='GET RELATED VISUALIZATION'), axis=1).tolist()\n",
    "png_urls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_urls = rfl_gdf.apply(lambda row: get_asset_url(row, asset='MASK'), axis=1).tolist()\n",
    "mask_urls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write these as a text file so we don't need to search again, although we will use the `rfl_gdf` GeoDataFrame later in the tutorial. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Save URL List\n",
    "# with open(f'{methane_dir}rfl_mask_urls.txt', 'w') as f:\n",
    "#     for line in mask_urls:\n",
    "#         f.write(f\"{line}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the mask files are not chunked, its quicker to download them to do the processing.\n",
    "\n",
    "Login with `earthaccess` and download these files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earthaccess.login(persist=True)\n",
    "# Get requests https Session using Earthdata Login Info\n",
    "fs = earthaccess.get_requests_https_session()\n",
    "# Retrieve granule asset ID from URL (to maintain existing naming convention)\n",
    "for url in mask_urls:\n",
    "    granule_asset_id = url.split('/')[-1]\n",
    "    # Define Local Filepath\n",
    "    fp = f'{methane_dir}{granule_asset_id}'\n",
    "    # Download the Granule Asset if it doesn't exist\n",
    "    if not os.path.isfile(fp):\n",
    "        with fs.get(url,stream=True) as src:\n",
    "            with open(fp,'wb') as dst:\n",
    "                for chunk in src.iter_content(chunk_size=64*1024*1024):\n",
    "                    dst.write(chunk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each of these scenes we want to open EMIT L2A Mask data, then subset spatially and select only the variable we want, in this case we'll use the `Cloud flag` from the `masks` dataarray. There are other flags, such as a cirrus mask, dilated cloud flag, spacecraft flag, water flag, and and aerosol optical depth. We could potentially use some of these as well to inform our decisions about plume detection, but will stick to the `Cloud flag` in this example for simplicity.\n",
    "\n",
    "As we do this, we will also use the GLT included in the mask file to orthorectify our RGB browse image. We can do this because the browse png files are in the native resolution and can be broadcast onto an orthorectified grid using the GLT. We will use the `reproject_match` function to reproject the data on a `common_grid`, which will automatically clip the data to the `common_grid` extent.\n",
    "\n",
    "First, get the filepaths for our downloaded mask data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List the downloaded files\n",
    "fps = glob.glob(f'{methane_dir}*.nc')\n",
    "fns = [os.path.basename(fp) for fp in fps]\n",
    "fns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a function to loop through our files, orthorectifying the mask and browse image, clipping and reprojecting to our predefined `common_grid`, and finally saving outputs as a COG."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_scenes(fns, outdir, common_grid):\n",
    "    \"\"\"\n",
    "    This function will process a list of EMIT Mask scenes, selecting the cloud flag, orthorectifying the mask and browse image, then reprojecting both to a common grid and saving an output.\n",
    "    \"\"\"\n",
    "    for fn in fns:\n",
    "        # Get Granule Asset ID for First Adjacent Scene (may only be one)\n",
    "        granule_asset_id = fn.split('.')[-2]\n",
    "        # Set Output Path\n",
    "        outpath_mask = f\"{outdir}{granule_asset_id}_cloud_flag.tif\"\n",
    "        outpath_browse = f\"{outdir}{granule_asset_id}_ortho_browse.tif\"\n",
    "        # Check if the file exists\n",
    "        if not os.path.isfile(outpath_mask):\n",
    "            # Open Mask Dataset\n",
    "            emit_ds = emit_xarray(f'{methane_dir}{fn}', ortho=False)\n",
    "            # Retrieve GLT, spatial_ref, and geotransform to use on browse image\n",
    "            glt = np.nan_to_num(np.stack([emit_ds[\"glt_x\"].data, emit_ds[\"glt_y\"].data], axis=-1),nan=0).astype(int)\n",
    "            spatial_ref = emit_ds.spatial_ref\n",
    "            gt = emit_ds.geotransform\n",
    "            # Select browse image url corresponding to the scene\n",
    "            png_url = [url for url in png_urls if fn.split('.')[-2].split('_')[-3] in url][0]\n",
    "            # Orthorectify browse and mask\n",
    "            rgb = ortho_browse(png_url, glt, spatial_ref, gt, white_background=True)\n",
    "            emit_ds = ortho_xr(emit_ds)\n",
    "            # Select only mask array and desired quality flag and reproject to match our chosen extent\n",
    "            mask_da = emit_ds['mask'].sel(mask_bands='Cloud flag')\n",
    "            # Drop elevation\n",
    "            mask_da = mask_da.drop_vars('elev')\n",
    "            mask_da.name = 'Cloud flag'\n",
    "            mask_da.data = np.nan_to_num(mask_da.data, nan=-9999)\n",
    "            mask_da = mask_da.rio.reproject_match(common_grid, nodata=-9999)\n",
    "            #mask_da.rio.write_nodata(np.nan, inplace=True)\n",
    "            # Reproject rgb\n",
    "            rgb = rgb.rio.reproject_match(common_grid, nodata=255) # 255 for white background\n",
    "            # Write cog outputs        \n",
    "            mask_da.rio.to_raster(outpath_mask,driver=\"COG\")\n",
    "            rgb.rio.to_raster(outpath_browse,driver=\"COG\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "process_scenes(fns, methane_dir, common_grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a list of the processed files to use in creation of a timeseries. We'll use a similar process to what we did for the plumes, adding a `time` variable to our datasets and concatenating."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_files = sorted(glob.glob(f'{methane_dir}*cloud_flag.tif'))\n",
    "mask_files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_files = sorted(glob.glob(f'{methane_dir}*ortho_browse.tif'))\n",
    "rgb_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a time index from the filenames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_index_from_filenames(file_names,datetime_pos):\n",
    "    \"\"\"\n",
    "    Helper function to create a pandas DatetimeIndex\n",
    "    \"\"\"\n",
    "    return [datetime.strptime(f.split('_')[datetime_pos], '%Y%m%dT%H%M%S') for f in file_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_time = xr.Variable('time', time_index_from_filenames(mask_files, -5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open and concatenate our datasets along the time dimension, then assign fill_values to `np.nan` to make those sections of the data transparent in our visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quality_ts_da = xr.concat([rxr.open_rasterio(f).squeeze('band', drop=True).rio.reproject_match(common_grid) for f in mask_files], dim=mask_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a plot object for the quality timeseries, first setting the quality mask values representing no clouds (0) or no data (-9999) to `np.nan` so they will be transparent in our visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quality_ts_da.data[quality_ts_da.data < 1] = np.nan\n",
    "quality_ts_map = quality_ts_da.hvplot.image(x='x',y='y',cmap='greys',groupby='time',clim=(0,1),geo=True,frame_height=400)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RGB images are a good way to add something more visually understandable than just the mask layers. Follow the same process as above to build an RGB timeseries, then plot it with the bounding box and plume extents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_ts_ds = xr.concat([rxr.open_rasterio(f).rio.reproject_match(common_grid) for f in rgb_files], dim=mask_time)\n",
    "rgb_ts_ds.data[rgb_ts_ds.data == -1] = 255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_ts_map = rgb_ts_ds.hvplot.rgb(x='x',y='y', bands='band',groupby='time',geo=True, frame_height=400, crs='EPSG:4326')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, add a new column to our plume geodataframe named `time` so we are using the same naming convention and can layer our plume polygons from our filtered search on top of our quality data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plm_gdf['time'] = pd.to_datetime(plm_gdf.loc[:,'_single_date_time'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we've set our RGB to display as white where there is no data, and our cloud mask where no clouds are present as transparent, we can identify any areas with no data over our ROI by white/transparent, and cloudy areas as black. The larger black lines/rectangles are representative of onboard cloud masking, where no data was downlinked due to a high volume of clouds detected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_ts_map*quality_ts_map*plm_gdf.hvplot(groupby='time', geo=True, line_color='red', fill_color=None)*max_conc_gdf.hvplot(color='red',crs='EPSG:4326',fill_color=None, line_color='yellow')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this information, we can build a dataframe assigning a category to each of the dates where there was no plume detected and add this information to our IME timeseries figure.\n",
    "\n",
    "Create a dataframe of dates where no plume was detected by finding the dates where there are no plumes in our `mask_time` arrays by removing timestamps from our `plm_time` array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_plm_time = np.setdiff1d(mask_time.data,plm_time.data)\n",
    "no_plm_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a dataframe out of these dates where there were no plumes detected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a dataframe\n",
    "no_plm_df = pd.DataFrame({'time':no_plm_time})\n",
    "no_plm_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, categorize them based on the visualizations we made and remove any times where there wasn't a good observation of our area of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add a category column to describe the observation \n",
    "no_plm_df['category'] = 'no_data'\n",
    "no_plm_df.loc[[0,4],'category'] = 'cloud'\n",
    "no_plm_df.loc[7,'category'] = 'no_plume'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_plm_df = no_plm_df[~no_plm_df['category'].str.contains('no_data')]\n",
    "no_plm_df['time'] = pd.to_datetime(no_plm_df['time'])\n",
    "no_plm_df.reset_index(drop=True, inplace=True)\n",
    "no_plm_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now merge our IME timeseries dataframe with this one, along the time column to add these observations with different categories to our timeseries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ime_df = pd.merge(ime_ts.to_dataframe(), no_plm_df, on=['time'], how='outer')\n",
    "ime_df = ime_df.sort_values(by='time')\n",
    "ime_df['category']=ime_df['category'].fillna('plume')\n",
    "ime_df['value'] = ime_df['value'].fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ime_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot this IME timeseries alongside our plume timeseries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ime_timeline = ime_df.hvplot.scatter(x='time',y='value', by='category', size=100, xlabel='Date Observed', ylabel='kg', title='Observed Methane IME over 2023', rot=90,\n",
    "                                     grid=True, frame_height=400, frame_width=800, ylim=(0,11000),xticks=list(ime_df.time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ime_timeline.opts(legend_position='bottom') + plm_ts_plot.opts(title='Methane Plumes Observed', frame_height=400, frame_width=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contact Info:\n",
    "Email: LPDAAC@usgs.gov  \n",
    "Voice: +1-866-573-3222  \n",
    "Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹  \n",
    "Website: https://lpdaac.usgs.gov/  \n",
    "Date last modified: 03-22-2024  \n",
    "\n",
    "¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.  "
   ]
  }
 ],
 "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}