From 68142626baa4d3ad192978d38dee9bb32213521a Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Tue, 27 Sep 2022 11:54:10 -0600 Subject: [PATCH 1/5] Comparison - ASO - Update outlier value Get the .py scripts in line with notebook values --- comparison/ASO/ASO_difference.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/comparison/ASO/ASO_difference.py b/comparison/ASO/ASO_difference.py index 65591cc0..667dc9a7 100755 --- a/comparison/ASO/ASO_difference.py +++ b/comparison/ASO/ASO_difference.py @@ -20,7 +20,7 @@ ) # Values above 25m are considered outlier -OUTLIER = 25 +OUTLIER = 12 def violin_plot(ax, data, color): @@ -43,6 +43,7 @@ def area_plot(file, title='', save=False): data = differences.band_values() data[data > OUTLIER] = np.ma.masked + data[data < -OUTLIER] = np.ma.masked data[np.isnan(data)] = np.ma.masked fig, (ax1, ax2) = plt.subplots( From 0494080c660d74745a9d97bcedf4e098cbf580a6 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Tue, 27 Sep 2022 11:55:05 -0600 Subject: [PATCH 2/5] Comparison - ASO - Update differencing calculations Use reviewer suggested: bias = model - observation --- comparison/ASO/ASO_Snobal_difference.sh | 6 ++-- comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb | 35 +++++++++++---------- 2 files changed, 21 insertions(+), 20 deletions(-) diff --git a/comparison/ASO/ASO_Snobal_difference.sh b/comparison/ASO/ASO_Snobal_difference.sh index a5ae4a01..b8c492ae 100755 --- a/comparison/ASO/ASO_Snobal_difference.sh +++ b/comparison/ASO/ASO_Snobal_difference.sh @@ -29,7 +29,7 @@ echo $WATER_YEAR SNOBAL_50M=${DIFFERENCE_DIR}/Snobal_thickness_${DAY}_50m.tif gdal_translate --optfile $GDAL_OPTS \ - NETCDF:"${SNOBAL_HOME}/erw_isnobal/wy${WATER_YEAR}/erw/run${DAY}/snow.nc":thickness \ + NETCDF:"${SNOBAL_HOME}/erw_isnobal/GMD/wy${WATER_YEAR}/erw/run${DAY}/snow.nc":thickness \ ${SNOBAL_50M} if [ $? != 0 ]; then exit -1; fi @@ -54,6 +54,6 @@ rm ${DIFF_FILE} 2> /dev/null gdal_calc.py --co="TILED=YES" --co="COMPRESS=LZW" --co="NUM_THREADS=ALL_CPUS" \ --calc="A-B" \ - -A ${ASO_50M_SG} \ - -B ${SNOBAL_50M} \ + -A ${SNOBAL_50M} \ + -B ${ASO_50M_SG} \ --outfile ${DIFF_FILE} \ No newline at end of file diff --git a/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb b/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb index 82de2eac..8811d8ef 100644 --- a/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb +++ b/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb @@ -18,7 +18,8 @@ "metadata": {}, "outputs": [], "source": [ - "zone_data = cbrfc_zones()" + "zone_data = cbrfc_zones()\n", + "ASO_DIR = ASO_DIR / 'Depth-Difference/GMD'" ] }, { @@ -30,7 +31,7 @@ "source": [ "def plot_hv(first, second=None):\n", " diff = RasterFile(\n", - " (ASO_DIR / f'Depth-Difference/Depth_difference_{first}_50m.tif').as_posix()\n", + " (ASO_DIR / f'Depth_difference_{first}_50m.tif').as_posix()\n", " ).band_values()\n", "\n", " first = hv.Violin(diff[zone_data == ALEC2HLF].compressed(), label='ALEC2HLF').opts(\n", @@ -153,21 +154,21 @@ " model_cax.set_yticks([0,1])\n", " model_cax.set_yticklabels([])\n", " \n", - " model_cax.annotate(\"Under\",\n", + " model_cax.annotate(\"Over\",\n", " xy=(x_pos, 1),\n", " xytext=(x_pos, .84),\n", " xycoords='axes fraction',\n", " va=\"center\", ha=\"center\",\n", " rotation=270,\n", - " arrowprops=dict(arrowstyle=\"simple\", fc='red')\n", + " arrowprops=dict(arrowstyle=\"simple\", fc='blue')\n", " )\n", - " model_cax.annotate(\"Over\",\n", + " model_cax.annotate(\"Under\",\n", " xy=(x_pos, 0.00),\n", " xytext=(x_pos, .15),\n", " xycoords='axes fraction',\n", " va=\"center\", ha=\"center\",\n", " rotation=270,\n", - " arrowprops=dict(arrowstyle=\"simple\", fc='blue')\n", + " arrowprops=dict(arrowstyle=\"simple\", fc='red')\n", " )\n", "\n", "\n", @@ -208,10 +209,10 @@ " \n", "def plot_differences_mpl(axes, first, second):\n", " first_data = RasterFile(\n", - " (ASO_DIR / f'Depth-Difference/Depth_difference_{first}_50m.tif').as_posix()\n", + " (ASO_DIR / f'Depth_difference_{first}_50m.tif').as_posix()\n", " ).band_values()\n", " second_data = RasterFile(\n", - " (ASO_DIR / f'Depth-Difference/Depth_difference_{second}_50m.tif').as_posix()\n", + " (ASO_DIR / f'Depth_difference_{second}_50m.tif').as_posix()\n", " ).band_values()\n", " \n", " # Filter difference above 12m and below -12m\n", @@ -279,7 +280,7 @@ " mpatches.Patch(facecolor='indigo', ec='grey', label='Upper', alpha=0.5),\n", " mpatches.Patch(facecolor='brown', ec='grey', label='All', alpha=0.5),\n", " ],\n", - " loc='upper left',\n", + " loc='lower left',\n", " fontsize=8,\n", " frameon=False,\n", " title='HRU',\n", @@ -302,11 +303,11 @@ " set_ticks(axes, True)\n", " figure_label='b)'\n", "\n", - " axes.set_ylim(top=7, bottom=-3)\n", - " axes.set_yticks(np.arange(-2, 7, 2))\n", + " axes.set_ylim(top=3, bottom=-7)\n", + " axes.set_yticks(np.arange(-6, 3, 2))\n", "\n", " zero_line(axes)\n", - " axes.annotate(figure_label, xy=(0.02, 0.03 ), xycoords='axes fraction', fontsize=14)\n", + " axes.annotate(figure_label, xy=(0.02, 0.93 ), xycoords='axes fraction', fontsize=14)\n", " axes.set_title(datetime.strptime(first, '%Y%m%d').strftime('%Y'), size=14)\n", " \n", " print_stats(data, first, second)" @@ -349,7 +350,7 @@ "\n", " for day in [first, second]:\n", " diff = RasterFile(\n", - " (ASO_DIR / f'Depth-Difference/Depth_difference_{day}_50m.tif').as_posix()\n", + " (ASO_DIR / f'Depth_difference_{day}_50m.tif').as_posix()\n", " ).band_values()\n", "\n", " # Filter difference above 12m and below -12m\n", @@ -379,7 +380,7 @@ " zorder=1\n", " )\n", " \n", - " axes[ax_index].set_ylim(top=2, bottom=-1.5)\n", + " axes[ax_index].set_ylim(top=1.5, bottom=-2)\n", " axes[ax_index].set_xlabel('Aspect')\n", " \n", " if ax_index == 1:\n", @@ -443,7 +444,7 @@ "source": [ "def plot_area(day):\n", " differences = RasterFile(\n", - " (ASO_DIR / f'Depth-Difference/Depth_difference_{day}_50m.tif').as_posix()\n", + " (ASO_DIR / f'Depth_difference_{day}_50m.tif').as_posix()\n", " )\n", " \n", " data = differences.band_values()\n", @@ -485,7 +486,7 @@ " \n", " x_pos=1.03\n", "\n", - " ax1.annotate(\"Under\",\n", + " ax1.annotate(\"Over\",\n", " xy=(x_pos, 1),\n", " xytext=(x_pos, .85),\n", " xycoords='axes fraction',\n", @@ -493,7 +494,7 @@ " rotation=270,\n", " arrowprops=dict(arrowstyle=\"simple\", fc='white')\n", " )\n", - " ax1.annotate(\"Over\",\n", + " ax1.annotate(\"Under\",\n", " xy=(x_pos, 0.00),\n", " xytext=(x_pos, .15),\n", " xycoords='axes fraction',\n", From 52db717f8498517da9739e00d64c6cee48dfc725 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Tue, 27 Sep 2022 12:00:09 -0600 Subject: [PATCH 3/5] Comparison - Update figure sizes --- comparison/CBRFC/CBRFC-Precip.ipynb | 4 ++-- comparison/Snobal-SWI_static.ipynb | 12 ++++++++++-- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/comparison/CBRFC/CBRFC-Precip.ipynb b/comparison/CBRFC/CBRFC-Precip.ipynb index e8d5d45c..873fa12c 100644 --- a/comparison/CBRFC/CBRFC-Precip.ipynb +++ b/comparison/CBRFC/CBRFC-Precip.ipynb @@ -81,7 +81,7 @@ "outputs": [], "source": [ "wy_precip = xr.open_mfdataset(\n", - " (SNOBAL_DIR / 'wy20*' / 'erw/*/smrf_20*.nc').as_posix(),\n", + " (SNOBAL_DIR / 'GMD/wy20*' / 'erw/*/smrf_20*.nc').as_posix(),\n", " parallel=True,\n", ").resample(**RESAMPLE_1_DAY_OPTS).sum()\n", "\n", @@ -260,7 +260,7 @@ " return aspect_data\n", "\n", "def plot_aspect_mpl(precip, year):\n", - " fig, (ax) = plt.subplots(ncols=1, sharey=True, dpi=300, gridspec_kw={'wspace': 0})\n", + " fig, (ax) = plt.subplots(ncols=1, sharey=True, dpi=300, gridspec_kw={'wspace': 0}, figsize=(7,2.5))\n", "\n", " ax.boxplot(\n", " drop_nan(precip),\n", diff --git a/comparison/Snobal-SWI_static.ipynb b/comparison/Snobal-SWI_static.ipynb index 0033c216..28fdba30 100644 --- a/comparison/Snobal-SWI_static.ipynb +++ b/comparison/Snobal-SWI_static.ipynb @@ -102,7 +102,7 @@ "outputs": [], "source": [ "snobal_em = xr.open_mfdataset(\n", - " (SNOBAL_DIR / 'wy*/erw/*/em.nc').as_posix(),\n", + " (SNOBAL_DIR / 'GMD/wy*/erw/*/em.nc').as_posix(),\n", " data_vars=['SWI'],\n", " parallel=True,\n", ")\n", @@ -173,7 +173,7 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(dpi=300, figsize=(9,3.3))\n", + "plt.figure(dpi=300, figsize=(9,2.5))\n", "ax = plt.gca()\n", "\n", "swi.rolling(time=7, center=True).mean().plot(ax=ax, label='iSnobal Runoff', color='peru', alpha=0.7, lw=1)\n", @@ -199,6 +199,14 @@ "wy2020 = dict(time=slice(\"2019-10-01\", \"2020-09-30\"))\n", "wy2021 = dict(time=slice(\"2020-10-01\", \"2021-09-30\"))" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a25a860b-cde2-42bd-b673-6cbacdb07636", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { From 5390300c0b806fde374882c1af5822b6459ecf20 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Tue, 27 Sep 2022 12:01:15 -0600 Subject: [PATCH 4/5] Comparison - SNOTEL - Expand difference stats Expand comparison to SNOTEL with median difference and calculate percentage values. --- comparison/SNOTEL/Snotel-Snobal-stats.ipynb | 59 ++++++++++++++++----- 1 file changed, 45 insertions(+), 14 deletions(-) diff --git a/comparison/SNOTEL/Snotel-Snobal-stats.ipynb b/comparison/SNOTEL/Snotel-Snobal-stats.ipynb index db0871ef..38948969 100644 --- a/comparison/SNOTEL/Snotel-Snobal-stats.ipynb +++ b/comparison/SNOTEL/Snotel-Snobal-stats.ipynb @@ -33,7 +33,15 @@ "metadata": {}, "outputs": [], "source": [ - "year = 2021\n", + "year = 2018" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "water_year = f\"wy{year}\"\n", "time=slice(f\"{year -1}-10-01\", f\"{year}-07-31\")" ] @@ -88,7 +96,7 @@ "outputs": [], "source": [ "wy_snow = xr.open_mfdataset(\n", - " (SNOBAL_DIR / f'{water_year}' / 'erw/*/snow.nc').as_posix(),\n", + " (SNOBAL_DIR / f'GMD/{water_year}' / 'erw/*/snow.nc').as_posix(),\n", " parallel=True,\n", ")\n", "\n", @@ -314,7 +322,7 @@ "\n", "\n", "def legend_italic(ax, offset):\n", - " legend = ax.legend(bbox_to_anchor=(offset, 0.5), loc='center right')\n", + " legend = ax.legend(bbox_to_anchor=(0.5, offset), loc='lower center', borderaxespad=0.15, ncol=3)\n", "\n", " for text in legend.get_texts():\n", " if text.get_text() in plot_data.keys():\n", @@ -367,7 +375,7 @@ " print_hash(depth_numbers(site_name, max_snobal, min_snobal))\n", " aso_scatter(ax, site_name, max_snobal, min_snobal, mean_snobal)\n", " elif title == 'SWE':\n", - " print_hash(swe_numbers(site_name, max_snobal))\n", + " print_hash(swe_numbers(site_name, max_snobal, min_snobal))\n", "\n", " ax.set_title(None)\n", " ax.set_ylabel(None)\n", @@ -419,19 +427,22 @@ "outputs": [], "source": [ "def melt_date_snobal(mean_snobal, site):\n", - " after_max_swe = mean_snobal.sel(time=slice(f\"{year}-03-01\", f\"{year}-07-30\")).where(mean_snobal.time > plot_data[site]['max_swe_snotel'], drop=True)\n", + " after_max_swe = mean_snobal.sel(time=slice(f\"{year}-03-01\", f\"{year}-07-30\")).where(mean_snobal.time >= plot_data[site]['max_swe_snotel'], drop=True)\n", " return after_max_swe.where(after_max_swe != 0, drop=True).time.max().values\n", "\n", "\n", "def plot_depth_diff(site, ax1):\n", " mean_snobal = plot_data[site]['snobal']['thickness'].coarsen(**COARSEN_OPTS).mean().squeeze(['x', 'y'], drop=True)\n", " melt_date = melt_date_snobal(mean_snobal, site)\n", + "\n", + " depth_diff = mean_snobal.values - plot_data[site]['snotel']['Depth(m)']\n", + "\n", + " pre_melt = depth_diff[depth_diff.index <= plot_data[site]['max_swe_snotel']]\n", + " post_melt = depth_diff[(depth_diff.index > plot_data[site]['max_swe_snotel']) & (depth_diff.index <= melt_date)]\n", " \n", - " depth_diff = plot_data[site]['snotel']['Depth(m)'] - mean_snobal.values\n", + " mean_snobal_pre = mean_snobal.where(mean_snobal.time <= plot_data[site]['max_swe_snotel'], drop=True)\n", + " mean_snobal_post = mean_snobal.where((mean_snobal.time > plot_data[site]['max_swe_snotel']) & (mean_snobal.time <= melt_date), drop=True)\n", " \n", - " pre_melt = depth_diff[depth_diff.index < plot_data[site]['max_swe_snotel']]\n", - " post_melt = depth_diff[(depth_diff.index > plot_data[site]['max_swe_snotel']) & (depth_diff.index <= melt_date)]\n", - "\n", " ax1.plot([], [], ' ', label=site)\n", " depth_diff.plot(ax=ax1, \n", " color=plot_data[site_name]['color_snotel'], \n", @@ -442,10 +453,30 @@ " color=plot_data[site_name]['color_snobal'], \n", " label=f'Start melt season'\n", " )\n", + " ax1.axhline(y=0, ls='--', lw=0.5, alpha=0.5, color='grey')\n", " \n", " print(f\"{site}\")\n", - " print(\" Pre-melt mean: {:.2f}\".format(pre_melt.mean()))\n", - " print(\" Post mean: {:.2f}\".format(post_melt.mean()))" + " print(\" Melt date: \" + str(plot_data[site]['max_swe_snotel'].date()))\n", + " print(\" Mean\")\n", + " print(\" Difference\")\n", + " print(\" Pre-melt: {:.2f} m\".format(pre_melt.mean()))\n", + " print(\" Post : {:.2f} m\".format(post_melt.mean()))\n", + " print(\" Depth\")\n", + " print(\" Pre-melt: {:.2f} m\".format(mean_snobal_pre.mean().values))\n", + " print(\" Post : {:.2f} m\".format(mean_snobal_post.mean().values))\n", + " print(\" Difference Percent\")\n", + " print(\" Pre-melt: {:.2f} %\".format(pre_melt.mean() / mean_snobal_pre.mean().values * 100))\n", + " print(\" Post : {:.2f} %\".format(post_melt.mean() / mean_snobal_post.mean().values * 100))\n", + " print(\" Median\")\n", + " print(\" Difference\")\n", + " print(\" Pre-melt: {:.2f} m\".format(pre_melt.median()))\n", + " print(\" Post : {:.2f} m\".format(post_melt.median()))\n", + " print(\" Depth\")\n", + " print(\" Pre-melt: {:.2f} m\".format(mean_snobal_pre.median().values))\n", + " print(\" Post : {:.2f} m\".format(mean_snobal_post.median().values))\n", + " print(\" Difference Percent\")\n", + " print(\" Pre-melt: {:.2f} %\".format(pre_melt.median() / mean_snobal_pre.median().values * 100))\n", + " print(\" Post : {:.2f} %\".format(post_melt.median() / mean_snobal_post.median().values * 100))" ] }, { @@ -485,10 +516,10 @@ "ax1.set_ylabel(r'$\\Delta$ Snow Depth (m)')\n", "ax1.xaxis.set_major_formatter(xTicks)\n", "ax1.set_xlim([plot_range[0], plot_range[-1]])\n", - "ax1.set_ylim(-1.2, 1.0)\n", - "ax1.set_yticks(np.arange(-1.2, 1.2, 0.2))\n", + "ax1.set_ylim(-1.0, 1.2)\n", + "ax1.set_yticks(np.arange(-1, 1.3, 0.25))\n", "ax1.set_xlabel(f'Water Year {plot_range[-1].year}')\n", - "legend_italic(ax1, 1.3)" + "legend_italic(ax1, 1.05)" ] }, { From c292210a89732675e0a1676ee2d266887a78378c Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Thu, 6 Oct 2022 09:38:52 -0600 Subject: [PATCH 5/5] Comparison - ASO - Add violin plots for 2020 --- comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb | 102 +++++++++++++++++++- 1 file changed, 101 insertions(+), 1 deletion(-) diff --git a/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb b/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb index 8811d8ef..dd3f6c11 100644 --- a/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb +++ b/comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb @@ -207,7 +207,13 @@ " print(f\" Second: {second}\")\n", " single_stat(data[5])\n", " \n", - "def plot_differences_mpl(axes, first, second):\n", + " print(\"All\")\n", + " print(f\" First: {first}\")\n", + " single_stat(data[6])\n", + " print(f\" Second: {second}\")\n", + " single_stat(data[7])\n", + " \n", + "def plot_differences_mpl(axes, first, second=None):\n", " first_data = RasterFile(\n", " (ASO_DIR / f'Depth_difference_{first}_50m.tif').as_posix()\n", " ).band_values()\n", @@ -326,6 +332,100 @@ "plot_differences_mpl(ax2, '20190407', '20190610')" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d57ad6c-c18e-469e-8381-de6a44ee5d84", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(5.5,3.5), dpi=300)\n", + "ax1 = fig.gca()\n", + "\n", + "raster_data = RasterFile(\n", + " (ASO_DIR / f'Depth_difference_20200214_50m.tif').as_posix()\n", + ").band_values()\n", + "raster_data[np.isnan(raster_data)] = np.ma.masked\n", + "\n", + "# Filter difference above 12m and below -12m\n", + "raster_data[raster_data > 12] = np.ma.masked\n", + "raster_data[raster_data < -12] = np.ma.masked\n", + "\n", + "violin_position = np.arange(0.5, 2.5, 0.5)\n", + "data = [\n", + " raster_data[zone_data == ALEC2HLF].compressed(),\n", + " raster_data[zone_data == ALEC2HMF].compressed(),\n", + " raster_data[zone_data == ALEC2HUF].compressed(),\n", + " raster_data.compressed(),\n", + "]\n", + "\n", + "vp = ax1.violinplot(\n", + " data,\n", + " positions=violin_position,\n", + " widths=0.4,\n", + " showmeans=False,\n", + " showmedians=True,\n", + " showextrema=False,\n", + " quantiles=[\n", + " [0.05, 0.95],\n", + " [0.05, 0.95],\n", + " [0.05, 0.95],\n", + " [0.05, 0.95],\n", + " ],\n", + ")\n", + "\n", + "color = ['gold', 'teal', 'indigo', 'brown']\n", + "ci = 0\n", + "for pc in vp['bodies']:\n", + " pc.set_facecolor(color[ci])\n", + " pc.set_edgecolor('Black')\n", + " pc.set_lw(0.75)\n", + " ci += 1\n", + "\n", + "for line in ['cquantiles', 'cmedians']:\n", + " vp[line].set_color('black')\n", + " vp[line].set_lw(0.75)\n", + "\n", + "vp['cmedians'].set_ls(':')\n", + "vp['cmedians'].set_lw(1)\n", + "\n", + "ax1.set_xticks(\n", + " violin_position, \n", + " labels=['Lower', 'Middle', 'Upper', 'All']\n", + ")\n", + "ax1.set_xlabel('HRU')\n", + "set_ticks(ax1)\n", + "\n", + "ax1.legend(\n", + " handles=[\n", + " Line2D([0], [0], color='black', linestyle=':', label='Median'),\n", + " Line2D([0], [0], color='black', label='Quantiles (95%, 5%)')\n", + " ],\n", + " loc='upper right',\n", + " fontsize=8,\n", + " frameon=False,\n", + " bbox_to_anchor=(1.05, 1.15),\n", + ")\n", + "model_axes(ax1, 1.03)\n", + "set_ticks(ax1, True)\n", + "\n", + "ax1.set_ylim(top=2, bottom=-7)\n", + "ax1.set_yticks(np.arange(-6, 3, 2))\n", + "\n", + "zero_line(ax1)\n", + "# ax1.annotate(figure_label, xy=(0.02, 0.93 ), xycoords='axes fraction', fontsize=14)\n", + "ax1.set_title(datetime.strptime('20200214', '%Y%m%d').strftime('%Y-%m-%d'), size=12);\n", + "\n", + "print(\"HLF\")\n", + "single_stat(data[0])\n", + "print(\"HMF\")\n", + "single_stat(data[1])\n", + "print(\"HUF\")\n", + "single_stat(data[2])\n", + "print(\"All\")\n", + "single_stat(data[3])" + ] + }, { "cell_type": "code", "execution_count": null,