{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Rivulet: U. S. NOAA CoastWatch and OceanWatch\n",
    "_by Michelle H Wilkerson and Adelmo Eloy_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## Purpose of this Notebook\n",
    "\n",
    "This notebook is focused on **Sea Level** and **Ocean Health** as phenomena. Many students, particularly those living in coastal areas, understand that there are rising sea levels. However, they may not be aware of the different causes for sea level rise which include but are not limited to ocean temperature and thermal expansion, ice melt, and local weather and wind conditions. Similarly, they may not know that sea level rise has different impacts in different coastal regions.\n",
    "\n",
    "This data tool allows users to connect to the NOAA CoastWatch and OceanWatch ERDDAP servers, which provide a wealth of oceanographic data including sea level, temperature, and chlorophyll levels. You are able to search for data streams in an area of interest, identify a date range of interest, and compare historical patterns to see how our oceans are changing over time. The datasets this notebook is constructed to fetch can serve as a launching point for examining the underlying mechanistic and compositional complexities of our marine environments.\n",
    "\n",
    "<details>\n",
    "    <summary>Click here for more information</summary>\n",
    "\n",
    "API stands for **Application Programming Interface**. Think of it as a language to communicate with data centers so you can search and get the data you need. For students, using APIs is like having a direct line to the most advanced scientific sensors on the planet. It allows us to work with the same **live, real-world data** that professional scientists use to track what's going on in the world around us.\n",
    "\n",
    "You are welcome to modify and adapt this script. You may find the CoastWatch documentation [here](https://coastwatch.noaa.gov/) and the `erddapy` documentation [here](https://ioos.github.io/erddapy/) helpful. We use a version of NOAA CoastWatch data that is quality controlled for scientific research in this notebook. But, it is possible to access real time data if you are interested in fetching a dataset that is related to a current event of interest by selecting a different DATASET_ID. A list of all datasets is here: https://coastwatch.noaa.gov/erddap/info/index.html?page=1&itemsPerPage=2000 (Note not all datasets are publicly accessible)\n",
    "\n",
    "This notebook was developed as part of NSF Grant 2445609 to support accessing and processing public datasets for middle and high school classroom activities. It's written to be relatively accessible to beginners, but if you have not interacted with computational notebooks or Python before, you may find navigating this tool difficult. (Check out the [Show Your Work](https://github.com/CalCoRE/show-your-work) project for a gentle introduction to computational notebooks for educators!)\n",
    "\n",
    "Our project is focused on supporting data analysis and mechanistic reasoning in science education. In other words, we want students to learn how data provides information about _how scientific mechanisms work_ and how understanding scientific mechanisms can help them to _explain and interpret patterns in data_. This builds on a long history of research on complex systems and agent-based modeling, and more closely connects that work to current expansions of data analysis across subjects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Part I: Setup\n",
    "\n",
    "First, you need to connect to the API and specify what region and time period you are interested in. This section will help you with that."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "### Connecting with ERDDAP\n",
    "\n",
    "ERDDAP is a \"data broker\" tool that helps you connect to a variety of datasets related to ocean health. Here, we'll start by making sure you've installed erddapy, the library tool that lets us send instructions to ERDDAP through this notebook, and ensure we can connect.\n",
    "\n",
    "#### **Note for Google CoLab users:**\n",
    "If you are running this notebook in Google CoLab, uncomment the first line of the code below. This will make it easier for you to download your data files. There will also be lines to uncomment in each \"Data Fetch\" section, once you are ready to get your data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install erddapy\n",
    "\n",
    "# from google.colab import files ## UNCOMMENT THIS LINE IF USING COLAB\n",
    "from erddapy import ERDDAP #this will get us the data\n",
    "import pandas as pd #this will help us work with the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "Now, we'll connect to the CoastWatch dataset to make sure it is ready. Most of the sections below use CoastWatch data, though some also add other data to help you explore relationships."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "### Customizing Your Location"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "Let's identify a \"home base\" location which will help you make sure the data you get is related to your students. This dataset divides locations into a \"grid\" of points, and reports data on each point in the grid. So, it's good to start small since any box is likely to have multiple \"gridpoints\" of data that we will be loading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_lat = 25.0\n",
    "max_lat = 27.0\n",
    "min_long = -85.0\n",
    "max_long = -83.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "Below, we'll create a map with the box you defined, to make sure you're capturing what you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install folium #install the mapping library\n",
    "import folium\n",
    "\n",
    "bbox = [[min_lat, min_long], [max_lat, max_long]]\n",
    "\n",
    "# Calculate the center of the box to position the map\n",
    "map_center = [(bbox[0][0] + bbox[1][0]) / 2, (bbox[0][1] + bbox[1][1]) / 2]\n",
    "\n",
    "# Create a Folium map object\n",
    "m = folium.Map(location=map_center, zoom_start=8)\n",
    "\n",
    "# Add a rectangle for the bounding box to the map\n",
    "folium.Rectangle(\n",
    "    bounds=bbox,\n",
    "    color=\"#ff0000\",        # Red border\n",
    "    fill=True,\n",
    "    fill_color=\"#ff7800\",   # Orange fill\n",
    "    fill_opacity=0.2\n",
    ").add_to(m)\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "### Specify Target Date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime, timedelta\n",
    "\n",
    "TARGET_DATE_STR = '2023-12-15' \n",
    "target_date = datetime.strptime(TARGET_DATE_STR, '%Y-%m-%d')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "# Part II: Fetch Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "## Comparing Sea Levels in Recent Decades"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Here, we are looking at sea levels as determined through \"Absolute Dynamic Topography\" (ADT). This is a measurement that adjusts for irregularities in the Earth's gravity that might pull on some parts of the ocean more than other, \"squeezing\" water and making height comparisons inaccurate. ADT is considered to be an especially good way to look at how thermal expansion and added water and contributing to overall sea level height in different parts of the world."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "### Connect to the Right Server"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "erddap_obj = ERDDAP(\n",
    "    server='https://erddap.aoml.noaa.gov/hdb/erddap', \n",
    "    protocol='griddap',\n",
    "    response='nc'\n",
    ")\n",
    "\n",
    "erddap_obj.get_info_url()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "If you ran the code above and saw a blue URL appear, that means have successfully identified and connected to the database. Let's get started!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "### Get Data for the Target Date\n",
    "\n",
    "First, let's use the target date you identified above as a starting point. The data we'll be looking at is reported daily, so we want enough days that some short-term event (like a cold snap) doesn't throw off our readings too much. Let's start with a window of 20 days, 10 before our target date and 10 after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "window_start = target_date - timedelta(days=10)\n",
    "window_end = target_date + timedelta(days=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "erddap_obj.dataset_id = f'SEA_SURFACE_HEIGHT_NRT_{str(target_date.year)}'\n",
    "\n",
    "erddap_obj.griddap_initialize()\n",
    "\n",
    "erddap_obj.constraints['time>='] = window_start\n",
    "erddap_obj.constraints['time<='] = window_end\n",
    "erddap_obj.constraints['latitude>='] = min_lat\n",
    "erddap_obj.constraints['latitude<='] = max_lat\n",
    "erddap_obj.constraints['longitude>='] = min_long\n",
    "erddap_obj.constraints['longitude<='] = max_long\n",
    "\n",
    "erddap_obj.variables = ['adt']\n",
    "\n",
    "latest_data = erddap_obj.to_pandas()\n",
    "\n",
    "latest_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "### Get Data From 10 Years Ago"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "Now let's fetch some data from exactly 10 years prior to the target date for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "ten_years_ago = target_date.year - 10\n",
    "\n",
    "ten_years_earlier = target_date.replace(year=ten_years_ago)\n",
    "\n",
    "window_start = ten_years_earlier - timedelta(days=10)\n",
    "window_end = ten_years_earlier + timedelta(days=10)\n",
    "\n",
    "DATASET_ID = f'SEA_SURFACE_HEIGHT_{str(ten_years_ago)}_v3'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "erddap_obj.dataset_id = DATASET_ID\n",
    "\n",
    "erddap_obj.griddap_initialize()\n",
    "\n",
    "erddap_obj.constraints['time>='] = window_start\n",
    "erddap_obj.constraints['time<='] = window_end\n",
    "erddap_obj.constraints['latitude>='] = min_lat\n",
    "erddap_obj.constraints['latitude<='] = max_lat\n",
    "erddap_obj.constraints['longitude>='] = min_long\n",
    "erddap_obj.constraints['longitude<='] = max_long\n",
    "erddap_obj.variables = ['adt']\n",
    "\n",
    "ten_year_data = erddap_obj.to_pandas()\n",
    "\n",
    "ten_year_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "The dataset starts in 1993, so for most recent-ish years, you can get data from 10 and 20 years ago. Below, we repeat the process one more time to fetch data from 20 years prior to the target date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "twenty_years_ago = target_date.year - 20\n",
    "erddap_obj.dataset_id = f'SEA_SURFACE_HEIGHT_{str(twenty_years_ago)}_v3'\n",
    "\n",
    "twenty_years_earlier = target_date.replace(year=twenty_years_ago)\n",
    "\n",
    "window_start = twenty_years_earlier - timedelta(days=5)\n",
    "window_end = twenty_years_earlier + timedelta(days=5)\n",
    "\n",
    "erddap_obj.griddap_initialize()\n",
    "\n",
    "erddap_obj.constraints['time>='] = window_start\n",
    "erddap_obj.constraints['time<='] = window_end\n",
    "erddap_obj.constraints['latitude>='] = min_lat\n",
    "erddap_obj.constraints['latitude<='] = max_lat\n",
    "erddap_obj.constraints['longitude>='] = min_long\n",
    "erddap_obj.constraints['longitude<='] = max_long\n",
    "erddap_obj.variables = ['adt']\n",
    "\n",
    "twenty_year_data = erddap_obj.to_pandas()\n",
    "twenty_year_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sla_groups = [\n",
    "    ('Latest (2025)', latest_data['adt (m)'].dropna()),\n",
    "    ('10 years ago (2015)', ten_year_data['adt (m)'].dropna()),\n",
    "    ('20 years ago (2005)', twenty_year_data['adt (m)'].dropna()),\n",
    "]\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharey=True)\n",
    "\n",
    "for ax, (label, values) in zip(axes, sla_groups):\n",
    "    ax.boxplot(values, patch_artist=True, boxprops=dict(facecolor='#8fbce6', color='#2a5d8f'))\n",
    "    ax.set_title(label)\n",
    "    ax.set_ylabel('SLA (m)' if ax is axes[0] else '')\n",
    "    ax.grid(True, axis='y', linestyle='--', alpha=0.5)\n",
    "\n",
    "fig.suptitle('Sea Level Anomaly Distributions: Latest, 10 Years Ago, 20 Years Ago')\n",
    "plt.tight_layout(rect=[0, 0, 1, 0.95])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29",
   "metadata": {},
   "source": [
    "Let's merge all the datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_adt = pd.concat(\n",
    "    [\n",
    "        latest_data,\n",
    "        ten_year_data,\n",
    "        twenty_year_data,\n",
    "    ],\n",
    "    ignore_index=True\n",
    ")\n",
    "\n",
    "merged_adt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_adt.to_csv('sealevelhistoric.csv')\n",
    "\n",
    "# files.download('sealevelhistoric.csv') ## UNCOMMENT THIS LINE IF USING COLAB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32",
   "metadata": {},
   "source": [
    "## Sea Surface Temperature (SST) and Thermal Expansion\n",
    "\n",
    "Why does the sea level rise? One of the primary mechanisms is **thermal expansion**. Just like the liquid in a thermometer rises as it gets warmer, the water in the ocean takes up more space as its temperature increases. \n",
    "\n",
    "By looking at Sea Surface Temperature (SST) alongside sea level data, students can investigate this relationship. First, let's get sea level data for the region and time period you are interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [],
   "source": [
    "window_start = target_date - timedelta(days=10)\n",
    "window_end = target_date + timedelta(days=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34",
   "metadata": {},
   "source": [
    "First, let's make sure we're connected to the right server to access the database we need for sea level data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [],
   "source": [
    "erddap_obj = ERDDAP(\n",
    "    server='https://erddap.aoml.noaa.gov/hdb/erddap', \n",
    "    protocol='griddap',\n",
    "    response='nc'\n",
    ")\n",
    "\n",
    "erddap_obj.get_info_url()\n",
    "\n",
    "erddap_obj.dataset_id = f'SEA_SURFACE_HEIGHT_NRT_{str(target_date.year)}'\n",
    "\n",
    "erddap_obj.griddap_initialize()\n",
    "\n",
    "erddap_obj.constraints['time>='] = window_start\n",
    "erddap_obj.constraints['time<='] = window_end\n",
    "erddap_obj.constraints['latitude>='] = min_lat\n",
    "erddap_obj.constraints['latitude<='] = max_lat\n",
    "erddap_obj.constraints['longitude>='] = min_long\n",
    "erddap_obj.constraints['longitude<='] = max_long\n",
    "\n",
    "erddap_obj.variables = ['adt']\n",
    "\n",
    "latest_data = erddap_obj.to_pandas()\n",
    "\n",
    "latest_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36",
   "metadata": {},
   "source": [
    "Now, we'll use a similar set of commands to connect to a different server where the sea surface temperature information is held, and get that data for the same region and time period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "erddap_sst = ERDDAP(\n",
    "    server='https://coastwatch.noaa.gov/erddap',\n",
    "    protocol='griddap',\n",
    "    response='nc'\n",
    ")\n",
    "\n",
    "erddap_sst.dataset_id = 'noaacwBLENDEDCsstDaily'\n",
    "erddap_sst.griddap_initialize()\n",
    "\n",
    "# Use the same spatial and temporal constraints as before\n",
    "erddap_sst.constraints['time>='] = window_start\n",
    "erddap_sst.constraints['time<='] = window_end\n",
    "erddap_sst.constraints['latitude>='] = min_lat\n",
    "erddap_sst.constraints['latitude<='] = max_lat\n",
    "erddap_sst.constraints['longitude>='] = min_long\n",
    "erddap_sst.constraints['longitude<='] = max_long\n",
    "\n",
    "sst_data = erddap_sst.to_pandas()\n",
    "sst_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38",
   "metadata": {},
   "source": [
    "Now, let's merge the datasets so we can explore the relationship between the sea surface temperature (SST) and the absolute dynamic topograph height (ADT)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [],
   "source": [
    "latest_for_merge = latest_data.copy()\n",
    "latest_for_merge[\"time_day\"] = pd.to_datetime(latest_for_merge[\"time (UTC)\"]).dt.normalize()\n",
    "\n",
    "sst_for_merge = sst_data[[\n",
    "    \"latitude (degrees_north)\",\n",
    "    \"longitude (degrees_east)\",\n",
    "    \"analysed_sst (degree_C)\",\n",
    "    \"analysis_error (kelvin)\",\n",
    "]].copy()\n",
    "sst_for_merge[\"time_day\"] = pd.to_datetime(sst_data[\"time (UTC)\"]).dt.normalize()\n",
    "\n",
    "merged_latest_sst = pd.merge(\n",
    "    latest_for_merge,\n",
    "    sst_for_merge,\n",
    "    on=[\"time_day\", \"latitude (degrees_north)\", \"longitude (degrees_east)\"],\n",
    "    how=\"inner\"\n",
    ").drop(columns=[\"time_day\"])\n",
    "\n",
    "# Keep all original columns from latest_data and add only the requested SST columns\n",
    "merged_latest_sst = merged_latest_sst[\n",
    "    list(latest_data.columns) + [\"analysed_sst (degree_C)\", \"analysis_error (kelvin)\"]\n",
    "]\n",
    "\n",
    "merged_latest_sst"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40",
   "metadata": {},
   "source": [
    "Let's take a quick look at the scatterplot of sea level height versus sea surface temperature. Remember, this dataset is for a grid of points within the location you identified above, with daily measurements of both temperature and sea level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 5))\n",
    "plt.scatter(\n",
    "    merged_latest_sst['analysed_sst (degree_C)'],\n",
    "    merged_latest_sst['adt (m)'],\n",
    "    alpha=0.6\n",
    ")\n",
    "plt.xlabel('Sea Surface Temperature (°C)')\n",
    "plt.ylabel('ADT (m)')\n",
    "plt.title('Sea Level Height vs. Sea Surface Temperature')\n",
    "plt.grid(True, linestyle='--', alpha=0.5)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42",
   "metadata": {},
   "source": [
    "If this looks good, you can download the file for further analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_latest_sst.to_csv('latest_sst_adt.csv')\n",
    "\n",
    "# files.download('latest_sst_adt.csv') ## UNCOMMENT THIS LINE IF USING COLAB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44",
   "metadata": {},
   "source": [
    "## Single dataset for multiple years\n",
    "\n",
    "Run the code below if you want to provide a get a single dataset with ADT and SST around a given data for multiple years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45",
   "metadata": {},
   "outputs": [],
   "source": [
    "import io\n",
    "import requests\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from erddapy import ERDDAP\n",
    "\n",
    "# ============================================================\n",
    "# USER PARAMETERS\n",
    "# ============================================================\n",
    "\n",
    "# Choose a month-day target. The year here is only used to get month/day.\n",
    "# Example: compare the same seasonal window around July 15 across many years.\n",
    "target_month_day = \"07-01\"\n",
    "\n",
    "# Window around the target month-day\n",
    "days_before = 30\n",
    "days_after = 30\n",
    "\n",
    "# Option A: explicitly choose years\n",
    "selected_years = [2003, 2013, 2023]\n",
    "\n",
    "# Option B: generate years automatically\n",
    "# Uncomment this if you prefer n years ending in a given year.\n",
    "# end_year = 2023\n",
    "# n_years = 21\n",
    "# selected_years = list(range(end_year - n_years + 1, end_year + 1))\n",
    "\n",
    "# Spatial bounds\n",
    "lat_min, lat_max = sorted([min_lat, max_lat])\n",
    "lon_min, lon_max = sorted([min_long, max_long])\n",
    "\n",
    "# SST grid rounding used for merging with ADT\n",
    "coord_rounding = 1\n",
    "\n",
    "\n",
    "# ============================================================\n",
    "# HELPER FUNCTIONS\n",
    "# ============================================================\n",
    "\n",
    "def make_target_window(year, target_month_day, days_before=15, days_after=15):\n",
    "    \"\"\"\n",
    "    Create a time window around the same month-day for a given year.\n",
    "    \"\"\"\n",
    "    target_date = pd.to_datetime(f\"{year}-{target_month_day}\")\n",
    "    window_start = target_date - pd.Timedelta(days=days_before)\n",
    "    window_end = target_date + pd.Timedelta(days=days_after)\n",
    "    return target_date, window_start, window_end\n",
    "\n",
    "\n",
    "def get_adt_dataset_id(year):\n",
    "    \"\"\"\n",
    "    AOML sea-surface-height datasets use different naming patterns:\n",
    "    historical years use SEA_SURFACE_HEIGHT_<year>_v3;\n",
    "    recent years use SEA_SURFACE_HEIGHT_NRT_<year>.\n",
    "    \"\"\"\n",
    "    if year < 2022:\n",
    "        return f\"SEA_SURFACE_HEIGHT_{year}_v3\"\n",
    "    else:\n",
    "        return f\"SEA_SURFACE_HEIGHT_NRT_{year}\"\n",
    "\n",
    "\n",
    "def download_adt_for_year(year, window_start, window_end,\n",
    "                          lat_min, lat_max, lon_min, lon_max):\n",
    "    \"\"\"\n",
    "    Download ADT from AOML ERDDAP for one year/window.\n",
    "    \"\"\"\n",
    "    erddap_obj = ERDDAP(\n",
    "        server=\"https://erddap.aoml.noaa.gov/hdb/erddap\",\n",
    "        protocol=\"griddap\",\n",
    "        response=\"nc\"\n",
    "    )\n",
    "\n",
    "    dataset_id = get_adt_dataset_id(year)\n",
    "    print(f\"Downloading ADT: {dataset_id}\")\n",
    "\n",
    "    erddap_obj.dataset_id = dataset_id\n",
    "\n",
    "    erddap_obj.constraints[\"time>=\"] = window_start.strftime(\"%Y-%m-%dT%H:%M:%SZ\")\n",
    "    erddap_obj.constraints[\"time<=\"] = window_end.strftime(\"%Y-%m-%dT%H:%M:%SZ\")\n",
    "    erddap_obj.constraints[\"latitude>=\"] = lat_min\n",
    "    erddap_obj.constraints[\"latitude<=\"] = lat_max\n",
    "    erddap_obj.constraints[\"longitude>=\"] = lon_min\n",
    "    erddap_obj.constraints[\"longitude<=\"] = lon_max\n",
    "\n",
    "    erddap_obj.variables = [\"adt\"]\n",
    "\n",
    "    adt_data = erddap_obj.to_pandas()\n",
    "\n",
    "    adt_data[\"time (UTC)\"] = pd.to_datetime(\n",
    "        adt_data[\"time (UTC)\"],\n",
    "        errors=\"coerce\"\n",
    "    )\n",
    "\n",
    "    return adt_data\n",
    "\n",
    "\n",
    "def erddap_time(x):\n",
    "    \"\"\"\n",
    "    Format datetime for direct ERDDAP griddap URL.\n",
    "    \"\"\"\n",
    "    return pd.to_datetime(x).strftime(\"%Y-%m-%dT%H:%M:%SZ\")\n",
    "\n",
    "\n",
    "def download_sst_for_window(window_start, window_end,\n",
    "                            lat_min, lat_max, lon_min, lon_max,\n",
    "                            time_stride=1, lat_stride=1, lon_stride=1):\n",
    "    \"\"\"\n",
    "    Download SST from NOAA CoastWatch ERDDAP using direct CSV URL.\n",
    "    This avoids erddapy metadata initialization issues for this dataset.\n",
    "    \"\"\"\n",
    "    server_url = \"https://coastwatch.noaa.gov/erddap/griddap\"\n",
    "    dataset_id = \"noaacwBLENDEDCsstDaily\"\n",
    "\n",
    "    brackets = (\n",
    "        f\"[({erddap_time(window_start)}):{time_stride}:({erddap_time(window_end)})]\"\n",
    "        f\"[({lat_min}):{lat_stride}:({lat_max})]\"\n",
    "        f\"[({lon_min}):{lon_stride}:({lon_max})]\"\n",
    "    )\n",
    "\n",
    "    variables = [\"analysed_sst\", \"analysis_error\"]\n",
    "    query = \",\".join([f\"{var}{brackets}\" for var in variables])\n",
    "\n",
    "    url = f\"{server_url}/{dataset_id}.csv?{query}\"\n",
    "\n",
    "    print(\"Downloading SST\")\n",
    "\n",
    "    response = requests.get(\n",
    "        url,\n",
    "        headers={\"User-Agent\": \"Mozilla/5.0\"},\n",
    "        timeout=120\n",
    "    )\n",
    "    response.raise_for_status()\n",
    "\n",
    "    sst_data = pd.read_csv(io.StringIO(response.text), skiprows=[1])\n",
    "\n",
    "    sst_data = sst_data.rename(columns={\n",
    "        \"time\": \"time (UTC)\",\n",
    "        \"latitude\": \"latitude (degrees_north)\",\n",
    "        \"longitude\": \"longitude (degrees_east)\",\n",
    "        \"analysed_sst\": \"analysed_sst (degree_C)\",\n",
    "        \"analysis_error\": \"analysis_error (kelvin)\"\n",
    "    })\n",
    "\n",
    "    if \"analysis_error (kelvin)\" not in sst_data.columns:\n",
    "        sst_data[\"analysis_error (kelvin)\"] = np.nan\n",
    "\n",
    "    sst_data[\"time (UTC)\"] = pd.to_datetime(\n",
    "        sst_data[\"time (UTC)\"],\n",
    "        errors=\"coerce\"\n",
    "    )\n",
    "\n",
    "    return sst_data\n",
    "\n",
    "\n",
    "def merge_adt_sst(adt_data, sst_data, coord_rounding=1):\n",
    "    \"\"\"\n",
    "    Merge ADT and SST by day and rounded coordinate.\n",
    "    \"\"\"\n",
    "    adt_for_merge = adt_data.copy()\n",
    "    sst_for_merge = sst_data.copy()\n",
    "\n",
    "    adt_for_merge[\"time_day\"] = pd.to_datetime(\n",
    "        adt_for_merge[\"time (UTC)\"],\n",
    "        errors=\"coerce\"\n",
    "    ).dt.normalize()\n",
    "\n",
    "    sst_for_merge[\"time_day\"] = pd.to_datetime(\n",
    "        sst_for_merge[\"time (UTC)\"],\n",
    "        errors=\"coerce\"\n",
    "    ).dt.normalize()\n",
    "\n",
    "    adt_for_merge[\"lat_round\"] = adt_for_merge[\n",
    "        \"latitude (degrees_north)\"\n",
    "    ].round(coord_rounding)\n",
    "\n",
    "    adt_for_merge[\"lon_round\"] = adt_for_merge[\n",
    "        \"longitude (degrees_east)\"\n",
    "    ].round(coord_rounding)\n",
    "\n",
    "    sst_for_merge[\"lat_round\"] = sst_for_merge[\n",
    "        \"latitude (degrees_north)\"\n",
    "    ].round(coord_rounding)\n",
    "\n",
    "    sst_for_merge[\"lon_round\"] = sst_for_merge[\n",
    "        \"longitude (degrees_east)\"\n",
    "    ].round(coord_rounding)\n",
    "\n",
    "    sst_for_merge[\"analysed_sst (degree_C)\"] = pd.to_numeric(\n",
    "        sst_for_merge[\"analysed_sst (degree_C)\"],\n",
    "        errors=\"coerce\"\n",
    "    )\n",
    "\n",
    "    sst_for_merge[\"analysis_error (kelvin)\"] = pd.to_numeric(\n",
    "        sst_for_merge[\"analysis_error (kelvin)\"],\n",
    "        errors=\"coerce\"\n",
    "    )\n",
    "\n",
    "    # Aggregate SST because several SST pixels may round to one ADT grid cell\n",
    "    sst_for_merge = (\n",
    "        sst_for_merge\n",
    "        .dropna(subset=[\n",
    "            \"time_day\",\n",
    "            \"lat_round\",\n",
    "            \"lon_round\",\n",
    "            \"analysed_sst (degree_C)\"\n",
    "        ])\n",
    "        .groupby([\"time_day\", \"lat_round\", \"lon_round\"], as_index=False)\n",
    "        .agg({\n",
    "            \"analysed_sst (degree_C)\": \"mean\",\n",
    "            \"analysis_error (kelvin)\": \"mean\"\n",
    "        })\n",
    "    )\n",
    "\n",
    "    merged = pd.merge(\n",
    "        adt_for_merge,\n",
    "        sst_for_merge,\n",
    "        on=[\"time_day\", \"lat_round\", \"lon_round\"],\n",
    "        how=\"left\"\n",
    "    )\n",
    "\n",
    "    merged = merged.drop(columns=[\"time_day\", \"lat_round\", \"lon_round\"])\n",
    "\n",
    "    merged = merged[\n",
    "        list(adt_data.columns) +\n",
    "        [\"analysed_sst (degree_C)\", \"analysis_error (kelvin)\"]\n",
    "    ]\n",
    "\n",
    "    return merged\n",
    "\n",
    "\n",
    "# ============================================================\n",
    "# MAIN LOOP: DOWNLOAD AND COMBINE ALL YEARS\n",
    "# ============================================================\n",
    "\n",
    "all_years = []\n",
    "\n",
    "for year in selected_years:\n",
    "    print(\"\\n\" + \"=\" * 60)\n",
    "    print(f\"Processing year: {year}\")\n",
    "\n",
    "    target_date, window_start, window_end = make_target_window(\n",
    "        year=year,\n",
    "        target_month_day=target_month_day,\n",
    "        days_before=days_before,\n",
    "        days_after=days_after\n",
    "    )\n",
    "\n",
    "    print(\"Target date:\", target_date.date())\n",
    "    print(\"Window:\", window_start.date(), \"to\", window_end.date())\n",
    "\n",
    "    adt_data = download_adt_for_year(\n",
    "        year=year,\n",
    "        window_start=window_start,\n",
    "        window_end=window_end,\n",
    "        lat_min=lat_min,\n",
    "        lat_max=lat_max,\n",
    "        lon_min=lon_min,\n",
    "        lon_max=lon_max\n",
    "    )\n",
    "\n",
    "    sst_data = download_sst_for_window(\n",
    "        window_start=window_start,\n",
    "        window_end=window_end,\n",
    "        lat_min=lat_min,\n",
    "        lat_max=lat_max,\n",
    "        lon_min=lon_min,\n",
    "        lon_max=lon_max,\n",
    "        time_stride=1\n",
    "    )\n",
    "\n",
    "    merged_year = merge_adt_sst(\n",
    "        adt_data=adt_data,\n",
    "        sst_data=sst_data,\n",
    "        coord_rounding=coord_rounding\n",
    "    )\n",
    "\n",
    "    merged_year[\"target_year\"] = year\n",
    "    merged_year[\"target_month_day\"] = target_month_day\n",
    "    merged_year[\"target_date\"] = target_date\n",
    "    merged_year[\"window_start\"] = window_start\n",
    "    merged_year[\"window_end\"] = window_end\n",
    "\n",
    "    merged_year[\"day_of_year\"] = merged_year[\"time (UTC)\"].dt.dayofyear\n",
    "    merged_year[\"month_day\"] = merged_year[\"time (UTC)\"].dt.strftime(\"%m-%d\")\n",
    "\n",
    "    print(\"Rows:\", len(merged_year))\n",
    "    print(\n",
    "        \"Rows with SST:\",\n",
    "        merged_year[\"analysed_sst (degree_C)\"].notna().sum()\n",
    "    )\n",
    "\n",
    "    all_years.append(merged_year)\n",
    "\n",
    "\n",
    "all_years_adt_sst = pd.concat(all_years, ignore_index=True)\n",
    "\n",
    "# Keep only the final columns and remove rows with missing ADT/SST values\n",
    "final_columns = [\n",
    "    \"time (UTC)\",\n",
    "    \"latitude (degrees_north)\",\n",
    "    \"longitude (degrees_east)\",\n",
    "    \"adt (m)\",\n",
    "    \"analysed_sst (degree_C)\",\n",
    "    \"analysis_error (kelvin)\",\n",
    "    \"target_year\"\n",
    "]\n",
    "\n",
    "all_years_adt_sst = (\n",
    "    all_years_adt_sst[final_columns]\n",
    "    .dropna()\n",
    "    .reset_index(drop=True)\n",
    ")\n",
    "\n",
    "print(\"Combined dataframe created: all_years_adt_sst\")\n",
    "print(\"Rows after dropping NaNs:\", len(all_years_adt_sst))\n",
    "print(\"Date range:\")\n",
    "print(\n",
    "    all_years_adt_sst[\"time (UTC)\"].min(),\n",
    "    \"to\",\n",
    "    all_years_adt_sst[\"time (UTC)\"].max()\n",
    ")\n",
    "\n",
    "all_years_adt_sst.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46",
   "metadata": {},
   "source": [
    "## Credits\n",
    "\n",
    "Data provided by NOAA CoastWatch, OceanWatch, and the AOML ERDDAP servers.\n",
    "\n",
    "This notebook was developed for the Rivulet Project to support climate and oceanographic inquiry. \n",
    "\n",
    "Simons, R.A. (2019) ERDDAP. https://coastwatch.pfeg.noaa.gov/erddap/index.html. NOAA/NMFS/SWFSC/ERD, Monterey, CA."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}