Skip to content

Commit

Permalink
Merge pull request #48 from UofU-Cryosphere/GMD_revision_1
Browse files Browse the repository at this point in the history 
GMD submission revision updates
  • Loading branch information
@jomey
jomey authored Oct 6, 2022
2 parents b7c84cf + c292210 commit cfc9ee3
Show file tree
Hide file tree
Showing 6 changed files with 181 additions and 40 deletions.
6 changes: 3 additions & 3 deletions comparison/ASO/ASO_Snobal_difference.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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}
3 changes: 2 additions & 1 deletion comparison/ASO/ASO_difference.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
)

# Values above 25m are considered outlier
OUTLIER = 25
OUTLIER = 12


def violin_plot(ax, data, color):
Expand All @@ -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(
Expand Down
4 changes: 2 additions & 2 deletions comparison/CBRFC/CBRFC-Precip.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down
137 changes: 119 additions & 18 deletions comparison/CBRFC/CBRFC-zones-ASO-diff.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,8 @@
"metadata": {},
"outputs": [],
"source": [
"zone_data = cbrfc_zones()"
"zone_data = cbrfc_zones()\n",
"ASO_DIR = ASO_DIR / 'Depth-Difference/GMD'"
]
},
{
Expand All @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -206,12 +207,18 @@
" 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/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",
Expand Down Expand Up @@ -279,7 +286,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",
Expand All @@ -302,11 +309,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)"
Expand All @@ -325,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,
Expand All @@ -349,7 +450,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",
Expand Down Expand Up @@ -379,7 +480,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",
Expand Down Expand Up @@ -443,7 +544,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",
Expand Down Expand Up @@ -485,15 +586,15 @@
" \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",
" va=\"center\", ha=\"center\",\n",
" 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",
Expand Down
Loading

0 comments on commit cfc9ee3

Please sign in to comment.