yield_seasonality.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5e804417",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For the calculation\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "\n",
    "# For the animation plot\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "\n",
    "# Data at:\n",
    "#   gs://ml-training-zarrs/pv/pvlive_gsp_new_boundaries_2019-2025.zarr\n",
    "# but also on donatello at the path used below\n",
    "ds_gen = (\n",
    "    xr.open_zarr(\"/mnt/storage_u2_30tb_a/ml_training_zarrs/pv/pvlive_gsp_new_boundaries_2019-2025.zarr\")\n",
    "    .sel(gsp_id=0, drop=True)\n",
    "    .compute()\n",
    ")\n",
    "\n",
    "da_yield = (ds_gen.generation_mw / ds_gen.capacity_mwp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e2d144e7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 366/366 [00:22<00:00, 16.28it/s]\n"
     ]
    }
   ],
   "source": [
    "def slice_to_dayofyear_window(\n",
    "    da: xr.DataArray, \n",
    "    central_day_num: int, \n",
    "    window_size: int\n",
    ") -> xr.DataArray:\n",
    "    \"\"\"Slice the xarray object to timestamps which are within N days of the central day\n",
    "\n",
    "    Args:\n",
    "        da: The xarray object (must have datetime_gmt dimension)\n",
    "        central_day: The day of year of the central day in the window\n",
    "        window_size: The window size\n",
    "    \"\"\"\n",
    "\n",
    "    daysofyear = da.datetime_gmt.dt.dayofyear\n",
    "\n",
    "    # Find the window of days\n",
    "    # If the window size is even, then there will be one more day to the left of the central day\n",
    "    # than to to the right.\n",
    "    min_day = central_day_num - window_size//2\n",
    "    max_day = central_day_num + (window_size-1)//2\n",
    "    \n",
    "    mask = (min_day <= daysofyear) & (daysofyear <= max_day)\n",
    "\n",
    "    # In January we may need to include days from December\n",
    "    if min_day < 1:\n",
    "        mask = mask | (daysofyear >= 366 + min_day)\n",
    "    \n",
    "    # In December we may need to include days from January\n",
    "    if max_day > 366:\n",
    "        mask = mask | (daysofyear <= (max_day - 366))\n",
    "    \n",
    "    return da.where(mask, drop=True)\n",
    "\n",
    "\n",
    "def windowed_dayofyear_quantiles(\n",
    "    da: xr.DataArray, \n",
    "    window_size: int,\n",
    "    quantiles: list[float]\n",
    ") -> xr.DataArray:\n",
    "    \"\"\"Calculate quantiles of input data for each day-of-year and time\n",
    "\n",
    "    Args:\n",
    "        da: The xarray object (must have datetime_gmt dimension)\n",
    "        window_size: The window size\n",
    "        quantiles: The quantiles to calculate\n",
    "    \"\"\"\n",
    "    \n",
    "    unique_days_of_year = np.unique(da.datetime_gmt.dt.dayofyear.values)\n",
    "\n",
    "    results = []\n",
    "    for day_num in tqdm(unique_days_of_year):\n",
    "        # Slice data to values within date within window of days\n",
    "        da_day = slice_to_dayofyear_window(da, day_num, window_size)\n",
    "        \n",
    "        da_res = (\n",
    "            da_day.groupby(da_day.datetime_gmt.dt.time)\n",
    "            .quantile(quantiles, dim=\"datetime_gmt\")\n",
    "            .expand_dims(dayofyear=[day_num])\n",
    "        )\n",
    "        results.append(da_res)\n",
    "\n",
    "    return xr.concat(results, dim=\"dayofyear\")\n",
    "\n",
    "\n",
    "da_res = windowed_dayofyear_quantiles(\n",
    "    da_yield,\n",
    "    window_size=30,\n",
    "    quantiles=[0.1, 0.25, 0.5, 0.75, 0.9]\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "32648a23",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_animation(\n",
    "    da: xr.DataArray, \n",
    "    figsize: tuple[int, int] = (6,6),\n",
    ") -> FuncAnimation:\n",
    "    \"\"\"Make animation of yield seasonality\"\"\"\n",
    "\n",
    "    # Set up the figure\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "\n",
    "    # Convert the times to float hour-of-day\n",
    "    x = [(t.hour + t.minute/60) for t in da.time.values]\n",
    "\n",
    "    # Initialize lines (one for each quantile)\n",
    "    lines = []\n",
    "    for q in da[\"quantile\"].values:\n",
    "        y = da.isel(dayofyear=0).sel(quantile=q).values\n",
    "        line, = ax.plot(x, y, label=f'q={q}')\n",
    "        lines.append(line)\n",
    "\n",
    "    # Labels, legends, and titles\n",
    "    ax.set_xlabel(\"Time\")\n",
    "    ax.set_ylabel(\"Hour of day\")\n",
    "    ax.set_ylim(0, da.max()*1.2)\n",
    "    ax.legend(loc=\"upper right\")\n",
    "    title = ax.set_title(\"Day: 1\")\n",
    "\n",
    "    def update(frame_idx: int):\n",
    "\n",
    "        # Select the day of year\n",
    "        daily_slice = da.isel(dayofyear=frame_idx)\n",
    "        \n",
    "        # Update all the quantile lines\n",
    "        for q, line in zip(da[\"quantile\"].values, lines):\n",
    "            y = daily_slice.sel(quantile=q).values\n",
    "            line.set_data(x, y)\n",
    "        \n",
    "        # Update the title\n",
    "        current_doy = da.dayofyear[frame_idx].item()\n",
    "        title.set_text(f\"Day of Year: {current_doy}\")\n",
    "        \n",
    "        return lines + [title]\n",
    "\n",
    "    anim = FuncAnimation(\n",
    "        fig, \n",
    "        update, \n",
    "        frames=len(da.dayofyear), \n",
    "        interval=50, \n",
    "        blit=True\n",
    "    )\n",
    "    plt.close()\n",
    "\n",
    "    return anim\n",
    "\n",
    "anim = make_animation(da_res)\n",
    "\n",
    "# Show in jupyter notebook\n",
    "# HTML(anim.to_jshtml())\n",
    "\n",
    "# Save\n",
    "anim.save(\"yield_seasonality.gif\", writer=\"pillow\", fps=30)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "james-pvnet",
   "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.11.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}