Skip to content

Commit

Permalink
add current wet snow test material
Browse files Browse the repository at this point in the history
  • Loading branch information
@egagli
egagli committed Sep 6, 2023
1 parent 2a91691 commit 8258cf9
Show file tree
Hide file tree
Showing 3 changed files with 1,078 additions and 0 deletions.
299 changes: 299 additions & 0 deletions contrib/gagliano/wet_snow_get_rmse_all_sites.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,299 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import numpy as np\n",
"from pathlib import Path\n",
"from glob import glob\n",
"from datetime import datetime\n",
"import sys\n",
"sys.path.append('../../../spicy-snow/')\n",
"\n",
"from spicy_snow.processing.snow_index import calc_delta_cross_ratio, calc_delta_gamma, \\\n",
" clip_delta_gamma_outlier, calc_snow_index, calc_snow_index_to_snow_depth\n",
"from spicy_snow.processing.wet_snow import id_newly_wet_snow, id_wet_negative_si, \\\n",
" id_newly_frozen_snow, flag_wet_snow\n",
"from spicy_snow.retrieval import retrieval_from_parameters\n",
"\n",
"from dask.distributed import Client\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"client = Client(local_directory='/tmp', processes=False)\n",
"client"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"files = sorted(glob('spicy_s1_stacks/*.nc'))\n",
"\n",
"\n",
"f = files[1]\n",
"ds_name = f.split('stacks/')[-1].split('.')[0]\n",
"print(datetime.now(), f' -- starting {ds_name}')\n",
"\n",
" # Open dataset \n",
"ds_ = xr.open_dataset(f).load()\n",
"dataset = ds_[['s1','deltaVV','ims','fcf','lidar-sd']]\n",
"td = abs(pd.to_datetime(dataset.time) - pd.to_datetime(dataset.attrs['lidar-flight-time']))\n",
"closest_ts_idx = np.where(td == td.min())[0][0]\n",
"closest_ts = dataset.time[closest_ts_idx]\n",
"\n",
"a = 2.5\n",
"b = 0.2\n",
"c = 0.55"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds = retrieval_from_parameters(dataset,A=a,B=b,C=c,wet_SI_thresh=2,freezing_snow_thresh=1,wet_snow_thres=-2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds['wet_snow'].plot(col='time',col_wrap=10)\n",
"#ds['wet_flag'].plot(col='time',col_wrap=10)\n",
"#ds['alt_wet_flag'].plot(col='time',col_wrap=10)\n",
"#ds['freeze_flag'].plot(col='time',col_wrap=10)\n",
"#ds['perma_wet'].plot(col='time',col_wrap=10)\n",
"\n",
"#ds['snow_index'].plot(col='time',col_wrap=10)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mask_wet = ~(ds['lidar-sd'].isnull() | ds['snow_depth'].sel(time=closest_ts).isnull() | ds['wet_snow'].sel(time=closest_ts).astype(bool))\n",
"diff_wet = ds['lidar-sd'].where(mask_wet) - ds['snow_depth'].sel(time=closest_ts).where(mask_wet)\n",
"rmse_wet = float(np.sqrt((diff_wet**2).sum()/len(diff_wet.values.flatten())))\n",
"print(f'RMSE with wet snow masked out = {rmse_wet:0.2f}')\n",
"#rmse_wet_flag.loc[a, b, c,wst] = rmse_wet\n",
"# Compare snow depths - no wet snow mask\n",
"mask = ~(ds['lidar-sd'].isnull() | pd.isnull(ds['snow_depth'].sel(time=closest_ts)))\n",
"diff = ds['lidar-sd'].where(mask) - ds['snow_depth'].sel(time=closest_ts).where(mask)\n",
"rmse = float(np.sqrt((diff**2).sum()/len(diff.values.flatten())))\n",
"print(f'Full RMSE = {rmse:0.2f}')\n",
"#rmse_no_flag.loc[a,b,c,wst] = rmse\n",
"#valid_pixels.loc[a,b,c,wst] = mask_wet.sum() / mask.sum()\n",
"print(f'Frac valid pixels = {mask_wet.sum()/ mask.sum():0.2f}')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f,ax=plt.subplots(1,2,figsize=(10,4))\n",
"mask.plot(ax=ax[0])\n",
"mask_wet.plot(ax=ax[1])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"files = sorted(glob('spicy_s1_stacks/*.nc'))\n",
"\n",
"# Create parameter space\n",
"a = 2.5\n",
"b = 0.2\n",
"c = 0.55\n",
"\n",
"wet_snow_thresh = np.arange(-3, -0.9, 0.1)\n",
"freeze_snow_thresh = np.arange(1, 3.1, 0.1)\n",
"SI_thresh = [0,-100]\n",
"\n",
"total_count = len(wet_snow_thresh)*len(freeze_snow_thresh)*len(SI_thresh)\n",
"\n",
"for f in files:\n",
" ds_name = f.split('stacks/')[-1].split('.')[0]\n",
" print(datetime.now(), f' -- starting {ds_name}')\n",
"\n",
" if Path(f'rmse_test_wet_snow/{ds_name}_wet_flag.nc').is_file():\n",
" print('This file already exists, continuing.')\n",
" continue\n",
"\n",
" # Open dataset \n",
" ds_ = xr.open_dataset(f).load()\n",
" dataset = ds_[['s1','deltaVV','ims','fcf','lidar-sd']]\n",
" td = abs(pd.to_datetime(dataset.time) - pd.to_datetime(dataset.attrs['lidar-flight-time']))\n",
" closest_ts_idx = np.where(td == td.min())[0][0]\n",
" closest_ts = dataset.time[closest_ts_idx]\n",
"\n",
" if 'Frasier_2020-02-11' in ds_name:\n",
" closest_ts = '2020-02-16T13:09:43.000000000'\n",
" \n",
" # Initialize RMSE arrays\n",
" rmse_wet_flag = xr.DataArray(np.empty((len(SI_thresh),len(wet_snow_thresh), len(freeze_snow_thresh)))*np.nan,\n",
" coords=(SI_thresh,wet_snow_thresh,freeze_snow_thresh), dims=('SI_thresh','wet_snow_thresh','freeze_snow_thresh'))\n",
" rmse_no_flag = xr.DataArray(np.empty((len(SI_thresh),len(wet_snow_thresh), len(freeze_snow_thresh)))*np.nan,\n",
" coords=(SI_thresh,wet_snow_thresh,freeze_snow_thresh), dims=('SI_thresh','wet_snow_thresh','freeze_snow_thresh'))\n",
" valid_pixels = xr.DataArray(np.empty((len(SI_thresh),len(wet_snow_thresh), len(freeze_snow_thresh)))*np.nan,\n",
" coords=(SI_thresh,wet_snow_thresh,freeze_snow_thresh), dims=('SI_thresh','wet_snow_thresh','freeze_snow_thresh'))\n",
" # Brute-force loop\n",
" for wst in wet_snow_thresh:\n",
" for fst in freeze_snow_thresh:\n",
" for sit in SI_thresh:\n",
" print(f'sit={sit:0.2f}, wst={wst:0.2f}; fst={fst:0.2f}')\n",
"\n",
" ds = retrieval_from_parameters(dataset,A=a,B=b,C=c,wet_SI_thresh=sit,freezing_snow_thresh=fst,wet_snow_thres=wst)\n",
"\n",
" mask_wet = ~(ds['lidar-sd'].isnull() | ds['snow_depth'].sel(time=closest_ts).isnull() | ds['wet_snow'].sel(time=closest_ts).astype(bool))\n",
" diff_wet = ds['lidar-sd'].where(mask_wet) - ds['snow_depth'].sel(time=closest_ts).where(mask_wet)\n",
" rmse_wet = float(np.sqrt((diff_wet**2).sum()/len(diff_wet.values.flatten())))\n",
" print(f'RMSE with wet snow masked out = {rmse_wet:0.2f}')\n",
" rmse_wet_flag.loc[sit,wst,fst] = rmse_wet\n",
" # Compare snow depths - no wet snow mask\n",
" mask = ~(ds['lidar-sd'].isnull() | pd.isnull(ds['snow_depth'].sel(time=closest_ts)))\n",
" diff = ds['lidar-sd'].where(mask) - ds['snow_depth'].sel(time=closest_ts).where(mask)\n",
" rmse = float(np.sqrt((diff**2).sum()/len(diff.values.flatten())))\n",
" print(f'Full RMSE = {rmse:0.2f}')\n",
" rmse_no_flag.loc[sit,wst,fst] = rmse\n",
" valid_pixels.loc[sit,wst,fst] = mask_wet.sum() / mask.sum()\n",
" print(f'Frac valid pixels = {mask_wet.sum()/ mask.sum():0.2f}')\n",
"\n",
"\n",
" # After loop, save RMSE results per file\n",
" rmse_wet_flag.to_netcdf(f'rmse_test_wet_snow/{ds_name}_wet_flag.nc')\n",
" rmse_no_flag.to_netcdf(f'rmse_test_wet_snow/{ds_name}_no_flag.nc')\n",
" valid_pixels.to_netcdf(f'rmse_test_wet_snow/{ds_name}_valid_pixels.nc')\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"directory = 'rmse_test_wet_snow'\n",
"\n",
"\n",
"which_site = 0\n",
"\n",
"results1 = sorted(glob(f'{directory}/*wet*.nc'))\n",
"results2 = sorted(glob(f'{directory}/*no*.nc'))\n",
"results3 = sorted(glob(f'{directory}/*valid*.nc'))\n",
"\n",
"wet_snow = xr.open_dataarray(results1[which_site])\n",
"all_snow = xr.open_dataarray(results2[which_site])\n",
"frac_valid = xr.open_dataarray(results3[which_site])\n",
"\n",
"all_rmse = xr.concat([wet_snow,all_snow],'wet_or_all')\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f,ax=plt.subplots(1,2)\n",
"wet_snow.sel(SI_thresh=0).plot(ax=ax[0])\n",
"frac_valid.sel(SI_thresh=0).plot(ax=ax[1])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"all_rmse"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sit = 0\n",
"f=all_rmse.sel(SI_thresh=sit).plot(hue='wet_or_all',col='wet_snow_thresh',add_legend=False)\n",
"for wst,ax in zip(wet_snow_thresh,f.axs[0]):\n",
" frac_ax = ax.twinx()\n",
" fv = frac_valid.sel(SI_thresh=sit,wet_snow_thresh=wst).plot(ax=frac_ax,color='green',label='dry pixel fraction')\n",
" frac_ax.set_title('')\n",
" ax.axvline(wet_snow.sel(SI_thresh=sit,wet_snow_thresh=wst).idxmin(),color='black',linestyle='--')\n",
" ax.set_title('')\n",
" dry_percent = 100*frac_valid.sel(SI_thresh=sit,wet_snow_thresh=wst,freeze_snow_thresh=float(wet_snow.sel(SI_thresh=sit,wet_snow_thresh=wst).idxmin()))\n",
" ax.set_title(f'sit={sit:0.1f}, wst={wst:0.1f}, \\n min(RMSE)={float(wet_snow.sel(SI_thresh=sit,wet_snow_thresh=wst).min()):0.2f} @ {float(wet_snow.sel(SI_thresh=sit,wet_snow_thresh=wst).idxmin()):0.2f}dB,\\n Dry={dry_percent:0.2f}%')\n",
"\n",
" \n",
"ax.legend(labels=['wet snow mask','no mask'], title= 'RMSE', loc='lower right')\n",
"frac_ax.legend(handles=fv,labels=['Dry pixel fraction'], loc='upper right')\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.3"
},
"vscode": {
"interpreter": {
"hash": "ce4bdd2387e2daa803a7d0f8b0d766d25a1c9eab6b20981c1c0786f34d7ccd75"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
275 changes: 275 additions & 0 deletions contrib/gagliano/wet_snow_plot_rmse.ipynb
View file

Large diffs are not rendered by default.

Loading

0 comments on commit 8258cf9

Please sign in to comment.