Skip to content

Commit

Permalink
Merge pull request #53 from UofU-Cryosphere/JH-paper
Browse files Browse the repository at this point in the history 
Updates with journal submission
  • Loading branch information
@jomey
jomey authored Sep 22, 2023
2 parents 30cb75f + aed0aa3 commit a962615
Show file tree
Hide file tree
Showing 5 changed files with 526 additions and 123 deletions.
322 changes: 307 additions & 15 deletions comparison/ASO/ASO-diff.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'"
]
},
{
Expand All @@ -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])"
]
},
{
Expand Down Expand Up @@ -70,7 +72,6 @@
" ASO_DIR / folder / f'Depth_difference_{DATES[1]}_50m.tif',\n",
" compressed\n",
" )\n",
"\n",
" \n",
" return data"
]
Expand Down Expand Up @@ -251,7 +252,9 @@
{
"cell_type": "markdown",
"id": "12c84783-0e81-45f6-8d2e-b25dd912e4a9",
"metadata": {},
"metadata": {
"tags": []
},
"source": [
"## Area Plot"
]
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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": {
Expand Down
Loading

0 comments on commit a962615

Please sign in to comment.