Merge branch 'main' of github.com:SnowEx/spicy-snow into main

contrib/gagliano/wet_snow_get_rmse_all_sites.ipynb

@@ -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
+}