satellite_insitu_matchups.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "58a2b4bb",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Matchups of in situ data with satellite data\n",
    "\n",
    "**Tutorial Leads:** Anna Windle (NASA, SSAI), James Allen (NASA, MSU)\n",
    "\n",
    "## Summary\n",
    "\n",
    "In this example we will conduct matchups of in situ AERONET-OC Rrs data with PACE OCI Rrs data. The Aerosol Robotic Network (AERONET) was developed to sustain atmospheric studies at various scales with measurements from worldwide distributed autonomous sun-photometers. This has been extended to support marine applications, called AERONET – Ocean Color [(AERONET-OC)](https://aeronet.gsfc.nasa.gov/new_web/ocean_levels_versions.html), and provides the additional capability of measuring the radiance emerging from the sea (i.e., water-leaving radiance) with modified sun-photometers installed on offshore platforms like lighthouses, oceanographic and oil towers. AERONET-OC is instrumental in satellite ocean color validation activities. \n",
    "\n",
    "In this tutorial, we will be collecting Rrs data from the  [Casablanca Platform](https://aeronet.gsfc.nasa.gov/cgi-bin/data_display_seaprism_v3?site=Casablanca_Platform&nachal=2&level=3&place_code=10) AERONET-OC site located at 40.7N, 1.4W in the western Mediterranean Sea which is typically characterized as oligotrophic/mesotrophic (ocean color signals tend to strongly covary with chlorophyll a).\n",
    "\n",
    "We will be collecting PACE OCI Rrs data in a 5x5 pixel window around the AERONET-OC site to compare to the AERONET-OC Rrs data. \n",
    "\n",
    "## Learning Objectives\n",
    "\n",
    "At the end of this notebook you will know:\n",
    "\n",
    "* How to access Rrs data from a specific AERONET-OC site and time\n",
    "* How to access PACE OCI Rrs data from a specific location and time\n",
    "* How to match in situ and satellite data\n",
    "* How to apply statistics and plot matchup data\n",
    "\n",
    "## Contents\n",
    "\n",
    "1. [Setup](#1.-Setup)\n",
    "2. [Process AERONET-OC data](#2.-Process-AERONET-OC-data)\n",
    "3. [Process PACE OCI data](#3.-Process-PACE-OCI-data)\n",
    "4. [Apply matchup code](#4.-Apply-matchup-code)\n",
    "5. [Make plots](#5.-Make-plots)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb8c795c",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Setup\n",
    "\n",
    "We begin by loading a set of utility functions that work behind the scenes to do the majority of the work for us. \n",
    "Collapse this for now and just run the cell, but feel free to dig into it and see how things work!\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\"> \n",
    "    \n",
    "Note that the `get_f0` function requires the thuillier2003_f0.nc file. For the Hackweek, this is being pulled from the shared-public directory. For others trying to run the notebook, the .nc file can be found [here](https://oceancolor.gsfc.nasa.gov/docs/rsr/f0.txt). \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "03ffe055",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Helper functions for PACE Hackweek Validation Tutorial.\n",
    "\n",
    "Authors:\n",
    "    James Allen and Anna Windle\n",
    "\"\"\"\n",
    "import datetime\n",
    "from pathlib import Path\n",
    "import re\n",
    "\n",
    "import earthaccess\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.style as style\n",
    "from matplotlib.ticker import FuncFormatter\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pvlib.solarposition as sunpos\n",
    "from scipy import stats, odr\n",
    "import seaborn as sns\n",
    "import xarray as xr\n",
    "\n",
    "# AERONET-OC Download Constants\n",
    "# Valid AERONET-OC site list\n",
    "AERONET_SITES = [\n",
    "    'AAOT', 'Abu_Al_Bukhoosh', 'ARIAKE_TOWER', 'Bahia_Blanca', 'Banana_River',\n",
    "    'Blyth_NOAH', 'Casablanca_Platform', 'Chesapeake_Bay', 'COVE_SEAPRISM',\n",
    "    'Galata_Platform', 'Gloria', 'GOT_Seaprism', 'Grizzly_Bay',\n",
    "    'Gustav_Dalen_Tower', 'Helsinki_Lighthouse', 'Ieodo_Station',\n",
    "    'Irbe_Lighthouse', 'Kemigawa_Offshore', 'Lake_Erie', 'Lake_Okeechobee',\n",
    "    'Lake_Okeechobee_N', 'LISCO', 'Lucinda', 'MVCO', 'Palgrunden',\n",
    "    'PLOCAN_Tower', 'RdP-EsNM', 'Sacramento_River', 'San_Marco_Platform',\n",
    "    'Section-7_Platform', 'Socheongcho', 'South_Greenbay', 'Thornton_C-power',\n",
    "    'USC_SEAPRISM', 'Venise', 'WaveCIS_Site_CSI_6', 'Zeebrugge-MOW1'\n",
    "    ]\n",
    "\n",
    "# Get subset of AERONET columns to make it a bit more manageable (also rename)\n",
    "AOC_KEEP_COLS = [\"AERONET_Site\", \"aoc_datetime\", \"Site_Latitude(Degrees)\",\n",
    "                 \"Site_Longitude(Degrees)\", \"Solar_Zenith_Angle[400nm]\"]\n",
    "COLUMN_RENAME = {\n",
    "    \"Site_Latitude(Degrees)\": \"aoc_latitude\",\n",
    "    \"Site_Longitude(Degrees)\": \"aoc_longitude\",\n",
    "    \"AERONET_Site\": \"aoc_site\",\n",
    "    \"Solar_Zenith_Angle[400nm]\": \"aoc_solar_zenith\"\n",
    "    }\n",
    "\n",
    "# Bland-Altman/Scatterplot Constants\n",
    "# Plot colors, font sizes\n",
    "COLOR_PALETTE = sns.color_palette(\"colorblind\")\n",
    "COLOR_SCATTER = COLOR_PALETTE[0]\n",
    "COLOR_LINE = \"black\"  # Was \"black\"\n",
    "COLOR_LOA = COLOR_PALETTE[2]  # Was \"green\"\n",
    "COLOR_FITLINE = COLOR_PALETTE[1]  # Was \"magenta\"\n",
    "SIZE_TITLE = 24\n",
    "SIZE_AXLABEL = 20\n",
    "SIZE_TEXTLABEL = 14\n",
    "SHOW_LEGEND = False\n",
    "\n",
    "# Update some defaults\n",
    "plt.rcParams.update({\"figure.dpi\": 300})\n",
    "sns.set_style(\"ticks\", rc={\"figure.dpi\": 300})\n",
    "sns.set_context(\"notebook\", font_scale=1.45)\n",
    "\n",
    "# Satellite Matchup Constants\n",
    "# Short names for earthaccess lookup\n",
    "SAT_LOOKUP = {\n",
    "    \"PACE\": \"PACE_OCI_L2_AOP_NRT\",\n",
    "    \"AQUA\": \"MODISA_L2_OC\",\n",
    "    \"TERRA\": \"MODIST_L2_OC\",\n",
    "    \"NOAA-20\": \"VIIRSJ1_L2_OC\",\n",
    "    \"NOAA-21\": \"VIIRSJ2_L2_OC\",\n",
    "    \"SUOMI-NPP\": \"VIIRSN_L2_OC\"\n",
    "    }\n",
    "\n",
    "# List l2 flags, then build them into a dict\n",
    "l2_flags_list = [\n",
    "    \"ATMFAIL\", \"LAND\", \"PRODWARN\", \"HIGLINT\", \"HILT\", \"HISATZEN\", \"COASTZ\",\n",
    "    \"SPARE\", \"STRAYLIGHT\", \"CLDICE\", \"COCCOLITH\", \"TURBIDW\", \"HISOLZEN\",\n",
    "    \"SPARE\", \"LOWLW\", \"CHLFAIL\", \"NAVWARN\", \"ABSAER\", \"SPARE\", \"MAXAERITER\",\n",
    "    \"MODGLINT\", \"CHLWARN\", \"ATMWARN\", \"SPARE\", \"SEAICE\", \"NAVFAIL\", \"FILTER\",\n",
    "    \"SPARE\", \"BOWTIEDEL\", \"HIPOL\", \"PRODFAIL\", \"SPARE\"]\n",
    "L2_FLAGS = {flag: 1 << idx for idx, flag in enumerate(l2_flags_list)}\n",
    "\n",
    "# Bailey and Werdell 2006 exclusion criteria\n",
    "EXCLUSION_FLAGS = [\"LAND\", \"HIGLINT\", \"HILT\", \"STRAYLIGHT\", \"CLDICE\",\n",
    "                   \"ATMFAIL\", \"LOWLW\", \"FILTER\", \"NAVFAIL\", \"NAVWARN\"]\n",
    "\n",
    "# --------------------------------------------------------------------------- #\n",
    "#                              General Utilities                              #\n",
    "# --------------------------------------------------------------------------- #\n",
    "\n",
    "\n",
    "def get_f0(wavelengths=None, obs_time=None, window_size=10, f0_file=None):\n",
    "    \"\"\"\n",
    "    Load the Thuillier2003 netCDF file and return F0.\n",
    "\n",
    "    Defaults to returning the full table. Input obs_time to correct for the\n",
    "    Earth-Sun distance.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    wavelengths : array-like, optional\n",
    "        Wavelengths at which to compute the average irradiance.\n",
    "        If None, returns the full wavelength and irradiance table.\n",
    "    obs_time : datetime.datetime or pd.Series, optional\n",
    "        Observation time(s) used to correct for the Earth-Sun distance.\n",
    "        If None, return the mean F0 values.\n",
    "    window_size : int, optional\n",
    "        Bandpass filter size for mean filtering to selected wavelengths, in nm.\n",
    "    f0_file : str or pathlib.Path\n",
    "        Path to the f0 netCDF file of the lookup table.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    tuple of np.ndarray\n",
    "        A tuple containing:\n",
    "        - f0_spectra : np.ndarray\n",
    "            The solar irradiance values.\n",
    "        - f0_wave : np.ndarray\n",
    "            The corresponding wavelengths.\n",
    "    \"\"\"\n",
    "    if f0_file is None:\n",
    "        f0_file = Path(\"/home/jovyan/shared/pace-hackweek-2024/thuillier2003_f0.nc\")\n",
    "    f0_file = Path(f0_file)\n",
    "\n",
    "    if not f0_file.is_file():\n",
    "        raise FileNotFoundError(f\"File not found: {f0_file}\")\n",
    "\n",
    "    ds_f0 = xr.load_dataset(f0_file)\n",
    "    wl = ds_f0[\"wavelength\"].values\n",
    "    f0 = ds_f0[\"irradiance\"].values\n",
    "\n",
    "    if wavelengths is not None:\n",
    "        f0_wave = np.array(wavelengths)\n",
    "        f0_spectra = bandpass_avg(f0, wl, window_size, f0_wave)\n",
    "    else:\n",
    "        f0_wave = wl\n",
    "        f0_spectra = f0\n",
    "\n",
    "    if obs_time is not None:\n",
    "        # Calculate Earth-Sun distance\n",
    "        es_distance = sunpos.nrel_earthsun_distance(obs_time).to_numpy()\n",
    "\n",
    "        # Deal with multiple input times\n",
    "        if len(pd.Series(obs_time)) > 1:\n",
    "            f0_spectra = f0_spectra[None, :] / es_distance[:, None]**2\n",
    "        else:\n",
    "            f0_spectra /= es_distance**2\n",
    "\n",
    "    return f0_spectra, f0_wave\n",
    "\n",
    "\n",
    "def bandpass_avg(data, input_wavelengths, window_size=10,\n",
    "                 target_wavelengths=None):\n",
    "    \"\"\"\n",
    "    Apply a band-pass filter to the data.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : np.ndarray\n",
    "        1D or 2D array containing the spectral data (samples x wavelengths).\n",
    "        If 1D, it's assumed to be a single sample.\n",
    "    input_wavelengths : np.ndarray\n",
    "        1D array of wavelength values corresponding to the columns of data.\n",
    "    window_size : int, optional\n",
    "        Size of the window to use for averaging. Default is 10 nm.\n",
    "    target_wavelengths : np.ndarray, optional\n",
    "        1D array of target wavelengths for filtered values.\n",
    "        If None, the input wavelengths are used.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    np.ndarray\n",
    "        1D or 2D array containing the band-pass filtered data.\n",
    "    \"\"\"\n",
    "    data = np.atleast_2d(data)\n",
    "    half_window = window_size / 2\n",
    "    num_samples, num_input_wavelengths = data.shape\n",
    "    if target_wavelengths is None:\n",
    "        target_wavelengths = input_wavelengths\n",
    "\n",
    "    filtered_data = np.empty((num_samples, len(target_wavelengths)))\n",
    "\n",
    "    for idx, target_wl in enumerate(target_wavelengths):\n",
    "        start = target_wl - half_window\n",
    "        end = target_wl + half_window\n",
    "        cols_in_range = np.where((input_wavelengths >= start)\n",
    "                                 & (input_wavelengths <= end))[0]\n",
    "        filtered_data[:, idx] = np.nanmean(data[:, cols_in_range], axis=1)\n",
    "\n",
    "    return filtered_data if num_samples > 1 else filtered_data.flatten()\n",
    "\n",
    "\n",
    "def get_column_prods(df, type_prefix):\n",
    "    \"\"\"Process a dataframe to create a dictionary of data products.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    df : pandas DataFrame\n",
    "        Extracted dataframes from read_extract_file\n",
    "    type_prefix : str\n",
    "        Prefix to identify the product columns, e.g. \"aoc\"\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    data_dict\n",
    "        dictionary mapping data product with their wavelengths and columns.\n",
    "    \"\"\"\n",
    "    data_dict = {}\n",
    "    pattern = rf'{type_prefix}_(\\w+?)(\\d*\\.?\\d+)?$'\n",
    "\n",
    "    for col in df.columns:\n",
    "        match = re.match(pattern, col)\n",
    "        if match:\n",
    "            product = match.group(1)\n",
    "            wavelength = match.group(2) if match.group(2) else None\n",
    "            if product not in data_dict:\n",
    "                data_dict[product] = {'wavelengths': [], 'columns': []}\n",
    "            data_dict[product]['columns'].append(col)\n",
    "            if wavelength:\n",
    "                if '.' in wavelength:\n",
    "                    data_dict[product]['wavelengths'].append(float(wavelength))\n",
    "                else:\n",
    "                    data_dict[product]['wavelengths'].append(int(wavelength))\n",
    "    return data_dict\n",
    "\n",
    "# --------------------------------------------------------------------------- #\n",
    "#                            AERONET_OC Utilities                             #\n",
    "# --------------------------------------------------------------------------- #\n",
    "\n",
    "\n",
    "def construct_url(aoc_site, data_level, start_date, end_date):\n",
    "    \"\"\"\n",
    "    Craft the AERONET-OC data URL.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    aoc_site : str, optional\n",
    "        Specific AERONET-OC site (else AAOT by default)\n",
    "    start_date : datetime object, optional\n",
    "        Beginning of Aeronet data to run. Defaults to 1 Mar 2024.\n",
    "    end_date : datetime object, optional\n",
    "        End of Aeronet data to run. Defaults to today.\n",
    "    data_level : int, {10, 15, 20}\n",
    "        data quality; 20 (default, highest quality), 15, or 10.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    str\n",
    "        url to API pull\n",
    "    \"\"\"\n",
    "    # Validate inputs\n",
    "    if aoc_site not in AERONET_SITES:\n",
    "        raise ValueError(f\"{aoc_site} is not an AERONET site. Available \"\n",
    "                         f\"sites are: {', '.join(AERONET_SITES)}\")\n",
    "\n",
    "    url = (\"https://aeronet.gsfc.nasa.gov/cgi-bin/print_web_data_v3?\"\n",
    "           f\"AVG=10&LWN{data_level}=1&year={start_date.year}\"\n",
    "           f\"&month={start_date.month}&day={start_date.day}\"\n",
    "           f\"&if_no_html=1&year2={end_date.year}&month2={end_date.month}\"\n",
    "           f\"&day2={end_date.day}&site={aoc_site}\")\n",
    "\n",
    "    return url\n",
    "\n",
    "\n",
    "def get_data_dict(df, search_str=None):\n",
    "    \"\"\"Process a dataframe to create a dict and ndarray of products and waves.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    df : pandas DataFrame\n",
    "        Extracted dataframes from read_extract_file\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    wavelengths\n",
    "        numpy array of wavelengths of the data\n",
    "    column_map\n",
    "        dict of the dataframe columns associated with each wavelength\n",
    "    \"\"\"\n",
    "    if search_str is None:\n",
    "        search_str = \"Lwn_IOP\"\n",
    "    wavelengths = []\n",
    "    column_map = {}\n",
    "    pattern = re.compile(rf'{search_str}\\[(\\d+)nm\\]')\n",
    "\n",
    "    for col in df.columns:\n",
    "        match = pattern.search(col)\n",
    "        if match:\n",
    "            wavelength = int(match.group(1))\n",
    "            wavelengths.append(wavelength)\n",
    "            column_map[wavelength] = col\n",
    "    return np.array(wavelengths), column_map\n",
    "\n",
    "\n",
    "def process_aeronet(aoc_site=\"AAOT\", start_date=\"2024-03-01\", end_date=None,\n",
    "                    data_level=15):\n",
    "    \"\"\"\n",
    "    Download and process AERONET-OC data for matchups.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    aoc_site : str, optional\n",
    "        Specific AERONET-OC site (else AAOT by default)\n",
    "    start_date : datetime or str, optional\n",
    "        Beginning of Aeronet data to run. Defaults to \"2024-03-01\"\n",
    "    end_date : datetime or str, optional\n",
    "        End of Aeronet data to run. Defaults to today.\n",
    "    data_level : int, {10, 15, 20}\n",
    "        data quality; 20 (highest, but fewest), 15 (autochecked), or 10.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    pandas DataFrame object\n",
    "        Dataframe of downloaded AERONET-OC data\n",
    "    \"\"\"\n",
    "    # Set up processing\n",
    "    if end_date is None:\n",
    "        end_date = datetime.now()\n",
    "    start_date = pd.to_datetime(start_date, errors='raise')\n",
    "    end_date = pd.to_datetime(end_date, errors='raise')\n",
    "\n",
    "    # Make url\n",
    "    url_aoc = construct_url(aoc_site, data_level, start_date, end_date)\n",
    "\n",
    "    # Download data (skip the 5 header rows)\n",
    "    try:\n",
    "        df_aoc_full = pd.read_csv(url_aoc, delimiter=\",\", na_values=-999,\n",
    "                                  skiprows=5)\n",
    "    except Exception as e:\n",
    "        raise Exception(f\"Could not download data. Try another station, reduce\"\n",
    "                        f\" the data_level, or expand the times. (Error: {e})\")\n",
    "\n",
    "    # Drop empty columns\n",
    "    df_aoc_full.dropna(axis=1, how='all', inplace=True)\n",
    "\n",
    "    # Parse datetimes\n",
    "    df_aoc_full['aoc_datetime'] = pd.to_datetime(\n",
    "        df_aoc_full[\"Date(dd-mm-yyyy)\"] + ' ' + df_aoc_full[\"Time(hh:mm:ss)\"],\n",
    "        format=\"%d:%m:%Y %H:%M:%S\"\n",
    "        ).dt.tz_localize(\"UTC\")\n",
    "\n",
    "    # Get subset of Lwn_f/Q columns (ignore the count columns)\n",
    "    # Alternatively, could pull Lwn_IOP for L11 BRDF\n",
    "    subset_lwn = [col for col in df_aoc_full.columns\n",
    "                  if \"Lwn_f/Q\" in col and \"N[Lwn_f/Q\" not in col]\n",
    "    lwn_iop = df_aoc_full[subset_lwn].values\n",
    "\n",
    "    # Now get array of wavelengths from columns\n",
    "    wavelengths, _ = get_data_dict(df_aoc_full[subset_lwn], \"Lwn_f/Q\")\n",
    "\n",
    "    # Lwn need to be normalized by F0, the mean solar irradiance at top of atm\n",
    "    # Note: Lwn_IOP already accounts for the Earth-Sun Distance, BRDF, and\n",
    "    # atmosphere transmittance\n",
    "    f0_spectra, _ = get_f0(wavelengths)\n",
    "\n",
    "    # Normalize to get Rrs\n",
    "    rrs = lwn_iop / f0_spectra[None, :]\n",
    "\n",
    "    # Generate new column names and make the rrs dataframe\n",
    "    aoc_rrs_cols = [f\"aoc_rrs{wavelength}\" for wavelength in wavelengths]\n",
    "    df_rrs = pd.DataFrame(rrs, columns=aoc_rrs_cols)\n",
    "\n",
    "    # Now combine with the subset of the full dataframe\n",
    "    df_aoc = pd.concat([df_aoc_full[AOC_KEEP_COLS], df_rrs], axis=1)\n",
    "\n",
    "    # Do some final cleanup\n",
    "    df_aoc.rename(columns=COLUMN_RENAME, inplace=True)\n",
    "\n",
    "    return df_aoc\n",
    "\n",
    "\n",
    "# --------------------------------------------------------------------------- #\n",
    "#                             Satellite Utilities                             #\n",
    "# --------------------------------------------------------------------------- #\n",
    "\n",
    "\n",
    "def parse_quality_flags(flag_value):\n",
    "    \"\"\"\n",
    "    Parse bitwise flag into a list of flag names.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    flag_value : int\n",
    "        The integer representing the combined bitwise quality flags.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    list of str\n",
    "        List of flag names that are set in the flag_value.\n",
    "    \"\"\"\n",
    "    return [flag_name for flag_name, value in L2_FLAGS.items()\n",
    "            if (flag_value & value) != 0]\n",
    "\n",
    "\n",
    "def get_fivebyfive(file, latitude, longitude, rrs_wavelengths):\n",
    "    \"\"\"\n",
    "    Get stats on a 5x5 box around station coordinates of a satellite granule.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    file : earthaccess granule object\n",
    "        Satellite granule from earthaccess.\n",
    "    latitude : float\n",
    "        In decimal degrees for Aeronet-OC site for matchups\n",
    "    longitude : float\n",
    "        In decimal degrees (negative West) for Aeronet-OC site for matchups\n",
    "    rrs_wavelengths ; numpy array\n",
    "        Rrs wavelengths (from wavelength_3d for OCI)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    None.\n",
    "    \"\"\"\n",
    "    with xr.open_dataset(file, group=\"navigation_data\") as ds_nav:\n",
    "        sat_lat = ds_nav['latitude'].values\n",
    "        sat_lon = ds_nav['longitude'].values\n",
    "\n",
    "    # Calculate the Euclidean distance for 2D lat/lon arrays\n",
    "    distances = np.sqrt((sat_lat - latitude)**2 + (sat_lon - longitude)**2)\n",
    "\n",
    "    # Find the index of the minimum distance\n",
    "    # Dimensions are (lines, pixels)\n",
    "    min_dist_idx = np.unravel_index(np.argmin(distances), distances.shape)\n",
    "    center_line, center_pixel = min_dist_idx\n",
    "\n",
    "    # Get indices for a 5x5 box around the center pixel\n",
    "    line_start = max(center_line - 2, 0)\n",
    "    line_end = min(center_line + 2 + 1, sat_lat.shape[0])\n",
    "    pixel_start = max(center_pixel - 2, 0)\n",
    "    pixel_end = min(center_pixel + 2 + 1, sat_lat.shape[1])\n",
    "\n",
    "    # Extract the data\n",
    "    with xr.open_dataset(file, group=\"geophysical_data\") as ds_data:\n",
    "        rrs_data = ds_data['Rrs'].isel(\n",
    "            number_of_lines=slice(line_start, line_end),\n",
    "            pixels_per_line=slice(pixel_start, pixel_end)\n",
    "            ).values\n",
    "        flags_data = ds_data['l2_flags'].isel(\n",
    "            number_of_lines=slice(line_start, line_end),\n",
    "            pixels_per_line=slice(pixel_start, pixel_end)\n",
    "            ).values\n",
    "\n",
    "    # Calculate the bitwise OR of all flags in EXCLUSION_FLAGS to get a mask\n",
    "    exclude_mask = sum(L2_FLAGS[flag] for flag in EXCLUSION_FLAGS)\n",
    "\n",
    "    # Create a boolean mask\n",
    "    # True means the flag value does not contain any of the EXCLUSION_FLAGS\n",
    "    valid_mask = np.bitwise_and(flags_data, exclude_mask) == 0\n",
    "\n",
    "    # Get stats and averages\n",
    "    if valid_mask.any():\n",
    "        rrs_valid = rrs_data[valid_mask]\n",
    "        rrs_std_initial = np.std(rrs_valid, axis=0)\n",
    "        rrs_mean_initial = np.mean(rrs_valid, axis=0)\n",
    "\n",
    "        # Exclude spectra > 1.5 stdevs away\n",
    "        std_mask = np.all(\n",
    "            np.abs(rrs_valid - rrs_mean_initial) <= 1.5 * rrs_std_initial,\n",
    "            axis=1)\n",
    "        rrs_std = np.std(rrs_valid[std_mask], axis=0)\n",
    "        rrs_mean = np.mean(rrs_valid[std_mask], axis=0).flatten()\n",
    "\n",
    "        # Matchup criteria uses cv as median of 405-570nm\n",
    "        rrs_cv = rrs_std / rrs_mean\n",
    "        rrs_cv_median = np.median(rrs_cv[(rrs_wavelengths >= 405)\n",
    "                                         & (rrs_wavelengths <= 570)])\n",
    "    else:\n",
    "        rrs_cv_median = np.nan\n",
    "        rrs_mean = np.nan * np.empty_like(rrs_wavelengths)\n",
    "\n",
    "    # Put in dictionary of the row\n",
    "    row = {\n",
    "        \"oci_datetime\": pd.to_datetime(file.granule[\"umm\"][\"TemporalExtent\"]\n",
    "                                       [\"RangeDateTime\"][\"BeginningDateTime\"]),\n",
    "        \"oci_cv\": rrs_cv_median,\n",
    "        \"oci_latitude\": sat_lat[center_line, center_pixel],\n",
    "        \"oci_longitude\": sat_lon[center_line, center_pixel],\n",
    "        \"oci_pixel_valid\": np.sum(valid_mask)\n",
    "    }\n",
    "\n",
    "    # Add mean spectra to the row dictionary\n",
    "    for wavelength, mean_value in zip(rrs_wavelengths, rrs_mean):\n",
    "        key = f'oci_rrs{int(wavelength)}'\n",
    "        row[key] = mean_value\n",
    "\n",
    "    return row\n",
    "\n",
    "\n",
    "def process_satellite(start_date, end_date, latitude, longitude, sat=\"PACE\",\n",
    "                      selected_dates=None):\n",
    "    \"\"\"\n",
    "    Download and process satellite data for matchups.\n",
    "\n",
    "    Caution: If the date or coordinates aren't formatted correctly, it might\n",
    "    pull a huge granule list and take forever to run. If it takes more than 45\n",
    "    seconds to print the number of granules, just kill the process.\n",
    "\n",
    "    Uses the earthaccess package. Defaults to the PACE OCI L2 IOP datasets,\n",
    "    but other satellites can be used if they have a corresponding short_name\n",
    "    in the SAT_LOOKUP dictionary.\n",
    "\n",
    "    Workflow:\n",
    "        1. Get list of matchup granules\n",
    "        2. Loop through each file and:\n",
    "            2a. Find closest pixel to station, extract 5x5 pixel box\n",
    "            2b. Exclude pixels based on l2_flags\n",
    "            2c. Filtered mean to get single spectra\n",
    "            2d. Compute statistics and save data row\n",
    "        3. Organize output pandas dataframe\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    start_date : datetime or str\n",
    "        Beginning of Aeronet data to run.\n",
    "    end_date : datetime or str, optional\n",
    "        End of Aeronet data to run.\n",
    "    latitude : float\n",
    "        In decimal degrees for Aeronet-OC site for matchups\n",
    "    longitude : float\n",
    "        In decimal degrees (negative West) for Aeronet-OC site for matchups\n",
    "    sat : str\n",
    "        Name of satellite to search. Must be in SAT_LOOKUP dict constant.\n",
    "    selected_dates : list of str, optional\n",
    "        If given, only pull granules if the dates are in this list\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    pandas DataFrame object\n",
    "        Flattened table of all satellite granule matchups.\n",
    "\n",
    "    \"\"\"\n",
    "    # Look up short name from constants\n",
    "    if sat not in SAT_LOOKUP.keys():\n",
    "        raise ValueError(f\"{sat} is not in the lookup dictionary. Available \"\n",
    "                         f\"sats are: {', '.join(SAT_LOOKUP)}\")\n",
    "    short_name = SAT_LOOKUP[sat]\n",
    "\n",
    "    # Format search parameters\n",
    "    time_bounds = (f\"{start_date}T00:00:00\", f\"{end_date}T23:59:59\")\n",
    "\n",
    "    # Run Earthaccess data search\n",
    "    results = earthaccess.search_data(point=(longitude, latitude),\n",
    "                                      temporal=time_bounds,\n",
    "                                      short_name=short_name)\n",
    "    if selected_dates is not None:\n",
    "        filtered_results = [\n",
    "            result for result in results\n",
    "            if result[\"umm\"][\"TemporalExtent\"][\"RangeDateTime\"][\"BeginningDateTime\"][:10]\n",
    "            in selected_dates\n",
    "            ]\n",
    "        print(f\"Filtered to {len(filtered_results)} Granules.\")\n",
    "        files = earthaccess.open(filtered_results)\n",
    "    else:\n",
    "        files = earthaccess.open(results)\n",
    "\n",
    "    # Pull out Rrs wavelengths for easier processing\n",
    "    with xr.open_dataset(files[0], group=\"sensor_band_parameters\") as ds_bands:\n",
    "        rrs_wavelengths = ds_bands[\"wavelength_3d\"].values\n",
    "\n",
    "    # Loop through files and process\n",
    "    sat_rows = []\n",
    "    for idx, file in enumerate(files):\n",
    "        granule_date = pd.to_datetime(file.granule[\"umm\"][\"TemporalExtent\"]\n",
    "                                      [\"RangeDateTime\"][\"BeginningDateTime\"])\n",
    "        print(f\"Running Granule: {granule_date}\")\n",
    "        row = get_fivebyfive(file, latitude, longitude, rrs_wavelengths)\n",
    "        sat_rows.append(row)\n",
    "\n",
    "    return pd.DataFrame(sat_rows)\n",
    "\n",
    "# --------------------------------------------------------------------------- #\n",
    "#                              Matchup Utilities                              #\n",
    "# --------------------------------------------------------------------------- #\n",
    "\n",
    "\n",
    "def match_data(df_sat, df_aoc, cv_max=0.15, senz_max=60.0,\n",
    "               min_percent_valid=55.0, max_time_diff=180, std_max=1.5):\n",
    "    \"\"\"Create matchup dataframe based on selection criteria.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    df_sat : pandas dataframe\n",
    "        Satellite data from flat validation file.\n",
    "    df_aoc : pandas dataframe\n",
    "        Field data from flat validation file.\n",
    "    cv_max : float, default 0.15\n",
    "        Maximum coefficient of variation (stdev/mean) for sat data.\n",
    "    senz_max : float, default 60.0\n",
    "        Maximum sensor zenith for sat data.\n",
    "    min_percent_valid : float, default 55.0\n",
    "        Minimum percentage of valid satellite pixels.\n",
    "    max_time_diff : int, default 180\n",
    "        Maximum time difference (minutes) between sat and field matchup.\n",
    "    std_max : float, default 1.5\n",
    "        If multiple valid field matchups, select within std_max stdevs of mean.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    pandas dataframe of matchups for product\n",
    "    \"\"\"\n",
    "    # Setup\n",
    "    time_window = pd.Timedelta(minutes=max_time_diff)\n",
    "    df_match_list = []\n",
    "\n",
    "    # Filter Field data based on Solar Zenith\n",
    "    df_aoc_filtered = df_aoc[df_aoc['aoc_solar_zenith'] <= senz_max]\n",
    "\n",
    "    # Filter satellite data based on cv threshold\n",
    "    df_sat_filtered = df_sat[df_sat['oci_cv'] <= cv_max]\n",
    "\n",
    "    # Filter satellite data based on percent good pixels\n",
    "    df_sat_filtered = df_sat_filtered[\n",
    "        df_sat_filtered['oci_pixel_valid'] >= min_percent_valid * 25 / 100]\n",
    "\n",
    "    for _, sat_row in df_sat_filtered.iterrows():\n",
    "        # Filter field data based on time difference and coordinates\n",
    "        time_diff = abs(df_aoc_filtered['aoc_datetime']-sat_row['oci_datetime'])\n",
    "        within_time = time_diff <= time_window\n",
    "        within_lat = 0.2 >= abs(\n",
    "            df_aoc_filtered['aoc_latitude'] - sat_row['oci_latitude'])\n",
    "        within_lon = 0.2 >= abs(\n",
    "            df_aoc_filtered['aoc_longitude'] - sat_row['oci_longitude'])\n",
    "        field_matches = df_aoc_filtered[within_time & within_lat & within_lon]\n",
    "\n",
    "        if field_matches.shape[0] > 5:\n",
    "            # Filter by Standard Deviation for rrs columns\n",
    "            rrs_cols = [col for col in field_matches.columns\n",
    "                        if col.startswith('aoc_rrs')]\n",
    "            if rrs_cols:\n",
    "                mean_spectra = field_matches[rrs_cols].mean(axis=0)\n",
    "                std_spectra = field_matches[rrs_cols].std(axis=0)\n",
    "                within_std = (abs(field_matches[rrs_cols] - mean_spectra)\n",
    "                              <= std_max * std_spectra)\n",
    "                field_matches = field_matches[within_std.all(axis=1)]\n",
    "\n",
    "        if not field_matches.empty:\n",
    "            # Select the best match based on time delta\n",
    "            time_diff = abs(\n",
    "                field_matches['aoc_datetime']-sat_row['oci_datetime'])\n",
    "            best_match = field_matches.loc[time_diff.idxmin()]\n",
    "            df_match_list.append({**best_match.to_dict(), **sat_row.to_dict()})\n",
    "\n",
    "    df_match = pd.DataFrame(df_match_list)\n",
    "    return df_match\n",
    "\n",
    "\n",
    "# --------------------------------------------------------------------------- #\n",
    "#                              Plotting Utilities                             #\n",
    "# --------------------------------------------------------------------------- #\n",
    "\n",
    "\n",
    "def compute_bland_altman_metrics(xx, yy, xx_unc_modl, yy_unc_modl):\n",
    "    \"\"\"\n",
    "    Compute metrics for Bland-Altman plot.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    xx : array\n",
    "        Array of X data values.\n",
    "    yy : array\n",
    "        Array of Y data values.\n",
    "    xx_unc_modl : float\n",
    "        Uncertainty in X.\n",
    "    yy_unc_modl : float\n",
    "        Uncertainty in Y.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of Bland-Altman metrics.\n",
    "    \"\"\"\n",
    "    jj = (xx + yy) / 2\n",
    "    kk = (yy - xx) / np.sqrt((xx_unc_modl**2) + (yy_unc_modl**2))\n",
    "\n",
    "    meanbias = np.mean(kk)\n",
    "    stdbias = np.std(kk)\n",
    "    LOAlow = meanbias - stdbias\n",
    "    LOAhgh = meanbias + stdbias\n",
    "\n",
    "    ba_stat, ba_p = stats.spearmanr(jj, kk)\n",
    "    ba_independ = ba_p > 0.05\n",
    "\n",
    "    return {\n",
    "        \"count\": kk.shape[0],\n",
    "        \"jj\": jj,\n",
    "        \"kk\": kk,\n",
    "        \"meanbias\": meanbias,\n",
    "        \"LOAlow\": LOAlow,\n",
    "        \"LOAhgh\": LOAhgh,\n",
    "        \"ba_stat\": ba_stat,\n",
    "        \"ba_p\": ba_p,\n",
    "        \"ba_independ\": ba_independ\n",
    "    }\n",
    "\n",
    "\n",
    "def compute_regression_metrics(xx, yy, is_type2=False):\n",
    "    \"\"\"\n",
    "    Compute regression metrics using specified type.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    xx : array\n",
    "        Array of X data values.\n",
    "    yy : array\n",
    "        Array of Y data values.\n",
    "    is_type2 : bool, optional\n",
    "        Whether to use Type 2 regression (orthogonal distance regression).\n",
    "        Default is False, for Type 1 regression (ordinary least squares).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of regression metrics.\n",
    "    \"\"\"\n",
    "    if is_type2:\n",
    "        # Perform Type 2 regression (orthogonal distance regression)\n",
    "        def linear_model(B, x):\n",
    "            \"\"\"\n",
    "            Linear function y = m*x + b.\n",
    "\n",
    "            B is a vector of the parameters.\n",
    "            x is an array of the current x values.\n",
    "            x is in the same format as the x passed to Data or RealData.\n",
    "            Return an array in the same format as y passed to Data or RealData.\n",
    "            \"\"\"\n",
    "            return B[0] * x + B[1]\n",
    "\n",
    "        # Create a model instance\n",
    "        linear = odr.Model(linear_model)\n",
    "\n",
    "        # Create a RealData object using the data\n",
    "        data = odr.RealData(xx, yy)\n",
    "\n",
    "        # Set up ODR with the model and data\n",
    "        odr_instance = odr.ODR(data, linear, beta0=[1., 0.])\n",
    "\n",
    "        # Run the regression\n",
    "        odr_result = odr_instance.run()\n",
    "        slope = odr_result.beta[0]\n",
    "        intercept = odr_result.beta[1]\n",
    "    else:\n",
    "        # Perform Type 1 regression (ordinary least squares)\n",
    "        regress_result = stats.linregress(xx, yy)\n",
    "        slope = regress_result.slope\n",
    "        intercept = regress_result.intercept\n",
    "\n",
    "    spearman_r = stats.spearmanr(xx, yy)\n",
    "    pearson_r = stats.pearsonr(xx, yy)\n",
    "    rmse_all = np.sqrt(np.mean((yy - xx) ** 2))\n",
    "    mae_all = np.mean(np.abs(yy - xx))\n",
    "\n",
    "    return {\n",
    "        \"count\": len(xx),\n",
    "        \"slope\": slope,\n",
    "        \"intercept\": intercept,\n",
    "        \"r_spear\": spearman_r.correlation,\n",
    "        \"r_pear\": pearson_r[0],\n",
    "        \"rmse\": rmse_all,\n",
    "        \"mae\": mae_all\n",
    "    }\n",
    "\n",
    "\n",
    "def add_text_annotations(ax, text_lines, position='top right',\n",
    "                         fontsize=SIZE_TEXTLABEL):\n",
    "    \"\"\"\n",
    "    Add text annotations to the plot.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ax : Axes\n",
    "        The axis to add text to.\n",
    "    text_lines : list of str\n",
    "        List of strings to be displayed as text.\n",
    "    position : str, default 'top right'\n",
    "        Position of the text on the plot.\n",
    "    fontsize : int, default 12\n",
    "        Font size of the text.\n",
    "    \"\"\"\n",
    "    if position == 'top right':\n",
    "        x = 0.95\n",
    "        y = 0.95\n",
    "        ha = 'right'\n",
    "        va = 'top'\n",
    "    elif position == 'top left':\n",
    "        x = 0.05\n",
    "        y = 0.95\n",
    "        ha = 'left'\n",
    "        va = 'top'\n",
    "    elif position == 'bottom left':\n",
    "        x = 0.05\n",
    "        y = 0.05\n",
    "        ha = 'left'\n",
    "        va = 'bottom'\n",
    "    elif position == 'bottom right':\n",
    "        x = 0.95\n",
    "        y = 0.05\n",
    "        ha = 'right'\n",
    "        va = 'bottom'\n",
    "\n",
    "    text = '\\n'.join(text_lines)\n",
    "    ax.text(\n",
    "        x, y, text, transform=ax.transAxes, fontsize=fontsize,\n",
    "        verticalalignment=va, horizontalalignment=ha,\n",
    "        bbox=dict(facecolor='white', alpha=0.6, edgecolor='none')\n",
    "        )\n",
    "\n",
    "\n",
    "def setup_plot(label):\n",
    "    \"\"\"\n",
    "    Set up the plot with titles and labels.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    label : str\n",
    "        Title of the plot.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    tuple\n",
    "        Figure and axes of the plot.\n",
    "    \"\"\"\n",
    "    style.use('seaborn-v0_8-whitegrid')\n",
    "    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6), layout=\"constrained\")\n",
    "    fig.suptitle(label, fontsize=22)\n",
    "    return fig, ax1, ax2\n",
    "\n",
    "\n",
    "def format_ticks(ax):\n",
    "    \"\"\"Format the tick labels on the axes to be more readable.\"\"\"\n",
    "    ax.xaxis.set_major_formatter(FuncFormatter(lambda x, _: f'{x:.3g}'))\n",
    "    ax.yaxis.set_major_formatter(FuncFormatter(lambda y, _: f'{y:.3g}'))\n",
    "    ax.tick_params(axis='both', which='major', width=2, length=6)\n",
    "    ax.spines['top'].set_linewidth(2)\n",
    "    ax.spines['right'].set_linewidth(2)\n",
    "    ax.spines['left'].set_linewidth(2)\n",
    "    ax.spines['bottom'].set_linewidth(2)\n",
    "\n",
    "\n",
    "def plot_bland_altman(ax1, metrics, binscale, scat, xx_unc_modl,\n",
    "                      x_label=\"x\", y_label=\"y\"):\n",
    "    \"\"\"\n",
    "    Plot Bland-Altman plot.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ax1 : Axes\n",
    "        Axis for the Bland-Altman plot.\n",
    "    metrics : dict\n",
    "        Bland-Altman metrics.\n",
    "    binscale : float\n",
    "        Scaling factor for bin size.\n",
    "    scat : bool\n",
    "        If False, plot as 2D histogram.\n",
    "    xx_unc_modl : float\n",
    "        Uncertainty in X.\n",
    "    x_label : string, default \"x\"\n",
    "        String for labels for x data\n",
    "    y_label : string, default \"y\"\n",
    "        String for labels for y data\n",
    "    \"\"\"\n",
    "    jj = metrics[\"jj\"]\n",
    "    kk = metrics[\"kk\"]\n",
    "    npoints = metrics[\"count\"]\n",
    "    meanbias = metrics[\"meanbias\"]\n",
    "    LOAlow = metrics[\"LOAlow\"]\n",
    "    LOAhgh = metrics[\"LOAhgh\"]\n",
    "    ba_independ = metrics[\"ba_independ\"]\n",
    "    ba_stat = metrics[\"ba_stat\"]\n",
    "\n",
    "    nbin = int(0.5 * binscale * np.sqrt(len(jj)))\n",
    "    min_kk = meanbias - 5 * np.std(kk)\n",
    "    max_kk = meanbias + 5 * np.std(kk)\n",
    "\n",
    "    gamma = 0.5\n",
    "    if scat:\n",
    "        min_jj = np.min(jj)\n",
    "        max_jj = np.max(jj)\n",
    "        lineclr, loaclr, fitclr = (COLOR_LINE, COLOR_LOA, COLOR_FITLINE)\n",
    "        ax1.scatter(jj, kk, color=COLOR_SCATTER)\n",
    "        ax1.set_xlim([min_jj, max_jj])\n",
    "        ax1.set_ylim([min_kk, max_kk])\n",
    "    else:\n",
    "        jj_sorted = np.sort(jj)\n",
    "        min_jj = jj_sorted[int(0.01 * len(jj))]\n",
    "        max_jj = jj_sorted[int(0.99 * len(jj))]\n",
    "        lineclr, loaclr, fitclr = ('white', 'yellow', 'cyan')\n",
    "        h = ax1.hist2d(jj, kk, bins=(nbin, nbin),\n",
    "                       norm=mcolors.PowerNorm(gamma), cmap=plt.cm.inferno,\n",
    "                       range=[[min_jj, max_jj], [min_kk, max_kk]])\n",
    "        plt.colorbar(h[3], ax=ax1)\n",
    "\n",
    "    ax1.set_title('Bland-Altman plot', fontsize=SIZE_TITLE)\n",
    "    ylabel = ('Uncertainty normalized bias' if xx_unc_modl != np.sqrt(0.5)\n",
    "              else f'Bias, ${y_label}-{x_label}$')\n",
    "    ax1.set_ylabel(ylabel, fontsize=SIZE_AXLABEL)\n",
    "    ax1.set_xlabel(f'Paired mean, $({x_label}+{y_label})/2$',\n",
    "                   fontsize=SIZE_AXLABEL)\n",
    "    ax1.plot([min_jj, max_jj], [0, 0],\n",
    "             color=lineclr, linestyle='solid', linewidth=4.0)\n",
    "\n",
    "    if ba_independ:\n",
    "        ax1.plot([min_jj, max_jj], [meanbias, meanbias],\n",
    "                 color=fitclr, linestyle='dashed', linewidth=3.0,\n",
    "                 label='Mean Bias')\n",
    "        ax1.plot([min_jj, max_jj], [LOAlow, LOAlow],\n",
    "                 color=loaclr, linestyle='dashed', linewidth=2.0,\n",
    "                 label='Lower LOA')\n",
    "        ax1.plot([min_jj, max_jj], [LOAhgh, LOAhgh],\n",
    "                 color=loaclr, linestyle='dashed', linewidth=2.0,\n",
    "                 label='Upper LOA')\n",
    "        ax1.fill_between([min_jj, max_jj], LOAlow, LOAhgh,\n",
    "                         color=loaclr, alpha=0.1)\n",
    "    else:\n",
    "        ba_regress_result = stats.linregress(jj, kk)\n",
    "        ba_min_fit_yy = (ba_regress_result.slope * min_jj\n",
    "                         + ba_regress_result.intercept)\n",
    "        ba_max_fit_yy = (ba_regress_result.slope * max_jj\n",
    "                         + ba_regress_result.intercept)\n",
    "        ax1.plot([min_jj, max_jj], [ba_min_fit_yy, ba_max_fit_yy],\n",
    "                 color=fitclr, linestyle='dashed', linewidth=3.0,\n",
    "                 label='Linear Fit')\n",
    "    if SHOW_LEGEND:\n",
    "        ax1.legend()\n",
    "    ax1.grid(True)\n",
    "    format_ticks(ax1)\n",
    "\n",
    "    text_lines = [\n",
    "        f\"Number of Points: {npoints}\",\n",
    "        f\"Mean Bias: {meanbias:.2e}\",\n",
    "        f\"Limits of Agreement: [{LOAlow:.2e}, {LOAhgh:.2e}]\",\n",
    "        f\"Rank Correlation: {ba_stat:.3f}\",\n",
    "        \"Bias Independent\" if ba_independ else \"Bias Dependent\"\n",
    "    ]\n",
    "    add_text_annotations(ax1, text_lines, position='bottom right')\n",
    "\n",
    "\n",
    "def plot_scatter(ax2, xx, yy, regress_metrics, binscale, scat,\n",
    "                 x_label=\"x\", y_label=\"y\"):\n",
    "    \"\"\"\n",
    "    Plot scatter plot with regression line.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ax2 : Axes\n",
    "        Axis for the scatter plot.\n",
    "    xx : array\n",
    "        Array of X data values.\n",
    "    yy : array\n",
    "        Array of Y data values.\n",
    "    regress_metrics : dict\n",
    "        Regression metrics.\n",
    "    binscale : float\n",
    "        Scaling factor for bin size.\n",
    "    scat : bool\n",
    "        If False, plot as 2D histogram.\n",
    "    x_label : string, default \"x\"\n",
    "        String for labels for x data\n",
    "    y_label : string, default \"y\"\n",
    "        String for labels for y data\n",
    "    \"\"\"\n",
    "    nbin = int(0.5 * binscale * np.sqrt(len(xx)))\n",
    "    min_val = min(np.min(xx), np.min(yy))\n",
    "    max_val = max(np.max(xx), np.max(yy))\n",
    "    gamma = 0.5\n",
    "\n",
    "    if scat:\n",
    "        ax2.scatter(xx, yy, color=COLOR_SCATTER)\n",
    "        ax2.set_xlim([min_val, max_val])\n",
    "        ax2.set_ylim([min_val, max_val])\n",
    "    else:\n",
    "        g = ax2.hist2d(xx, yy, bins=(nbin, nbin),\n",
    "                       norm=mcolors.PowerNorm(gamma), cmap=plt.cm.inferno,\n",
    "                       range=[[min_val, max_val], [min_val, max_val]])\n",
    "        plt.colorbar(g[3], ax=ax2)\n",
    "\n",
    "    ax2.set_title('Scatterplot', fontsize=SIZE_TITLE)\n",
    "    ax2.set_xlabel(f'${x_label}$', fontsize=SIZE_AXLABEL)\n",
    "    ax2.set_ylabel(f'${y_label}$', fontsize=SIZE_AXLABEL)\n",
    "    ax2.plot([min_val, max_val], [min_val, max_val],\n",
    "             color=COLOR_LINE, linestyle='solid', linewidth=4.0)\n",
    "\n",
    "    slope = regress_metrics[\"slope\"]\n",
    "    intercept = regress_metrics[\"intercept\"]\n",
    "    min_fit_yy = slope * min_val + intercept\n",
    "    max_fit_yy = slope * max_val + intercept\n",
    "    ax2.plot([min_val, max_val], [min_fit_yy, max_fit_yy],\n",
    "             color=COLOR_FITLINE, linestyle='dashed', linewidth=3.0,\n",
    "             label='Regression Line')\n",
    "    if SHOW_LEGEND:\n",
    "        ax2.legend()\n",
    "    ax2.grid(True)\n",
    "    format_ticks(ax2)\n",
    "\n",
    "    text_lines = [\n",
    "        f\"Slope: {slope:.3f}\",\n",
    "        f\"Intercept: {intercept:.2e}\",\n",
    "        f\"Linear Correlation: {regress_metrics['r_pear']:.3f}\",\n",
    "        f\"Rank Correlation: {regress_metrics['r_spear']:.3f}\",\n",
    "        f\"RMSE: {regress_metrics['rmse']:.2e}\",\n",
    "        f\"MAE: {regress_metrics['mae']:.2e}\"\n",
    "    ]\n",
    "    add_text_annotations(ax2, text_lines, position='bottom right')\n",
    "\n",
    "\n",
    "def plot_BAvsScat(x_input, y_input, label='',\n",
    "                  saveplot=None, scat=True, binscale=1.0,\n",
    "                  xx_unc_modl=np.sqrt(0.5), yy_unc_modl=np.sqrt(0.5),\n",
    "                  x_label=\"x\", y_label=\"y\", is_type2=True):\n",
    "    \"\"\"\n",
    "    Routine to plot paired data as Bland-Altman and scatter plot.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x_input : array-like\n",
    "        Array of X data values.\n",
    "    y_input : array-like\n",
    "        Corresponding array of Y data values.\n",
    "    label : string, default ''\n",
    "        Text label for plotting.\n",
    "    saveplot : string, default None\n",
    "        Set to save plot in ../output/ with the string as the filename.\n",
    "    scat : boolean, default True\n",
    "        Make a 2D histogram if False, regular scatter plot if True.\n",
    "    binscale : float, default 1.0\n",
    "        Scaling factor for how many bins to include in a 2D histogram.\n",
    "    xx_unc_modl : float, default np.sqrt(0.5)\n",
    "        Uncertainty in X.\n",
    "    yy_unc_modl : float, default np.sqrt(0.5)\n",
    "        Uncertainty in Y.\n",
    "    x_label : string, default \"x\"\n",
    "        String for labels for x data\n",
    "    y_label : string, default \"y\"\n",
    "        String for labels for y data\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of computed statistics.\n",
    "    \"\"\"\n",
    "    xx = np.asarray(x_input)\n",
    "    yy = np.asarray(y_input)\n",
    "    valid_indices = (np.isfinite(x_input) & np.isfinite(y_input)\n",
    "                     & (x_input != -999) & (y_input != -999))\n",
    "    xx = x_input[valid_indices]\n",
    "    yy = y_input[valid_indices]\n",
    "\n",
    "    ba_metrics = compute_bland_altman_metrics(xx, yy, xx_unc_modl, yy_unc_modl)\n",
    "    regress_metrics = compute_regression_metrics(xx, yy, is_type2=is_type2)\n",
    "\n",
    "    fig, ax1, ax2 = setup_plot(label)\n",
    "    plot_bland_altman(ax1, ba_metrics, binscale, scat, xx_unc_modl,\n",
    "                      x_label, y_label)\n",
    "    plot_scatter(ax2, xx, yy, regress_metrics, binscale, scat,\n",
    "                 x_label, y_label)\n",
    "\n",
    "    if saveplot is not None:\n",
    "        figpath = Path(\"../output\") / saveplot\n",
    "        fig.savefig(figpath)\n",
    "        print('Saved figure to:', figpath)\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "    return {\n",
    "        \"Number_of_Points\": ba_metrics[\"count\"],\n",
    "        \"Scale_Independence\": ba_metrics[\"ba_independ\"],\n",
    "        \"Mean_Bias\": ba_metrics[\"meanbias\"],\n",
    "        \"Limits_of_Agreement_low\": (ba_metrics[\"LOAlow\"]\n",
    "                                    if ba_metrics[\"ba_independ\"]\n",
    "                                    else float(\"nan\")),\n",
    "        \"Limits_of_Agreement_high\": (ba_metrics[\"LOAhgh\"]\n",
    "                                     if ba_metrics[\"ba_independ\"]\n",
    "                                     else float(\"nan\")),\n",
    "        \"Linear_Slope\": regress_metrics[\"slope\"],\n",
    "        \"Linear_Intercept\": regress_metrics[\"intercept\"],\n",
    "        \"Linear_Correlation\": regress_metrics[\"r_pear\"],\n",
    "        \"Rank_Correlation\": regress_metrics[\"r_spear\"],\n",
    "        \"RMSE\": regress_metrics[\"rmse\"],\n",
    "        \"MAE\": regress_metrics[\"mae\"]\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "485f411e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Process AERONET-OC data\n",
    "\n",
    "We will use the function `process_aeronet` to download and process AERONET-OC data from the 'Casablanca_Platform' site. We will filter Level 1.5 data from the dates June 1, 2024 to July 31, 2024. This function will output a pandas dataframe of every AERONET-OC record between the dates.\n",
    "\n",
    "There are three \"levels\" of AERONET-OC data in terms of data quality: 1, 1.5, and 2. If a complete measurement sequence with the instruments is able to be performed, it is collected and stored as Level 1. These data are then passed through an automated quality control system and stored as Level 1.5 if they pass all tests. Finally, Level 2 data are data from Level 1.5 that are subsequently screened by an experienced scientist and validated. We'll be using Level 1.5 data to pull as much good quality data as possible without the time lag for manual validation. More information on AERONET-OC levels can be found in [Zibordi et al., 2009.](https://doi.org/10.1175/2009JTECHO654.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "779dc1d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "aoc_cb = process_aeronet(aoc_site=\"Casablanca_Platform\", \n",
    "                start_date=\"2024-06-01\", end_date=\"2024-07-31\",\n",
    "                data_level=15)\n",
    "aoc_cb.head()              "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43487309",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3. Process PACE OCI data \n",
    "\n",
    "We will use the function `process_satellite` to search for `PACE_OCI_L2_AOP_NRT` data using `earthacces` within the specified time range and at the (lat,lon) coordinate of the Casablanca_Platform AERONET-OC site. This function finds the closest pixel and extracts all data within a 5x5 pixel window, excludes pixels based on L2 flags, calculates the mean to retrive a single Rrs spectra, and computes matchup statistics. The function outputs a pandas dataframe of every `PACE_OCI_L2_AOP_NRT` Rrs spectra for the specified time range. We'll also include an optional list of unique date strings from the AERONET-OC dataframe to \"skip\" the granules that don't have any field data associated with them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee43678f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pull out coordinates \n",
    "aoc_lat = aoc_cb[\"aoc_latitude\"][0]\n",
    "aoc_lon = aoc_cb[\"aoc_longitude\"][0]\n",
    "\n",
    "# Pull out unique days\n",
    "unique_days = aoc_cb[\"aoc_datetime\"].dt.date.unique()\n",
    "unique_days_str = [day.strftime('%Y-%m-%d') for day in unique_days]\n",
    "\n",
    "sat_cb = process_satellite(start_date=\"2024-06-01\", end_date=\"2024-07-31\",\n",
    "                  latitude=aoc_lat, longitude=aoc_lon, sat=\"PACE\",\n",
    "                  selected_dates=unique_days_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "049a4721",
   "metadata": {},
   "outputs": [],
   "source": [
    "sat_cb.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d79e0152",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3. Apply matchup code\n",
    "\n",
    "We will use the function `match_data` to create a matchup dataframe based on selection criteria. This function defaults to using the [Bailey and Werdell 2006](https://oceancolor.gsfc.nasa.gov/staff/jeremy/bailey_and_werdell_2006_rse.pdf) matchup criteria, which reduces the measurements made at a given station to one representative sample for validating against the satellite spectra. Data are filtered based on the solar zenith angle, their noise level, and the time difference (here 180 minutes from the satellite overpass). Potential satellite matchups are also reduced based on the signal to noise level of the 5x5 pixel aggregation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3616803",
   "metadata": {},
   "outputs": [],
   "source": [
    "?match_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83cec894",
   "metadata": {},
   "outputs": [],
   "source": [
    "matchups = match_data(sat_cb, aoc_cb, cv_max=0.60, senz_max=60.0, \n",
    "                      min_percent_valid=55.0, max_time_diff=180, std_max=1.5)\n",
    "matchups "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62eb0c22",
   "metadata": {},
   "source": [
    "Pull out wavelengths and Rrs data from matchups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a10a4fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_aoc = get_column_prods(matchups, \"aoc\")\n",
    "waves_aoc = np.array(dict_aoc[\"rrs\"][\"wavelengths\"])\n",
    "rrs_aoc = matchups[dict_aoc[\"rrs\"][\"columns\"]].to_numpy()\n",
    "\n",
    "dict_sat = get_column_prods(matchups, \"oci\")\n",
    "waves_sat = np.array(dict_sat[\"rrs\"][\"wavelengths\"])\n",
    "rrs_sat = matchups[dict_sat[\"rrs\"][\"columns\"]].to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43630057",
   "metadata": {},
   "source": [
    "## 4. Make plots\n",
    "\n",
    "We will use the function `plot_BAvsScat` to plot the paired matchup data as Bland_Altman and scatter plots. The Bland-Altman plots provide insights into the bias and precision of the satellite measurements compared to field measurements. A mean difference close to zero indicates good agreement, while the spread of differences (limits of agreement) puts the bias within the context of the variability of the field data. Additionally, a check is done to assess the scale dependency of the bias, such as errors increasing when the magnitude of the observations increases. If a scale dependency exists, the limits of agreement are replaced with a regression line showing its direction and magnitude. Scatter plots complement Bland-Altman plots by showing the strength of the linear relationship between the two datasets, with high correlation coefficients and low RMSE values indicating strong agreement and high accuracy of the satellite-derived measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "107f9e32",
   "metadata": {},
   "outputs": [],
   "source": [
    "?plot_BAvsScat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d24f0ed6",
   "metadata": {},
   "outputs": [],
   "source": [
    "MATCH_WAVES = np.array([400, 412, 443, 490, 510, 560, 620, 667])\n",
    "\n",
    "# Loop through matchup wavelengths\n",
    "stats_list = []\n",
    "for idx, match_wave in enumerate(MATCH_WAVES):\n",
    "    # Find matching OCI columns\n",
    "    idx_sat = np.where(np.abs(waves_sat - match_wave) <= 5)[0]\n",
    "    match_sat = np.nanmean(rrs_sat[:, idx_sat], axis=1)\n",
    "\n",
    "    # Find matching AOC columns\n",
    "    idx_aoc = np.where(np.abs(waves_aoc - match_wave) <= 5)[0]\n",
    "    match_aoc = np.nanmean(rrs_aoc[:, idx_aoc], axis=1)\n",
    "\n",
    "    valid_indices = np.isfinite(match_sat) & np.isfinite(match_aoc)\n",
    "    if np.sum(valid_indices) > 5:\n",
    "        fig_label = f\"Rrs({match_wave}), sr\" + r\"$\\mathregular{^{-1}}$\"\n",
    "        dict_stats = plot_BAvsScat(match_aoc[valid_indices],\n",
    "                                   match_sat[valid_indices],\n",
    "                                   label=fig_label,\n",
    "                                   saveplot=None,\n",
    "                                   x_label=\"AERONET\", y_label=\"PACE\",\n",
    "                                   is_type2=True)\n",
    "        dict_stats[\"wavelength\"] = match_wave\n",
    "        stats_list.append(dict_stats)\n",
    "\n",
    "# Organize stats DataFrame\n",
    "df_stats = pd.DataFrame(stats_list)\n",
    "df_stats.set_index('wavelength', inplace=True)\n",
    "df_stats = df_stats.fillna(-999)\n",
    "df_stats"
   ]
  }
 ],
 "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}