Skip to content

Commit

Permalink
Merge branch 'main' of github.com:NREL/gcm_eval
Browse files Browse the repository at this point in the history
  • Loading branch information
@grantbuster
grantbuster committed Mar 11, 2024
2 parents a0ab999 + 3d65c55 commit 43c0fdc
Show file tree
Hide file tree
Showing 2 changed files with 165 additions and 103 deletions.
120 changes: 52 additions & 68 deletions analysis/collect_skill_assess.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,10 @@
"cells": [
{
"cell_type": "markdown",
"id": "f2aa99e9-8272-4661-85ea-f062cb06ba63",
"id": "09196158-33a8-414f-80c9-6a2b9739b7b5",
"metadata": {},
"source": [
"# Collect Skill Assessments\n",
"This notebook collects the sup3r skill assessment outputs"
"# Collect sup3r Skill Outputs and Aggregate Per Region"
]
},
{
Expand Down Expand Up @@ -97,6 +96,51 @@
"fp_base = '/projects/alcaps/gcm_eval/{}_historical*/*.h5'"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d78b93a7-47f7-496f-9b2c-248134012d63",
"metadata": {},
"outputs": [],
"source": [
"def add_data_to_df(i, df, data, mask, base_feature, bias_feature, model, tag, rname):\n",
" \n",
" df.at[i, 'gcm'] = model\n",
" df.at[i, 'tag'] = tag\n",
" df.at[i, 'region'] = rname\n",
" df.at[i, 'feature'] = bias_feature\n",
"\n",
" df.at[i, 'hist_mean'] = data[f'base_{base_feature}_mean'][mask].mean()\n",
" df.at[i, 'hist_p1'] = data[f'base_{base_feature}_percentile_1'][mask].mean()\n",
" df.at[i, 'hist_p5'] = data[f'base_{base_feature}_percentile_5'][mask].mean()\n",
" df.at[i, 'hist_p50'] = data[f'base_{base_feature}_percentile_50'][mask].mean()\n",
" df.at[i, 'hist_p95'] = data[f'base_{base_feature}_percentile_95'][mask].mean()\n",
" df.at[i, 'hist_p99'] = data[f'base_{base_feature}_percentile_99'][mask].mean()\n",
" df.at[i, 'gcm_mean'] = data[f'bias_{bias_feature}_mean'][mask].mean()\n",
" df.at[i, 'gcm_p1'] = data[f'bias_{bias_feature}_percentile_1'][mask].mean()\n",
" df.at[i, 'gcm_p5'] = data[f'bias_{bias_feature}_percentile_5'][mask].mean()\n",
" df.at[i, 'gcm_p50'] = data[f'bias_{bias_feature}_percentile_50'][mask].mean()\n",
" df.at[i, 'gcm_p95'] = data[f'bias_{bias_feature}_percentile_95'][mask].mean()\n",
" df.at[i, 'gcm_p99'] = data[f'bias_{bias_feature}_percentile_99'][mask].mean()\n",
" \n",
" df.at[i, 'ks_stat'] = data[f'{bias_feature}_ks_stat'][mask].mean()\n",
" df.at[i, 'bias_mean'] = (data[f'bias_{bias_feature}_mean'][mask] - data[f'base_{base_feature}_mean'][mask]).mean()\n",
" df.at[i, 'bias_p1'] = (data[f'bias_{bias_feature}_percentile_1'][mask] - data[f'base_{base_feature}_percentile_1'][mask]).mean()\n",
" df.at[i, 'bias_p5'] = (data[f'bias_{bias_feature}_percentile_5'][mask] - data[f'base_{base_feature}_percentile_5'][mask]).mean()\n",
" df.at[i, 'bias_p50'] = (data[f'bias_{bias_feature}_percentile_50'][mask] - data[f'base_{base_feature}_percentile_50'][mask]).mean()\n",
" df.at[i, 'bias_p95'] = (data[f'bias_{bias_feature}_percentile_95'][mask] - data[f'base_{base_feature}_percentile_95'][mask]).mean()\n",
" df.at[i, 'bias_p99'] = (data[f'bias_{bias_feature}_percentile_99'][mask] - data[f'base_{base_feature}_percentile_99'][mask]).mean()\n",
"\n",
" df.at[i, 'percent_bias_mean'] = 100 * df.at[i, 'bias_mean'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p1'] = 100 * df.at[i, 'bias_p1'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p5'] = 100 * df.at[i, 'bias_p5'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p50'] = 100 * df.at[i, 'bias_p50'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p95'] = 100 * df.at[i, 'bias_p95'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p99'] = 100 * df.at[i, 'bias_p99'] / df.at[i, 'hist_mean']\n",
"\n",
" return df"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down Expand Up @@ -135,40 +179,12 @@
" \n",
" for rname, rstates in REGIONS.items():\n",
" mask = np.isin(meta[states_col].values.reshape(lat.shape), rstates)\n",
" \n",
" df.at[i, 'gcm'] = model\n",
" df.at[i, 'tag'] = tag\n",
" df.at[i, 'region'] = rname\n",
" df.at[i, 'feature'] = bias_feature\n",
"\n",
" df.at[i, 'hist_mean'] = data[f'base_{base_feature}_mean'][mask].mean()\n",
" df.at[i, 'hist_p1'] = data[f'base_{base_feature}_percentile_1'][mask].mean()\n",
" df.at[i, 'hist_p5'] = data[f'base_{base_feature}_percentile_5'][mask].mean()\n",
" df.at[i, 'hist_p50'] = data[f'base_{base_feature}_percentile_50'][mask].mean()\n",
" df.at[i, 'hist_p95'] = data[f'base_{base_feature}_percentile_95'][mask].mean()\n",
" df.at[i, 'hist_p99'] = data[f'base_{base_feature}_percentile_99'][mask].mean()\n",
" df.at[i, 'gcm_mean'] = data[f'bias_{bias_feature}_mean'][mask].mean()\n",
" df.at[i, 'gcm_p1'] = data[f'bias_{bias_feature}_percentile_1'][mask].mean()\n",
" df.at[i, 'gcm_p5'] = data[f'bias_{bias_feature}_percentile_5'][mask].mean()\n",
" df.at[i, 'gcm_p50'] = data[f'bias_{bias_feature}_percentile_50'][mask].mean()\n",
" df.at[i, 'gcm_p95'] = data[f'bias_{bias_feature}_percentile_95'][mask].mean()\n",
" df.at[i, 'gcm_p99'] = data[f'bias_{bias_feature}_percentile_99'][mask].mean()\n",
" \n",
" df.at[i, 'ks_stat'] = data[f'{bias_feature}_ks_stat'][mask].mean()\n",
" df.at[i, 'bias_mean'] = (data[f'bias_{bias_feature}_mean'][mask] - data[f'base_{base_feature}_mean'][mask]).mean()\n",
" df.at[i, 'bias_p1'] = (data[f'bias_{bias_feature}_percentile_1'][mask] - data[f'base_{base_feature}_percentile_1'][mask]).mean()\n",
" df.at[i, 'bias_p5'] = (data[f'bias_{bias_feature}_percentile_5'][mask] - data[f'base_{base_feature}_percentile_5'][mask]).mean()\n",
" df.at[i, 'bias_p50'] = (data[f'bias_{bias_feature}_percentile_50'][mask] - data[f'base_{base_feature}_percentile_50'][mask]).mean()\n",
" df.at[i, 'bias_p95'] = (data[f'bias_{bias_feature}_percentile_95'][mask] - data[f'base_{base_feature}_percentile_95'][mask]).mean()\n",
" df.at[i, 'bias_p99'] = (data[f'bias_{bias_feature}_percentile_99'][mask] - data[f'base_{base_feature}_percentile_99'][mask]).mean()\n",
"\n",
" df.at[i, 'percent_bias_mean'] = 100 * df.at[i, 'bias_mean'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p1'] = 100 * df.at[i, 'bias_p1'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p5'] = 100 * df.at[i, 'bias_p5'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p50'] = 100 * df.at[i, 'bias_p50'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p95'] = 100 * df.at[i, 'bias_p95'] / df.at[i, 'hist_mean']\n",
" df.at[i, 'percent_bias_p99'] = 100 * df.at[i, 'bias_p99'] / df.at[i, 'hist_mean']\n",
" df = add_data_to_df(i, df, data, mask, base_feature, bias_feature, model, tag, rname)\n",
" i += 1\n",
"\n",
" for rname in ['atlantic', 'pacific', 'gulf']:\n",
" mask = meta[rname].values.reshape(lat.shape)\n",
" df = add_data_to_df(i, df, data, mask, base_feature, bias_feature, model, tag, rname)\n",
" i += 1"
]
},
Expand Down Expand Up @@ -232,38 +248,6 @@
" dfp = df[mask].pivot(index='gcm', columns='feature', values=metric).T.loc[index]\n",
" dfp.to_csv(f'./skill_summaries/skill_summary_{rstr}_{metric}.csv')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ba69178e-f79c-45e5-adf4-e5fc4363d5d2",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "5a1773f8-33cb-4b42-beef-537c0227cb74",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "499cc21a-1953-45fb-af54-2e04bf636ace",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "4b42df7c-039d-4f33-972e-0ddd8d9c72dc",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down
148 changes: 113 additions & 35 deletions analysis/plot_trends.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,17 +2,17 @@
"cells": [
{
"cell_type": "markdown",
"id": "6e80956e-e918-41a7-9d43-878a43567e64",
"id": "811ce75d-6a00-4e8a-9613-8cc9b3ba86e9",
"metadata": {},
"source": [
"# Plot Trends\n",
"Use the projection summary files to make trend plots"
"# Make Trend Timeseries Plots\n",
"Plots will compare historical trends vs. GCM projections per region"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a95a6131-3d84-453e-b235-0eb01f4b6e64",
"id": "332bdc63-03a3-4a19-861b-0aa30684b141",
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -156,6 +156,21 @@
" 'd': 'diamond-tall',\n",
"}\n",
"\n",
"symbol_map_3d = {\n",
" 'cesm2': 'circle',\n",
" 'cesm2waccm': 'circle',\n",
" 'ecearth3cc': 'circle-open',\n",
" 'ecearth3': 'circle-open',\n",
" 'ecearth3veg': 'circle-open',\n",
" 'gfdlcm4': 'cross',\n",
" 'gfdlesm4': 'cross',\n",
" 'inmcm48': 'diamond',\n",
" 'inmcm50': 'diamond',\n",
" 'mpiesm12hr': 'diamond-open',\n",
" 'mriesm20': 'square',\n",
" 'noresm2mm': 'square-open',\n",
" 'taiesm1': 'x',}\n",
"\n",
"markers = {'cesm2': 'h',\n",
" 'cesm2waccm': 'H',\n",
" 'ecearth3cc': '^',\n",
Expand Down Expand Up @@ -266,12 +281,12 @@
" idf = idf.groupby(idf.index.date).max()\n",
" elif '_min_' in feature:\n",
" idf = idf.groupby(idf.index.date).min()\n",
" \n",
"\n",
" if df is None:\n",
" df = idf\n",
" else:\n",
" df = df.join(idf, how='outer')\n",
"\n",
" \n",
" # drop leap days (some GCMs have them, some do not)\n",
" mask = (df.index.month == 2) & (df.index.day == 29)\n",
" df = df[~mask]\n",
Expand All @@ -283,6 +298,12 @@
" if period is not None and option == 'min':\n",
" df = df.rolling(period, center=True).min()\n",
"\n",
" if tslice is not None:\n",
" df = df.iloc[tslice]\n",
"\n",
" if years is not None:\n",
" df = df[df.index.year.isin(years)]\n",
"\n",
" mask = df.index.year.isin(range(1980,2020))\n",
" if baseline is True:\n",
" for col in df.columns:\n",
Expand All @@ -309,12 +330,6 @@
" df[col] = df[col] - value0 + baseline\n",
"\n",
" df = df.dropna(axis=0, how='all')\n",
"\n",
" if tslice is not None:\n",
" df = df.iloc[tslice]\n",
"\n",
" if years is not None:\n",
" df = df[df.index.year.isin(years)]\n",
" \n",
" return df"
]
Expand Down Expand Up @@ -435,6 +450,74 @@
" fig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "92f648ee-a2a2-42ff-ad98-04657a028b80",
"metadata": {},
"outputs": [],
"source": [
"def scatter_3d(all_df, title, fp_out, show):\n",
" df = pd.DataFrame(columns=list(all_df.keys()))\n",
" symbols_3d = []\n",
" colors_3d = []\n",
" \n",
" for i, (model, tag) in enumerate(zip(MODELS, TAGS)):\n",
" if not np.isnan(all_df['temperature_2m'][model].values[-1]):\n",
" symbols_3d.append(symbol_map_3d.get(tag, 'circle'))\n",
" colors_3d.append(colors.get(tag, 'k'))\n",
" df.at[i, 'Model'] = model\n",
" for key, feat_df in all_df.items():\n",
" value = feat_df[model].values[-1] - feat_df[model].values[:7].mean()\n",
" if np.isnan(value):\n",
" value = 0.0\n",
" df.at[i, key] = np.round(value, 3)\n",
"\n",
" df['temperature_2m_degf'] = df['temperature_2m'] * 9/5\n",
" df['temperature_max_2m_degf'] = df['temperature_max_2m'] * 9/5\n",
" df['temperature_min_2m_degf'] = df['temperature_min_2m'] * 9/5\n",
" df['temperature_2m_degf'] = [np.round(x, 3) for x in df['temperature_2m_degf']]\n",
" df['temperature_max_2m_degf'] = [np.round(x, 3) for x in df['temperature_max_2m_degf']]\n",
" df['temperature_min_2m_degf'] = [np.round(x, 3) for x in df['temperature_min_2m_degf']]\n",
" \n",
" offset=1\n",
" fig = px.scatter_3d(df, \n",
" x='windspeed_100m', \n",
" y='pr', \n",
" z='temperature_2m', \n",
" color='Model', \n",
" color_discrete_sequence=colors_3d, \n",
" symbol='Model', \n",
" symbol_sequence=symbols_3d,\n",
" hover_data=list(df.columns),\n",
" range_x=(df['windspeed_100m'].min()-1, df['windspeed_100m'].max()+1),\n",
" range_y=(df['pr'].min()-1, df['pr'].max()+1),\n",
" range_z=(df['temperature_2m'].min()-1, df['temperature_2m'].max()+1),\n",
" width=600, height=400,\n",
" title=title,\n",
" )\n",
" \n",
" camera = dict(\n",
" up=dict(x=0, y=0, z=1),\n",
" center=dict(x=.15, y=0, z=-0.3),\n",
" eye=dict(x=1.25, y=1.25, z=1.25)\n",
" )\n",
" \n",
" fig.update_layout(scene_aspectmode='cube',\n",
" scene_camera=camera,\n",
" scene=dict(\n",
" xaxis_title='WS Change (%)',\n",
" yaxis_title='Precip Change (%)',\n",
" zaxis_title='Temp Change (°C)'),\n",
" margin=dict(r=0, b=0, l=0, t=50))\n",
"\n",
" if fp_out is not None:\n",
" fp_out = fp_out.replace('.png', '.html')\n",
" fig.write_html(fp_out)\n",
" if show:\n",
" fig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -443,7 +526,6 @@
"outputs": [],
"source": [
"all_regions = list(REGIONS.keys())\n",
"# all_regions = []\n",
"all_regions += ['atlantic', 'pacific', 'gulf']"
]
},
Expand All @@ -456,6 +538,8 @@
"source": [
"for region in all_regions:\n",
"# for region in ['conus']:\n",
" all_df_245 = {}\n",
" all_df_585 = {}\n",
" display(region)\n",
" region = region.lower().replace(' ', '_')\n",
" \n",
Expand All @@ -466,7 +550,7 @@
" 'relativehumidity_2m': dict(period=365*10, baseline=True, tslice=slice(None, None, 365*5), relative=False),\n",
" 'rsds': dict(period=365*10, baseline=True, tslice=slice(None, None, 365*5), relative=True, years=list(range(2005, 2055))),\n",
" 'pr': dict(period=365*10, baseline=True, tslice=slice(None, None, 365*5), relative=True),\n",
" 'windspeed_100m': dict(period=365*10, baseline=True, tslice=slice(None, None, 365*5), relative=True, years=list(range(2005, 2055))),\n",
" 'windspeed_100m': dict(period=365*10, baseline=True, tslice=slice(None, None, 365*5), relative=True),\n",
" }\n",
" \n",
" plot_kwargs = {\n",
Expand Down Expand Up @@ -494,38 +578,32 @@
" }\n",
"\n",
" for feature in FEATURES:\n",
" # for feature in ['rsds']:\n",
" df_245 = get_trend_df(region, 'ssp245', feature, **data_kwargs[feature])\n",
" df_585 = get_trend_df(region, 'ssp585', feature, **data_kwargs[feature])\n",
" all_df_245[feature] = df_245\n",
" all_df_585[feature] = df_585\n",
" \n",
" plot(df_245, df_585, show=False, **plot_kwargs[feature])\n",
" plotly(df_245, df_585, show=False, **plot_kwargs[feature])\n",
"\n",
" if feature.startswith('temperature'):\n",
" df_245_degf = df_245.copy()\n",
" df_585_degf = df_585.copy()\n",
" plot_kwargs[feature]['ylabel'] = plot_kwargs[feature]['ylabel'].replace('$^\\circ$C', '$^\\circ$F')\n",
" plot_kwargs[feature]['fp_out'] = plot_kwargs[feature]['fp_out'].replace('.png', '_degf.png')\n",
" df_245 *= 9/5\n",
" df_585 *= 9/5\n",
" df_245_degf *= 9/5\n",
" df_585_degf *= 9/5\n",
" if '_max_' in feature or '_min_' in feature:\n",
" df_245 += 32\n",
" df_585 += 32\n",
" plotly(df_245, df_585, show=False, **plot_kwargs[feature])"
" df_245_degf += 32\n",
" df_585_degf += 32\n",
" plotly(df_245_degf, df_585_degf, show=False, **plot_kwargs[feature])\n",
" \n",
" scatter_3d(all_df_245, 'SSP2 4.5', f'./trend_plots/{region}_scatter_ssp245.png', show=False)\n",
" scatter_3d(all_df_585, 'SSP5 8.5', f'./trend_plots/{region}_scatter_ssp585.png', show=False)\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ae035fac-b3be-46a8-8721-84b371a1704d",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "23b14bc1-bccb-44a2-91d7-f1b6764c100c",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down

0 comments on commit 43c0fdc

Please sign in to comment.