Feat: add superimpose debug app

wavey/__main__.py

@@ -20,9 +20,6 @@
 from wavey.map import DEFAULT_ARROW_LENGTH, RESOLUTION, Map
 from wavey.nwfs import download_forecast, get_most_recent_forecast
 
-# Force non-interactive backend to keep consistency between local and github actions
-matplotlib.rcParams["backend"] = "agg"
-
 LOG = logging.getLogger(__name__)
 
 # Location of Breakwater (aka San Carlos Beach)
@@ -99,6 +96,9 @@ def main(
         dither: Dithering increases image quality, but also increases storage size.
     """
 
+    # Force non-interactive backend to keep consistency between local and github actions
+    matplotlib.rcParams["backend"] = "agg"
+
     if resolution != "f":
         LOG.warning("Not drawing full resolution coastlines. Use the flag '--resolution f'")
 

wavey/apps/superimpose.py

@@ -0,0 +1,93 @@
+import datetime
+import logging
+from typing import Literal
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pygrib
+
+from wavey.__main__ import (
+    BREAKWATER_LAT_IDX,
+    BREAKWATER_LON_IDX,
+    DPI,
+    MONASTERY_LAT_IDX,
+    MONASTERY_LON_IDX,
+    utc_to_pt,
+)
+from wavey.common import FEET_PER_METER, setup_logging
+from wavey.grib import NUM_DATA_POINTS, ForecastData, ForecastType, read_forecast_data
+from wavey.nwfs import download_forecast, get_all_forecasts
+
+LOG = logging.getLogger(__name__)
+
+
+def main(
+    time: Literal["00", "06"] = "06",
+) -> None:
+    """
+    Superimpose forecasts from multiple times on the same graph. This is useful for visualizing how reliable forecasts
+    are with different forecast horizons.
+
+    Args:
+        time: Forecast time in UTC, either "00" or "06".
+    """
+
+    all_forecasts = get_all_forecasts(time=time)
+
+    # Download forecast data
+    wave_height_forecasts: list[ForecastData] = []
+    for forecast in all_forecasts:
+        grib_path = download_forecast(forecast)
+
+        with pygrib.open(grib_path) as grbs:
+            wave_height_forecasts.append(read_forecast_data(grbs, ForecastType.WaveHeight))
+
+    LOG.info("Plotting graph")
+    _, ax = plt.subplots(1, 1, figsize=(8, 4), dpi=DPI)
+
+    for i, wave_height_forecast in enumerate(wave_height_forecasts):
+        wave_height_ft = wave_height_forecast.data * FEET_PER_METER
+        analysis_date_pacific = utc_to_pt(wave_height_forecast.analysis_date_utc)
+
+        offset = i * 24  # each forecast is 24 hours into the past
+
+        bw_wave_height_ft = wave_height_ft[offset:, BREAKWATER_LAT_IDX, BREAKWATER_LON_IDX]
+        assert not np.ma.is_masked(bw_wave_height_ft), "Unexpected: Breakwater data contains masked points"
+
+        mon_wave_height_ft = wave_height_ft[offset:, MONASTERY_LAT_IDX, MONASTERY_LON_IDX]
+        assert not np.ma.is_masked(mon_wave_height_ft), "Unexpected: Monastery data contains masked points"
+
+        # NOTE: need to erase timezone info for mlpd3 to plot local times correctly
+        time0 = analysis_date_pacific.replace(tzinfo=None)
+        times = [time0 + datetime.timedelta(hours=hour_i) for hour_i in range(offset, NUM_DATA_POINTS)]
+
+        # plot styles
+        labels: tuple[str | None, str | None]
+        if offset == 0:
+            labels = ("Breakwater", "Monastery")
+            linestyle = "-"
+            alpha = 1.0
+        else:
+            labels = (None, None)
+            linestyle = "--"
+            alpha = 0.5 ** ((offset / 24) - 1)
+
+        ax.plot(times, bw_wave_height_ft, label=labels[0], color="blue", linestyle=linestyle, alpha=alpha)  # type: ignore[arg-type]
+        ax.plot(times, mon_wave_height_ft, label=labels[1], color="red", linestyle=linestyle, alpha=alpha)  # type: ignore[arg-type]
+
+    ax.set_ylim(0)
+    ax.set_ylabel("Significant wave height (ft)")
+    ax.yaxis.label.set_fontsize(14)
+    ax.xaxis.set_ticks_position("bottom")
+    ax.legend(loc="upper right")
+    ax.grid(True, linestyle=":", alpha=0.7)
+
+    plt.tight_layout()
+    plt.show()
+
+
+if __name__ == "__main__":
+    import tyro
+
+    setup_logging()
+    tyro.cli(main)

wavey/nwfs.py

@@ -23,6 +23,9 @@ def get_most_recent_forecast() -> str:
 
     Returns:
         URL to the GRIB file.
+
+    Raises:
+        HTTPError: If accessing website returns error.
     """
 
     # 1. List dates with forecasts
@@ -55,17 +58,7 @@ def get_most_recent_forecast() -> str:
         "Unexpected: could not find any forecasts for Monterey bay."
     )
 
-    # parse date and time
-    date_match = re.search(r"\d{8}", most_recent_date)
-    assert date_match, f"Unexpected date: {most_recent_date}"
-    time_match = re.search(r"\d{2}", most_recent_time)
-    assert time_match, f"Unexpected time: {most_recent_time}"
-
-    yyyymmdd = date_match.group(0)
-    hh = time_match.group(0)
-    filename = f"{_MTR}_nwps_{_CG3}_{yyyymmdd}_{hh}00.grib2"
-
-    url = os.path.join(_BASE_URL, most_recent_date, _MTR, most_recent_time, _CG3, filename)
+    url = _get_url(date=most_recent_date, time=most_recent_time)
     LOG.info(f"Found most recent forecast: {url}")
     return url
 
@@ -103,6 +96,9 @@ def _list_dates() -> list[str]:
 
     Returns:
         List of strings like "wr.YYYYMMDD/"; sorted (most recent first).
+
+    Raises:
+        HTTPError: If accessing website returns error.
     """
 
     url = _BASE_URL
@@ -121,7 +117,7 @@ def _list_times(date: str) -> list[str]:
         List of strings like "HH/"; sorted (most recent first).
 
     Raises:
-        HTTPError: If no forecasts for the given date.
+        HTTPError: If no forecasts for Monterey on the given date.
     """
 
     url = os.path.join(_BASE_URL, date, _MTR)
@@ -146,6 +142,54 @@ def _check_time(date: str, time: str) -> bool:
     return r.ok
 
 
+def _get_url(date: str, time: str) -> str:
+    """
+    Given date and time, get URL to the GRIB file.
+
+    Args:
+        date: A string like "wr.YYYYMMDD/".
+        time: A string like "HH/".
+
+    Returns:
+        URL to the GRIB file.
+    """
+
+    date_match = re.search(r"\d{8}", date)
+    assert date_match, f"Unexpected date: {date}"
+    time_match = re.search(r"\d{2}", time)
+    assert time_match, f"Unexpected time: {time}"
+
+    yyyymmdd = date_match.group(0)
+    hh = time_match.group(0)
+    filename = f"{_MTR}_nwps_{_CG3}_{yyyymmdd}_{hh}00.grib2"
+
+    return os.path.join(_BASE_URL, date, _MTR, time, _CG3, filename)
+
+
+def get_all_forecasts(time: str = "06") -> list[str]:
+    """
+    Get all NWFS forecast data for Monterey bay.
+
+    Args:
+        time: A string like "HH/".
+
+    Returns:
+        List of URLS to GRIB files.
+
+    Raises:
+        HTTPError: If accessing website returns error.
+    """
+
+    # 1. List dates with forecasts
+    dates = _list_dates()
+    LOG.info(f"Found NWFS forecasts: {dates}")
+
+    # 2. For each date, check for Monterey bay forecast on the given time
+    good_dates = [date for date in dates if _check_time(date=date, time=time)]
+
+    return [_get_url(date=date, time=time) for date in good_dates]
+
+
 def download_forecast(url: str, dir: Path | None = None) -> Path:
     """
     Download NWFS forecast data to disk.
@@ -157,6 +201,9 @@ def download_forecast(url: str, dir: Path | None = None) -> Path:
 
     Returns:
         Path to the GRIB file.
+
+    Raises:
+        HTTPError: If error encountered during download.
     """
 
     if dir is None: