From 744605b35e413c07e76c18a17b5f5f9a26172a0c Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Fri, 22 Sep 2023 11:09:16 -0600 Subject: [PATCH 1/4] SMESHR - Scripts - Use temporary file for extraction Breaks up the single piped command into two, using a temporary file. Reading from STDIN with GDAL was not reliably working. --- scripts/smeshr/dswrf_for_topo.sh | 7 ++++++- scripts/smeshr/tcdc_for_topo.sh | 9 +++++++-- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/scripts/smeshr/dswrf_for_topo.sh b/scripts/smeshr/dswrf_for_topo.sh index ded5a4a0..fed22f02 100755 --- a/scripts/smeshr/dswrf_for_topo.sh +++ b/scripts/smeshr/dswrf_for_topo.sh @@ -22,11 +22,16 @@ combine_day() { echo $DAY pushd $DAY > /dev/null - cat ${HRRR_PATTERN} | wgrib2 - -match "DSWRF:surface" -inv /dev/null -grib - | \ + TMPF=$DAY.grib2 + + cat ${HRRR_PATTERN} | wgrib2 - -match "DSWRF:surface" -inv /dev/null -grib $TMPF + hrrr_param_for_topo --topo "${TOPO_FILE}" \ + --hrrr_in $TMPF \ --nc_out "${NC_OUT_PREFIX}_${DAY}.dswrf.nc" \ --add-shading + rm $TMPF popd > /dev/null } diff --git a/scripts/smeshr/tcdc_for_topo.sh b/scripts/smeshr/tcdc_for_topo.sh index f8d1e924..23af563d 100755 --- a/scripts/smeshr/tcdc_for_topo.sh +++ b/scripts/smeshr/tcdc_for_topo.sh @@ -36,10 +36,15 @@ get_day() { echo $DAY pushd $DAY > /dev/null - cat ${HRRR_PATTERN} | wgrib2 - -match "TCDC:entire atmosphere" -inv /dev/null -grib - | \ + TMPF=$DAY.grib2 + + cat ${HRRR_PATTERN} | wgrib2 - -match "TCDC:entire atmosphere" -inv /dev/null -grib $TMPF + hrrr_param_for_topo --topo "${TOPO_FILE}" \ - --nc_out "${NC_OUT_PREFIX}_${DAY}.tcdc.nc" \ + --hrrr_in $TMPF \ + --nc_out "${NC_OUT_PREFIX}${DAY}.tcdc.nc" \ + rm $TMPF popd > /dev/null } From dca8fa52d27d89680d24179101b1e5fc62f0ac59 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Fri, 22 Sep 2023 11:19:32 -0600 Subject: [PATCH 2/4] Comparison - RasterFile - Improve handling of no data Improves handling of rasters that don't have a value for no-data set in the metadata. --- comparison/raster_file.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/comparison/raster_file.py b/comparison/raster_file.py index e89d7b81..13b56820 100644 --- a/comparison/raster_file.py +++ b/comparison/raster_file.py @@ -111,9 +111,10 @@ def band_values(self, **kwargs): band_number = kwargs.get('band_number', self.band_number) band = self.file.GetRasterBand(band_number) + no_data = band.GetNoDataValue() values = np.ma.masked_values( gdal_array.BandReadAsArray(band), - band.GetNoDataValue(), + no_data if no_data is not None else np.nan, copy=False ) From dcb24ad11dd82dfc73ea338048a992ffd3002901 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Fri, 22 Sep 2023 11:23:05 -0600 Subject: [PATCH 3/4] Comparison - ASO - Add figure wit no snow in model Adds figure that compares no snow in the model and measured depth with ASO --- comparison/ASO/ASO-diff.ipynb | 322 ++++++++++++++++++++++++++++++++-- 1 file changed, 307 insertions(+), 15 deletions(-) diff --git a/comparison/ASO/ASO-diff.ipynb b/comparison/ASO/ASO-diff.ipynb index 36c3d734..02ab60dc 100644 --- a/comparison/ASO/ASO-diff.ipynb +++ b/comparison/ASO/ASO-diff.ipynb @@ -26,7 +26,8 @@ "ASO_DIR = ASO_DIR / 'Depth-Difference/'\n", "DATES = ['20220421', '20220518']\n", "TIME_DECAY = 'SMRF-2022'\n", - "HRRR_MODIS = 'HRRR-MODIS'" + "HRRR_MODIS = 'HRRR-MODIS'\n", + "ASO_DEPTH = 'ASO-snow-depth'" ] }, { @@ -37,7 +38,8 @@ "outputs": [], "source": [ "ERW = gp.read_file(str(DATA_DIR / 'Boundaries/Upper_Gunnison.geojson'))\n", - "ERW_extent=(ERW.bounds.values.flatten()[0],ERW.bounds.values.flatten()[2], ERW.bounds.values.flatten()[1], ERW.bounds.values.flatten()[3])" + "ERW_bounds = ERW.bounds.values.flatten()\n", + "ERW_extent=(ERW_bounds[0], ERW_bounds[2], ERW_bounds[1], ERW_bounds[3])" ] }, { @@ -70,7 +72,6 @@ " ASO_DIR / folder / f'Depth_difference_{DATES[1]}_50m.tif',\n", " compressed\n", " )\n", - "\n", " \n", " return data" ] @@ -251,7 +252,9 @@ { "cell_type": "markdown", "id": "12c84783-0e81-45f6-8d2e-b25dd912e4a9", - "metadata": {}, + "metadata": { + "tags": [] + }, "source": [ "## Area Plot" ] @@ -379,15 +382,6 @@ "metadata": {}, "outputs": [], "source": [ - "def read_data(path, compressed):\n", - " data = RasterFile(path.as_posix()).band_values()\n", - " \n", - " # Filter difference above 12m and below -12m\n", - " data[data > 12] = np.ma.masked\n", - " data[data < -12] = np.ma.masked\n", - " \n", - " return data[~np.isnan(data)].compressed() if compressed else data\n", - "\n", "def read_days(compressed=False):\n", " data = {}\n", " \n", @@ -665,13 +659,311 @@ ");" ] }, + { + "cell_type": "markdown", + "id": "9b7c27a5-0b9e-4fe3-8832-67132dc5ac42", + "metadata": {}, + "source": [ + "## ASO > 0; HRRR-iSnobal == 0" + ] + }, + { + "cell_type": "markdown", + "id": "b80415e5-336a-4458-b70a-159a58e7d7aa", + "metadata": {}, + "source": [ + "### May" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7aeba9c-662c-413f-811f-61de1137db95", + "metadata": {}, + "outputs": [], + "source": [ + "ASO_SD_0518 = read_data(\n", + " ASO_DIR / HRRR_MODIS / f'ASO_50M_SD_USCOGE_{DATES[1]}_iSnobal_grid.tif',\n", + " False,\n", + ")\n", + "\n", + "HM_SD_0518 = read_data(\n", + " ASO_DIR / HRRR_MODIS / f'Snobal_thickness_{DATES[1]}_50m.tif',\n", + " False,\n", + ")\n", + "\n", + "TD_SD_0518 = read_data(\n", + " ASO_DIR / TIME_DECAY / f'Snobal_thickness_{DATES[1]}_50m.tif',\n", + " False,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d36ae61-8472-4a6b-9ffb-3a51ba89b194", + "metadata": {}, + "outputs": [], + "source": [ + "ASO_yes_model_no_HM = (ASO_SD_0518 > 0).data & (HM_SD_0518 == 0).data & ~np.isnan(aso_data['HRRR-MODIS']['20220421']).mask\n", + "ASO_yes_model_no_TD = (ASO_SD_0518 > 0).data & (TD_SD_0518 == 0).data & ~np.isnan(aso_data['SMRF-2022']['20220421']).mask" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "95a54d2e-33d3-4361-9e89-7e747f7c7c3c", + "metadata": {}, + "outputs": [], + "source": [ + "snow_aso_HM = np.ma.masked_where(~ASO_yes_model_no_HM, ASO_SD_0518)\n", + "snow_aso_TD = np.ma.masked_where(~ASO_yes_model_no_TD, ASO_SD_0518)" + ] + }, + { + "cell_type": "markdown", + "id": "0bc044ae-867c-450a-bb23-2c11dd5df0a2", + "metadata": {}, + "source": [ + "## Vegetation analysis " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea65b0b1-5962-444e-9717-09c598c21c5c", + "metadata": {}, + "outputs": [], + "source": [ + "topo_file = xr.open_dataset(\n", + " (DATA_DIR / 'project-data/iSnobal/ERW/topo/basin_setup/topo.nc').as_posix()\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e48d287f-3d06-4c26-9267-25d9e7fac97b", + "metadata": {}, + "outputs": [], + "source": [ + "veg_height_td_error = np.ma.masked_where(~ASO_yes_model_no_TD, topo_file.veg_height.values)\n", + "veg_height_hm_error = np.ma.masked_where(~ASO_yes_model_no_HM, topo_file.veg_height.values)\n", + "\n", + "td_values, td_count = np.unique(veg_height_td_error.compressed(), return_counts=True)\n", + "hm_values, hm_count = np.unique(veg_height_hm_error.compressed(), return_counts=True)\n", + "\n", + "td_count = np.append(td_count, 0) # Time-Decay has no values for 37.5 veg height" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5abdcade-e7d3-465f-8a15-8daef34069d8", + "metadata": {}, + "outputs": [], + "source": [ + "td_values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44d65aa1-6afd-4831-8663-d5f510ee609d", + "metadata": {}, + "outputs": [], + "source": [ + "hm_values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0edb00cc-a0b1-4f3d-966e-4a53e0538ea2", + "metadata": {}, + "outputs": [], + "source": [ + "hm_count" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "944b37e2-a602-4035-a64c-d84ca5a4b326", + "metadata": {}, + "outputs": [], + "source": [ + "td_count" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6b44a36d-578a-43a4-8e6e-c536b1b991ca", + "metadata": {}, + "outputs": [], + "source": [ + "veg_positions = np.arange(len(hm_values))\n", + "veg_labels = [\"%.2f\" % x for x in hm_values]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c3bd9be-54d4-45d8-95ab-9677e52a5cf1", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(dpi=300, figsize=(8, 4))\n", + "ax = fig.gca()\n", + "\n", + "ax.barh(\n", + " veg_positions - 0.15, \n", + " hm_count,\n", + " 0.3,\n", + " color='indigo', \n", + " ec='black',\n", + " alpha=0.2,\n", + " label=LABELS[1],\n", + ")\n", + "ax.barh(\n", + " veg_positions + 0.15, \n", + " td_count,\n", + " 0.3,\n", + " color='teal', \n", + " ec='black',\n", + " alpha=0.2,\n", + " label=LABELS[0],\n", + ")\n", + "ax.set_ylabel('Vegetation Height (m)')\n", + "ax.set_xlabel('Pixel Count')\n", + "ax.set_yticks(veg_positions)\n", + "ax.set_yticklabels(veg_labels)\n", + "ax.legend();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3c3af9a-29da-479b-b91c-55ce38481a98", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(dpi=300, figsize=(8, 4))\n", + "sub_figures = fig.subfigures(1, 2, width_ratios=[2, 1], wspace=0.0)\n", + "\n", + "grid = ImageGrid(\n", + " sub_figures[0], 111,\n", + " nrows_ncols=(1, 2),\n", + " axes_pad=(0.1, 0.3),\n", + " cbar_location=\"bottom\",\n", + " cbar_mode=\"single\",\n", + " cbar_pad=\"-4%\",\n", + " cbar_size='3%'\n", + ")\n", + "\n", + "plot_area(\n", + " snow_aso_TD, \n", + " grid[0], \n", + " title=LABELS[0], \n", + " ylabel='18 May'\n", + ")\n", + "ERW.boundary.plot(color='grey', ax=grid[0], lw=0.5)\n", + "\n", + "cax = grid.cbar_axes[0]\n", + "cbar = cax.colorbar(\n", + " plot_area(snow_aso_HM, grid[1], title=LABELS[1]),\n", + " ticks=np.arange(0, 1, 0.1),\n", + " extend='max'\n", + ")\n", + "cbar.set_label(\n", + " label=r'Snow Depth (m)',\n", + ")\n", + "ERW.boundary.plot(color='grey', ax=grid[1], lw=0.5)\n", + "\n", + "left_side = sub_figures[1].subplots(1,1)\n", + "\n", + "left_side.barh(\n", + " veg_positions - 0.15, \n", + " hm_count,\n", + " 0.3,\n", + " color='indigo', \n", + " ec='black',\n", + " alpha=0.2,\n", + " label=LABELS[1],\n", + ")\n", + "left_side.barh(\n", + " veg_positions + 0.15, \n", + " td_count,\n", + " 0.3,\n", + " color='teal', \n", + " ec='black',\n", + " alpha=0.2,\n", + " label=LABELS[0],\n", + ")\n", + "left_side.set_ylabel('Vegetation Height (m)')\n", + "left_side.set_xlabel('Pixel Count')\n", + "left_side.set_yticks(veg_positions)\n", + "left_side.set_yticklabels(veg_labels)\n", + "left_side.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "292c79f5-c89c-4122-b565-20dca979073e", + "metadata": {}, + "source": [ + "## Histogram" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4dc75df8-46f1-4cce-83a5-8918004042cb", + "metadata": {}, + "outputs": [], + "source": [ + "def print_stats(model_run, snow, pixel_count):\n", + " print(model_run)\n", + " print(f\" Median: {np.nanmedian(snow):.2} m\")\n", + " print(f\" Mean: {np.nanmean(snow):.2} m\")\n", + " print(f\" Std: {np.nanstd(snow):.2} m\")\n", + " print(f\" Max depth: {np.nanmax(snow):.2} m\")\n", + " print(\" Pixel count: \" + str(np.count_nonzero(snow)))\n", + " print(f\" Area Pixel %: {np.count_nonzero(snow)/pixel_count:.2%}\")" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "bce33a43-55dd-4bd0-9161-319633d072ee", + "id": "80ebc349-b5d1-4156-970d-9435e4f8103a", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "print_stats(LABELS[0], snow_aso_TD.compressed(), snow_aso_TD.size)\n", + "print_stats(LABELS[1], snow_aso_HM.compressed(), snow_aso_HM.size)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d47fec2-4b04-4e7b-9724-68ce70b7d73b", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(dpi=300, figsize=(8, 4))\n", + "\n", + "bins = np.arange(0, 2.1, 0.1)\n", + "ax = fig.gca()\n", + "ax.hist(snow_aso_HM.compressed(), bins, histtype='step', label=LABELS[1])\n", + "ax.hist(snow_aso_TD.compressed(), bins, histtype='step', label=LABELS[0])\n", + "\n", + "ax.set_ylabel('Pixel Count')\n", + "ax.set_xlabel('ASO Snow Depth (m)')\n", + "ax.set_xticks(np.arange(0, 1.2, 0.2))\n", + "ax.set_xlim(left=-0.02, right=1.2)\n", + "ax.legend();" + ] } ], "metadata": { From aed0aa35d8c3b54b4c2bc305a19f5d2336e224e6 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Fri, 22 Sep 2023 11:26:06 -0600 Subject: [PATCH 4/4] Comparison - Irwin albedo - Updates for journal submission. This is the version intially submitted with Journal of Hydrology. --- comparison/HRRR/Net_solar_SNOTEL-exact.ipynb | 308 ++++++++++++------- 1 file changed, 204 insertions(+), 104 deletions(-) diff --git a/comparison/HRRR/Net_solar_SNOTEL-exact.ipynb b/comparison/HRRR/Net_solar_SNOTEL-exact.ipynb index ba301b5c..b8e069dd 100644 --- a/comparison/HRRR/Net_solar_SNOTEL-exact.ipynb +++ b/comparison/HRRR/Net_solar_SNOTEL-exact.ipynb @@ -35,7 +35,7 @@ "metadata": {}, "outputs": [], "source": [ - "year = 2021\n", + "year = 2022\n", "water_year = f'wy{year}'" ] }, @@ -427,7 +427,8 @@ "import matplotlib.dates as mdates\n", "import matplotlib.patches as mpatches\n", "import matplotlib.ticker as mticker\n", - "from matplotlib.offsetbox import AnchoredText" + "from matplotlib.offsetbox import AnchoredText\n", + "from matplotlib.lines import Line2D" ] }, { @@ -437,14 +438,27 @@ "metadata": {}, "outputs": [], "source": [ - "def plot_site(data, ax, site_name, y_limit, add_label=False):\n", + "def plot_site(data, ax, site_name, y_limit, add_label=False, add_mean=True):\n", " print(site_name)\n", " \n", + " if add_mean:\n", + " mean_ax = ax.twinx()\n", + " else:\n", + " mean_ax = None\n", + " \n", " for key in data: \n", " values = np.nan_to_num(data[key].data.flatten(), copy=False)\n", " \n", - " print(f' {key} Mean: {np.nanmean(values):.2f}')\n", + " mean = np.nanmean(values)\n", + " print(f' {key} Mean: {mean:.2f}')\n", " \n", + " if add_mean:\n", + " mean_ax.axhline(\n", + " mean, \n", + " ls=(0, (1, 1)), lw=1, \n", + " color=COLORS[key], alpha=0.9, \n", + " )\n", + "\n", " ax.plot(\n", " data[key].time,\n", " values,\n", @@ -452,6 +466,9 @@ " color=COLORS[key], \n", " alpha=0.9, lw=1\n", " )\n", + "\n", + " if add_mean:\n", + " ax.plot([], [], color='dimgrey', ls=(0, (1, 1)), lw=1, label='Season Means')\n", " \n", " ax.set_ylabel(r'Net Solar ($W/m^2$)')\n", "\n", @@ -462,8 +479,15 @@ " ax.set_ylim(bottom=-0.1, top=y_limit)\n", " ax.yaxis.set_minor_locator(mticker.MultipleLocator(10))\n", " \n", + " if add_mean:\n", + " mean_max = 100 if year == 2021 else 140\n", + " mean_ax.set_ylim(bottom=0, top=mean_max)\n", + " mean_ax.set_ylabel(r'Mean Net Solar ($W/m^2$)')\n", + " \n", " if add_label:\n", " site_label(site_name, ax)\n", + " \n", + " return mean_ax\n", "\n", "def site_label(site_name, ax):\n", " at = AnchoredText(\n", @@ -494,14 +518,6 @@ " axes.xaxis.set_major_formatter(mticker.NullFormatter())\n", " axes.xaxis.set_minor_formatter(xTicks)\n", " axes.tick_params(axis='x', which='minor', pad=-0.5)\n", - " \n", - " axes.legend(\n", - " frameon=False,\n", - " bbox_to_anchor=(0.75, -0.125),\n", - " ncol=4,\n", - " borderaxespad=0.15, \n", - " fontsize=8\n", - " )\n", "\n", "def set_title(axes):\n", " axes.add_artist(\n", @@ -516,101 +532,23 @@ " )" ] }, - { - "cell_type": "markdown", - "id": "36c08a21-91a5-4ebd-b351-20b1bfd2b0be", - "metadata": {}, - "source": [ - "## Net Solar " - ] - }, { "cell_type": "code", "execution_count": null, - "id": "cdbb19ea-9619-4671-8559-584fde590d27", + "id": "27d51346-8af0-41ca-b1ec-f142517b611c", "metadata": {}, "outputs": [], "source": [ - "figure_opts = dict(figsize=(7,6), dpi=300,)\n", "plot_range = pd.date_range(start=SNOW_ONSET[str(year)]['Schofield Pass'], periods=10, freq='MS')\n", - "xTicks = mdates.DateFormatter('%b')\n", - "\n", - "fig, axes = plt.subplots(2, 1, sharex=True, **figure_opts)\n", - "plt.subplots_adjust(hspace=0.1)\n", - "\n", - "MAX_NS = 400 if year == 2021 else 450\n", - "\n", - "plot_site(\n", - " {\n", - " 'HRRR-MODIS': with_snow_at_site('Butte', butte_snobal_hrrr_modis),\n", - " 'HRRR-SC': with_snow_at_site('Butte', butte_snobal_hrrr),\n", - " 'Time-Decay': with_snow_at_site('Butte', butte_snobal),\n", - " },\n", - " axes[0],\n", - " 'Butte',\n", - " MAX_NS,\n", - " True\n", - ")\n", - "\n", - "plot_site(\n", - " {\n", - " 'HRRR-MODIS': with_snow_at_site('Schofield Pass', schofield_snobal_hrrr_modis),\n", - " 'HRRR-SC': with_snow_at_site('Schofield Pass', schofield_snobal_hrrr),\n", - " 'Time-Decay': with_snow_at_site('Schofield Pass', schofield_snobal),\n", - " },\n", - " axes[1],\n", - " 'Schofield Pass',\n", - " MAX_NS,\n", - " True\n", - ")\n", - "\n", - "style_axes(axes[1])\n", - "set_title(axes[0])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3d8a53d6-0564-4884-a277-70462f23a3cd", - "metadata": {}, - "outputs": [], - "source": [ - "figure_opts = dict(figsize=(7,3.5), dpi=300,)\n", - "fig, axes = plt.subplots(1, 1, sharex=True, **figure_opts)\n", - "\n", - "MAX_NS = 640 if year == 2021 else 620\n", - "\n", - "plot_site(\n", - " {\n", - " 'HRRR-MODIS': with_snow_at_site('Irwin Study Plot', irwin_snobal_hrrr_modis),\n", - " 'HRRR-SC': with_snow_at_site('Irwin Study Plot', irwin_snobal_hrrr),\n", - " 'Time-Decay': with_snow_at_site('Irwin Study Plot', irwin_snobal),\n", - " },\n", - " axes,\n", - " 'Irwin Study Plot',\n", - " MAX_NS,\n", - ")\n", - "\n", - "irwin_snow = df_with_snow_at_site('Irwin Study Plot', irwin_MODIS)\n", - "\n", - "axes.plot(\n", - " irwin_snow.index,\n", - " irwin_snow['net_solar'],\n", - " label=STATION_LABEL, \n", - " color=COLORS[STATION_LABEL], \n", - " alpha=0.9, lw=1\n", - ")\n", - "\n", - "print(f\" Station Mean: {irwin_snow['net_solar'].mean():.2f}\")\n", - "\n", - "style_axes(axes)\n", - "subplot_letter(('e' if year == 2021 else 'f'), axes)" + "xTicks = mdates.DateFormatter('%b')" ] }, { "cell_type": "markdown", "id": "703745bb-b7ac-4b98-9e9a-97c00d1a7916", - "metadata": {}, + "metadata": { + "tags": [] + }, "source": [ "## Incoming Solar " ] @@ -637,16 +575,17 @@ "figure_opts = dict(figsize=(7,3.5), dpi=300,)\n", "fig, axes = plt.subplots(1, 1, sharex=True, **figure_opts)\n", "\n", - "MAX_SUN = 1600\n", + "MAX_SUN = 1200\n", "\n", - "plot_site(\n", + "mean_ax = plot_site(\n", " {\n", " 'HRRR-SC': with_snow_at_site('Irwin Study Plot', irwin_snobal_hrrr_solar),\n", - " 'Time-Decay': with_snow_at_site('Irwin Study Plot', irwin_snobal_solar),\n", + " # 'Time-Decay': with_snow_at_site('Irwin Study Plot', irwin_snobal_solar),\n", " },\n", " axes,\n", " 'Irwin Study Plot',\n", " MAX_SUN,\n", + " add_mean=True\n", ")\n", "\n", "irwin_snow = df_with_snow_at_site('Irwin Study Plot', irwin_MODIS)\n", @@ -658,6 +597,12 @@ " color=COLORS[STATION_LABEL], \n", " alpha=0.9, lw=1\n", ")\n", + "mean_ax.axhline(\n", + " irwin_snow['solar'].mean(), \n", + " ls=(0, (1, 1)), lw=1, \n", + " color=COLORS[STATION_LABEL], alpha=0.9, \n", + ")\n", + "mean_ax.set_ylim(bottom=0, top=400)\n", "\n", "print(f\" Station Mean: {irwin_snow['solar'].mean():.2f}\")\n", "\n", @@ -668,13 +613,22 @@ "\n", "style_axes(axes)\n", "set_title(axes)\n", + "axes.legend(\n", + " frameon=False,\n", + " bbox_to_anchor=(0.78, -0.125),\n", + " ncol=3,\n", + " borderaxespad=0.15, \n", + " fontsize=8\n", + ")\n", "subplot_letter(('a' if year == 2021 else 'b'), axes)" ] }, { "cell_type": "markdown", "id": "a981fd7a-049f-4cbd-8eab-2787ca79a3d6", - "metadata": {}, + "metadata": { + "tags": [] + }, "source": [ "## Albedo " ] @@ -733,6 +687,118 @@ "subplot_letter(('c' if year == 2021 else 'd'), ax)" ] }, + { + "cell_type": "markdown", + "id": "36c08a21-91a5-4ebd-b351-20b1bfd2b0be", + "metadata": { + "tags": [] + }, + "source": [ + "## Net Solar " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdbb19ea-9619-4671-8559-584fde590d27", + "metadata": {}, + "outputs": [], + "source": [ + "figure_opts = dict(figsize=(7,6), dpi=300,)\n", + "\n", + "fig, axes = plt.subplots(2, 1, sharex=True, **figure_opts)\n", + "plt.subplots_adjust(hspace=0.1)\n", + "\n", + "MAX_NS = 400 if year == 2021 else 450\n", + "\n", + "plot_site(\n", + " {\n", + " 'HRRR-MODIS': with_snow_at_site('Butte', butte_snobal_hrrr_modis),\n", + " 'HRRR-SC': with_snow_at_site('Butte', butte_snobal_hrrr),\n", + " 'Time-Decay': with_snow_at_site('Butte', butte_snobal),\n", + " },\n", + " axes[0],\n", + " 'Butte',\n", + " MAX_NS,\n", + " True\n", + ")\n", + "\n", + "plot_site(\n", + " {\n", + " 'HRRR-MODIS': with_snow_at_site('Schofield Pass', schofield_snobal_hrrr_modis),\n", + " 'HRRR-SC': with_snow_at_site('Schofield Pass', schofield_snobal_hrrr),\n", + " 'Time-Decay': with_snow_at_site('Schofield Pass', schofield_snobal),\n", + " },\n", + " axes[1],\n", + " 'Schofield Pass',\n", + " MAX_NS,\n", + " True\n", + ")\n", + "\n", + "style_axes(axes[1])\n", + "axes[1].legend(\n", + " frameon=False,\n", + " bbox_to_anchor=(0.9, -0.125),\n", + " ncol=4,\n", + " borderaxespad=0.15, \n", + " fontsize=8\n", + ")\n", + "set_title(axes[0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3d8a53d6-0564-4884-a277-70462f23a3cd", + "metadata": {}, + "outputs": [], + "source": [ + "figure_opts = dict(figsize=(7,3.5), dpi=300,)\n", + "fig, axes = plt.subplots(1, 1, sharex=True, **figure_opts)\n", + "\n", + "MAX_NS = 640 if year == 2021 else 620\n", + "\n", + "mean_ax = plot_site(\n", + " {\n", + " 'HRRR-MODIS': with_snow_at_site('Irwin Study Plot', irwin_snobal_hrrr_modis),\n", + " 'HRRR-SC': with_snow_at_site('Irwin Study Plot', irwin_snobal_hrrr),\n", + " # 'Time-Decay': with_snow_at_site('Irwin Study Plot', irwin_snobal),\n", + " },\n", + " axes,\n", + " 'Irwin Study Plot',\n", + " MAX_NS,\n", + " add_mean=True,\n", + ")\n", + "\n", + "irwin_snow = df_with_snow_at_site('Irwin Study Plot', irwin_MODIS)\n", + "\n", + "axes.plot(\n", + " irwin_snow.index,\n", + " irwin_snow['net_solar'],\n", + " label=STATION_LABEL, \n", + " color=COLORS[STATION_LABEL], \n", + " alpha=0.9, lw=1\n", + ")\n", + "mean_ax.axhline(\n", + " irwin_snow['net_solar'].mean(), \n", + " ls=(0, (1, 1)), lw=1, \n", + " color=COLORS[STATION_LABEL], alpha=0.9, \n", + ")\n", + "mean_ax.set_ylim(bottom=0, top=140)\n", + "\n", + "print(f\" Station Mean: {irwin_snow['net_solar'].mean():.2f}\")\n", + "\n", + "style_axes(axes)\n", + "axes.legend(\n", + " frameon=False,\n", + " bbox_to_anchor=(0.88, -0.125),\n", + " ncol=4,\n", + " borderaxespad=0.15, \n", + " fontsize=8\n", + ")\n", + "subplot_letter(('e' if year == 2021 else 'f'), axes)" + ] + }, { "cell_type": "markdown", "id": "43c5fdb7-60ba-4bcc-9b56-1dd1dc6acb2b", @@ -755,7 +821,7 @@ "data = {\n", " 'HRRR-MODIS': df_with_snow_at_site('Irwin Study Plot', hrrr_modis_diff),\n", " 'HRRR-SC': df_with_snow_at_site('Irwin Study Plot', hrrr_sc_diff),\n", - " 'Time-Decay': df_with_snow_at_site('Irwin Study Plot', time_decay_diff),\n", + " # 'Time-Decay': df_with_snow_at_site('Irwin Study Plot', time_decay_diff),\n", "}" ] }, @@ -790,6 +856,13 @@ ")\n", "\n", "style_axes(axes)\n", + "axes.legend(\n", + " frameon=False,\n", + " bbox_to_anchor=(0.78, -0.125),\n", + " ncol=3,\n", + " borderaxespad=0.15, \n", + " fontsize=8\n", + ")\n", "axes.set_ylabel(r'Δ Net Solar ($W/m^2$)')\n", "axes.set_ylim(bottom=-200, top=500)\n", "axes.set_yticks(np.arange(-200, 500+1, 50))\n", @@ -839,6 +912,27 @@ "x_lim = (np.datetime64(f'{year-1}-10-01'), np.datetime64(f'{year}-06-01'))" ] }, + { + "cell_type": "markdown", + "id": "cd253c72-36ce-4336-9951-19eaa2330d7d", + "metadata": {}, + "source": [ + "Fix pandas 0.14.6 error for time series with no name\n", + "Error message: `DataError: Having a non-string as a column name in a DataFrame is not supported.`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a931b7bd-5fac-4f63-a1fe-43fdbb4b8d33", + "metadata": {}, + "outputs": [], + "source": [ + "hrrr_modis_diff.name = 'HRRR MODIS Difference'\n", + "hrrr_sc_diff.name = 'HRRR-SC Difference'\n", + "time_decay_diff.name = 'Time-Decay Difference'" + ] + }, { "cell_type": "code", "execution_count": null, @@ -952,14 +1046,20 @@ "metadata": {}, "outputs": [], "source": [ - "(\n", + "vis = (\n", " irwin_pyra['VIS_IN'].rolling(ROLL_WINDOW).mean().interpolate('time') - \n", " irwin_pyra['VIS_OUT'].rolling(ROLL_WINDOW).mean().interpolate('time')\n", - ").hvplot(label='VIS', **HV_PLOT_OPTS) * \\\n", - "(\n", + ")\n", + "vis.name = 'VIS Difference'\n", + "\n", + "nir = (\n", " irwin_pyra['NIR_IN'].rolling(ROLL_WINDOW).mean().interpolate('time') - \n", " irwin_pyra['NIR_OUT'].rolling(ROLL_WINDOW).mean().interpolate('time')\n", - ").hvplot(title='Irwin Study Plot', label='NIR', **HV_PLOT_OPTS)" + ")\n", + "nir.name = 'NIR Difference'\n", + "\n", + "vis.hvplot(label='VIS', **HV_PLOT_OPTS) * \\\n", + "nir.hvplot(title='Irwin Study Plot', label='NIR', **HV_PLOT_OPTS)" ] }, {