diff --git a/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb b/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb index a7a9f61..f7ca899 100644 --- a/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb +++ b/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb @@ -65,10 +65,11 @@ "\n", "climo_nyears = 35\n", "grid_file = \"/glade/campaign/cesm/community/omwg/grids/tx2_3v2_grid.nc\"\n", + "path_nsidc = (\n", + " \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_obs_data/ice/\"\n", + ")\n", "\n", - "lc_kwargs = {}\n", - "\n", - "hist = 0" + "lc_kwargs = {}" ] }, { @@ -141,7 +142,7 @@ " CESM_output_dir\n", " + \"/\"\n", " + case_name\n", - " + \"/ice/proc/tseries/month_1/\"\n", + " + \"/ice/proc/tseries/\"\n", " + case_name\n", " + \".cice.h.\"\n", " + \"*.nc\",\n", @@ -153,7 +154,7 @@ " CESM_output_dir\n", " + \"/\"\n", " + base_case_name\n", - " + \"/ice/proc/tseries/month_1/\"\n", + " + \"/ice/proc/tseries/\"\n", " + base_case_name\n", " + \".cice.h.\"\n", " + \"*.nc\",\n", diff --git a/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb b/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb index 33dedb2..7a9e0e2 100644 --- a/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb +++ b/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb @@ -38,7 +38,7 @@ "import yaml\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", - "import nc_time_axis" + "import matplotlib.lines as mlines" ] }, { @@ -72,9 +72,11 @@ "grid_file = \"/glade/campaign/cesm/community/omwg/grids/tx2_3v2_grid.nc\"\n", "path_nsidc = \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_obs_data/ice/analysis_datasets/hemispheric_data/NSIDC_SeaIce_extent\"\n", "\n", - "lc_kwargs = {}\n", + "model_path = \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_model_data/ice/\"\n", + "path_lens1 = model_path + \"cesm_lens1\"\n", + "path_lens2 = model_path + \"cesm_lens2\"\n", "\n", - "hist = 0" + "lc_kwargs = {}" ] }, { @@ -176,13 +178,18 @@ "ds2_ann = ds2.resample(time=\"YS\").mean(dim=\"time\")\n", "\n", "with open(\"cice_masks.yml\", \"r\") as file:\n", - " cice_masks = yaml.safe_load(file)" + " cice_masks = yaml.safe_load(file)\n", + "\n", + "first_year = int(start_date.split(\"-\")[0])\n", + "base_first_year = int(base_start_date.split(\"-\")[0])\n", + "end_year = int(end_date.split(\"-\")[0])\n", + "base_end_year = int(base_end_date.split(\"-\")[0])" ] }, { "cell_type": "code", "execution_count": null, - "id": "173299ea-d64e-48ad-aa9a-40663dd5c7e9", + "id": "31c8d5bf-b76f-46b2-bb5b-a654fa2150ed", "metadata": { "editable": true, "slideshow": { @@ -194,52 +201,106 @@ }, "outputs": [], "source": [ - "### Read in the CESM1 and CESM2 LE sea ice area\n", + "### Read in the CESM LENS historical data\n", + "\n", + "ds_cesm1_aicetot_nh = xr.open_dataset(path_lens1 + \"/LE_aicetot_nh_1920-2100.nc\")\n", + "ds_cesm1_hitot_nh = xr.open_dataset(path_lens1 + \"/LE_hitot_nh_1920-2100.nc\")\n", + "ds_cesm1_hstot_nh = xr.open_dataset(path_lens1 + \"/LE_hstot_nh_1920-2100.nc\")\n", + "\n", + "ds_cesm1_aicetot_sh = xr.open_dataset(path_lens1 + \"/LE_aicetot_sh_1920-2100.nc\")\n", + "ds_cesm1_hitot_sh = xr.open_dataset(path_lens1 + \"/LE_hitot_sh_1920-2100.nc\")\n", + "ds_cesm1_hstot_sh = xr.open_dataset(path_lens1 + \"/LE_hstot_sh_1920-2100.nc\")\n", + "\n", + "cesm1_aicetot_nh_ann = ds_cesm1_aicetot_nh[\"aice_monthly\"].mean(dim=\"nmonth\")\n", + "cesm1_hitot_nh_ann = ds_cesm1_hitot_nh[\"hi_monthly\"].mean(dim=\"nmonth\")\n", + "cesm1_hstot_nh_ann = ds_cesm1_hstot_nh[\"hs_monthly\"].mean(dim=\"nmonth\")\n", + "\n", + "cesm1_aicetot_sh_ann = ds_cesm1_aicetot_sh[\"aice_monthly\"].mean(dim=\"nmonth\")\n", + "cesm1_hitot_sh_ann = ds_cesm1_hitot_sh[\"hi_monthly\"].mean(dim=\"nmonth\")\n", + "cesm1_hstot_sh_ann = ds_cesm1_hstot_sh[\"hs_monthly\"].mean(dim=\"nmonth\")\n", + "\n", + "cesm1_years = xr.cftime_range(start=\"1920\", end=\"2100\", freq=\"YS\", calendar=\"noleap\")\n", + "\n", + "cesm1_aicetot_nh_month = ds_cesm1_aicetot_nh[\"aice_monthly\"][:, 60:95, :].mean(\n", + " dim=\"nyr\"\n", + ")\n", + "cesm1_hitot_nh_month = ds_cesm1_hitot_nh[\"hi_monthly\"][:, 60:95, :].mean(dim=\"nyr\")\n", + "cesm1_hstot_nh_month = ds_cesm1_hstot_nh[\"hs_monthly\"][:, 60:95, :].mean(dim=\"nyr\")\n", + "cesm1_aicetot_sh_month = ds_cesm1_aicetot_sh[\"aice_monthly\"][:, 60:95, :].mean(\n", + " dim=\"nyr\"\n", + ")\n", + "cesm1_hitot_sh_month = ds_cesm1_hitot_sh[\"hi_monthly\"][:, 60:95, :].mean(dim=\"nyr\")\n", + "cesm1_hstot_sh_month = ds_cesm1_hstot_sh[\"hs_monthly\"][:, 60:95, :].mean(dim=\"nyr\")\n", "\n", - "# cesm1_path = \"/glade/campaign/cesm/collections/cesmLE/CESM-CAM5-BGC-LE/ice/proc/tseries/monthly/aice/\"\n", + "ds_cesm2_aicetot_nh = xr.open_dataset(path_lens2 + \"/LE2_aicetot_nh_1870-2100.nc\")\n", + "ds_cesm2_hitot_nh = xr.open_dataset(path_lens2 + \"/LE2_hitot_nh_1870-2100.nc\")\n", + "ds_cesm2_hstot_nh = xr.open_dataset(path_lens2 + \"/LE2_hstot_nh_1870-2100.nc\")\n", "\n", - "# ds_cesm1_nh = xr.open_dataset(cesm1_path + \"b.e11.B1850C5CN.f09_g16.005.cice.h.aice_nh.210001-220012.nc\")\n", - "# ds_cesm1_sh = xr.open_dataset(cesm1_path + \"b.e11.B1850C5CN.f09_g16.005.cice.h.aice_sh.210001-220012.nc\")\n", + "ds_cesm2_aicetot_sh = xr.open_dataset(path_lens2 + \"/LE2_aicetot_sh_1870-2100.nc\")\n", + "ds_cesm2_hitot_sh = xr.open_dataset(path_lens2 + \"/LE2_hitot_sh_1870-2100.nc\")\n", + "ds_cesm2_hstot_sh = xr.open_dataset(path_lens2 + \"/LE2_hstot_sh_1870-2100.nc\")\n", "\n", - "# cesm2_path = \"/glade/campaign/collections/cmip/CMIP6/timeseries-cmip6/b.e21.B1850.f09_g17.CMIP6-piControl.001/ice/proc/tseries/month_1/\"\n", - "# cesm2_ds = xr.open_dataset(cesm2_path + \"b.e21.B1850.f09_g17.CMIP6-piControl.001.cice.h.aice.190001-200012.nc\")\n", + "cesm2_aicetot_nh_ann = ds_cesm2_aicetot_nh[\"aice_monthly\"].mean(dim=\"nmonth\")\n", + "cesm2_hitot_nh_ann = ds_cesm2_hitot_nh[\"hi_monthly\"].mean(dim=\"nmonth\")\n", + "cesm2_hstot_nh_ann = ds_cesm2_hstot_nh[\"hs_monthly\"].mean(dim=\"nmonth\")\n", "\n", - "# nj_nh = ds_cesm1_nh['nj']\n", - "# nj_sh = ds_cesm2_sh['nh']\n", + "cesm2_aicetot_sh_ann = ds_cesm2_aicetot_sh[\"aice_monthly\"].mean(dim=\"nmonth\")\n", + "cesm2_hitot_sh_ann = ds_cesm2_hitot_sh[\"hi_monthly\"].mean(dim=\"nmonth\")\n", + "cesm2_hstot_sh_ann = ds_cesm2_hstot_sh[\"hs_monthly\"].mean(dim=\"nmonth\")\n", "\n", - "# aice_cesm1 = cesm2_ds['aice']\n", + "cesm2_years = xr.cftime_range(start=\"1870\", end=\"2100\", freq=\"YS\", calendar=\"noleap\")\n", "\n", - "# aice_cesm1 = 0.0\n", - "# aice_cesm1[:,0:nj_sh-1,:] = ds_cesm1_sh['aice']\n", - "# aice_cesm1[:,383-nj_nh:383,:] = ds_cesm1_nh['aice']" + "cesm2_aicetot_nh_month = ds_cesm2_aicetot_nh[\"aice_monthly\"][:, 110:145, :].mean(\n", + " dim=\"nyr\"\n", + ")\n", + "cesm2_hitot_nh_month = ds_cesm2_hitot_nh[\"hi_monthly\"][:, 110:145, :].mean(dim=\"nyr\")\n", + "cesm2_hstot_nh_month = ds_cesm2_hstot_nh[\"hs_monthly\"][:, 110:145, :].mean(dim=\"nyr\")\n", + "cesm2_aicetot_sh_month = ds_cesm2_aicetot_sh[\"aice_monthly\"][:, 110:145, :].mean(\n", + " dim=\"nyr\"\n", + ")\n", + "cesm2_hitot_sh_month = ds_cesm2_hitot_sh[\"hi_monthly\"][:, 110:145, :].mean(dim=\"nyr\")\n", + "cesm2_hstot_sh_month = ds_cesm2_hstot_sh[\"hs_monthly\"][:, 110:145, :].mean(dim=\"nyr\")" ] }, { "cell_type": "code", "execution_count": null, "id": "4f8d57eb-8411-4b70-a200-35b99d3152f2", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, "outputs": [], "source": [ "# Two functions to help draw plots even if only one year of data is available\n", "# If one year is present, horizontal line will be dotted instead of solid\n", "\n", "\n", - "def da_plot_len_time_might_be_one(da_in, alt_time):\n", + "def da_plot_len_time_might_be_one(da_in, alt_time, color):\n", " # If da_in.time only has 1 value, draw horizontal line across range of alt_time\n", " if len(da_in.time) > 1:\n", - " da_in.plot()\n", + " da_in.plot(color=color)\n", " else:\n", " time_arr = [alt_time.data[0], alt_time.data[-1]]\n", - " plt.plot(time_arr, [da_in.data[0], da_in.data[0]], linestyle=\":\")\n", + " plt.plot(time_arr, [da_in.data[0], da_in.data[0]], linestyle=\":\", color=color)\n", "\n", "\n", - "def plt_plot_len_x_might_be_one(da_in, x_in, alt_x):\n", + "def plt_plot_len_x_might_be_one(da_in, x_in, alt_x, color):\n", " # If x_in only has one value, draw horizontal line across range of alt_x\n", " if len(x_in) > 1:\n", - " plt.plot(x_in, da_in)\n", + " plt.plot(x_in, da_in, color=color)\n", " else:\n", - " plt.plot([alt_x[0], alt_x[-1]], [da_in.data[0], da_in.data[0]], linestyle=\":\")" + " plt.plot(\n", + " [alt_x[0], alt_x[-1]],\n", + " [da_in.data[0], da_in.data[0]],\n", + " linestyle=\":\",\n", + " color=color,\n", + " )" ] }, { @@ -266,35 +327,70 @@ "ds1_vhs_ann = (tarea * ds1_ann[\"hs\"]).where(TLAT > 0).sum(dim=[\"nj\", \"ni\"]) * 1.0e-13\n", "ds2_vhs_ann = (tarea * ds2_ann[\"hs\"]).where(TLAT > 0).sum(dim=[\"nj\", \"ni\"]) * 1.0e-13\n", "\n", + "p1 = mlines.Line2D([], [], color=\"lightgrey\", label=\"CESM1-LENS\")\n", + "p2 = mlines.Line2D([], [], color=\"lightblue\", label=\"CESM2-LENS\")\n", + "p3 = mlines.Line2D([], [], color=\"red\", label=case_name)\n", + "p4 = mlines.Line2D([], [], color=\"blue\", label=base_case_name)\n", + "\n", "fig = plt.figure(figsize=(10, 10), tight_layout=True)\n", "\n", "ax = fig.add_subplot(3, 1, 1)\n", - "da_plot_len_time_might_be_one(ds1_vhi_ann, alt_time=ds2_ann.time)\n", - "da_plot_len_time_might_be_one(ds2_vhi_ann, alt_time=ds1_ann.time)\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " cesm1_years[0:95],\n", + " cesm1_hitot_nh_ann[i, 0:95] * 1.0e-13,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " cesm2_years[0:145], cesm2_hitot_nh_ann[i, 0:145] * 1.0e-13, color=\"lightblue\"\n", + " )\n", + "da_plot_len_time_might_be_one(ds1_vhi_ann, alt_time=ds2_ann.time, color=\"red\")\n", + "da_plot_len_time_might_be_one(ds2_vhi_ann, alt_time=ds1_ann.time, color=\"blue\")\n", "\n", "plt.title(\"NH Annual Mean Integrated Timeseries\")\n", "plt.ylim((0, 5))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"NH Annual Mean Sea Ice Volume $m x 10^{13}$\")\n", - "plt.legend([case_name, base_case_name])\n", + "plt.legend(handles=[p1, p2, p3, p4])\n", "\n", "ax = fig.add_subplot(3, 1, 2)\n", - "da_plot_len_time_might_be_one(ds1_vhs_ann, alt_time=ds2_ann.time)\n", - "da_plot_len_time_might_be_one(ds2_vhs_ann, alt_time=ds1_ann.time)\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " cesm1_years[0:95],\n", + " cesm1_hstot_nh_ann[i, 0:95] * 1.0e-13,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " cesm2_years[0:145], cesm2_hstot_nh_ann[i, 0:145] * 1.0e-13, color=\"lightblue\"\n", + " )\n", + "da_plot_len_time_might_be_one(ds1_vhs_ann, alt_time=ds2_ann.time, color=\"red\")\n", + "da_plot_len_time_might_be_one(ds2_vhs_ann, alt_time=ds1_ann.time, color=\"blue\")\n", "\n", "plt.ylim((0, 0.5))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"NH Annual Mean Snow Volume $m x 10^{13}$\")\n", - "plt.legend([case_name, base_case_name])\n", + "plt.legend(handles=[p1, p2, p3, p4])\n", "\n", "ax = fig.add_subplot(3, 1, 3)\n", - "da_plot_len_time_might_be_one(ds1_area_ann, alt_time=ds2_ann.time)\n", - "da_plot_len_time_might_be_one(ds2_area_ann, alt_time=ds1_ann.time)\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " cesm1_years[0:95],\n", + " cesm1_aicetot_nh_ann[i, 0:95] * 1.0e-12,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " cesm2_years[0:145], cesm2_aicetot_nh_ann[i, 0:145] * 1.0e-12, color=\"lightblue\"\n", + " )\n", + "da_plot_len_time_might_be_one(ds1_area_ann, alt_time=ds2_ann.time, color=\"red\")\n", + "da_plot_len_time_might_be_one(ds2_area_ann, alt_time=ds1_ann.time, color=\"blue\")\n", "\n", - "plt.ylim((0, 25))\n", + "plt.ylim((0, 15))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"NH Annual Mean Sea Ice Area $m x 10^{12}$\")\n", - "plt.legend([case_name, base_case_name])" + "plt.legend(handles=[p1, p2, p3, p4])" ] }, { @@ -324,32 +420,62 @@ "fig = plt.figure(figsize=(10, 10), tight_layout=True)\n", "\n", "ax = fig.add_subplot(3, 1, 1)\n", - "da_plot_len_time_might_be_one(ds1_vhi_ann, alt_time=ds2_ann.time)\n", - "da_plot_len_time_might_be_one(ds2_vhi_ann, alt_time=ds1_ann.time)\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " cesm1_years[0:95],\n", + " cesm1_hitot_sh_ann[i, 0:95] * 1.0e-13,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " cesm2_years[0:145], cesm2_hitot_sh_ann[i, 0:145] * 1.0e-13, color=\"lightblue\"\n", + " )\n", + "da_plot_len_time_might_be_one(ds1_vhi_ann, alt_time=ds2_ann.time, color=\"red\")\n", + "da_plot_len_time_might_be_one(ds2_vhi_ann, alt_time=ds1_ann.time, color=\"blue\")\n", "\n", "plt.title(\"SH Annual Mean Integrated Timeseries\")\n", - "plt.ylim((0, 5))\n", + "plt.ylim((0, 3))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"SH Annual Mean Sea Ice Volume $m x 10^{13}$\")\n", - "plt.legend([case_name, base_case_name])\n", + "plt.legend(handles=[p1, p2, p3, p4])\n", "\n", "ax = fig.add_subplot(3, 1, 2)\n", - "da_plot_len_time_might_be_one(ds1_vhs_ann, alt_time=ds2_ann.time)\n", - "da_plot_len_time_might_be_one(ds2_vhs_ann, alt_time=ds1_ann.time)\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " cesm1_years[0:95],\n", + " cesm1_hstot_sh_ann[i, 0:95] * 1.0e-13,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " cesm2_years[0:145], cesm2_hstot_sh_ann[i, 0:145] * 1.0e-13, color=\"lightblue\"\n", + " )\n", + "da_plot_len_time_might_be_one(ds1_vhs_ann, alt_time=ds2_ann.time, color=\"red\")\n", + "da_plot_len_time_might_be_one(ds2_vhs_ann, alt_time=ds1_ann.time, color=\"blue\")\n", "\n", "plt.ylim((0, 0.5))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"SH Annual Mean Snow Volume $m x 10^{13}$\")\n", - "plt.legend([case_name, base_case_name])\n", + "plt.legend(handles=[p1, p2, p3, p4])\n", "\n", "ax = fig.add_subplot(3, 1, 3)\n", - "da_plot_len_time_might_be_one(ds1_area_ann, alt_time=ds2_ann.time)\n", - "da_plot_len_time_might_be_one(ds2_area_ann, alt_time=ds1_ann.time)\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " cesm1_years[0:95],\n", + " cesm1_aicetot_sh_ann[i, 0:95] * 1.0e-12,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " cesm2_years[0:145], cesm2_aicetot_sh_ann[i, 0:145] * 1.0e-12, color=\"lightblue\"\n", + " )\n", + "da_plot_len_time_might_be_one(ds1_area_ann, alt_time=ds2_ann.time, color=\"red\")\n", + "da_plot_len_time_might_be_one(ds2_area_ann, alt_time=ds1_ann.time, color=\"blue\")\n", "\n", - "plt.ylim((0, 25))\n", + "plt.ylim((0, 20))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"SH Annual Mean Sea Ice Area $m x 10^{12}$\")\n", - "plt.legend([case_name, base_case_name])" + "plt.legend(handles=[p1, p2, p3, p4])" ] }, { @@ -432,7 +558,7 @@ "nov_ext = nov_nsidc.iloc[:, 4].values\n", "dec_ext = dec_nsidc.iloc[:, 4].values\n", "\n", - "nsidc_clim = [\n", + "nsidc_clim_ext = [\n", " np.nanmean(jan_ext[0:35]),\n", " np.nanmean(feb_ext[0:35]),\n", " np.nanmean(mar_ext[0:35]),\n", @@ -447,6 +573,21 @@ " np.nanmean(dec_ext[0:35]),\n", "]\n", "\n", + "nsidc_clim_area = [\n", + " np.nanmean(jan_area[0:35]),\n", + " np.nanmean(feb_area[0:35]),\n", + " np.nanmean(mar_area[0:35]),\n", + " np.nanmean(apr_area[0:35]),\n", + " np.nanmean(may_area[0:35]),\n", + " np.nanmean(jun_area[0:35]),\n", + " np.nanmean(jul_area[0:35]),\n", + " np.nanmean(aug_area[0:35]),\n", + " np.nanmean(sep_area[0:35]),\n", + " np.nanmean(oct_area[0:35]),\n", + " np.nanmean(nov_area[0:35]),\n", + " np.nanmean(dec_area[0:35]),\n", + "]\n", + "\n", "# plt.plot(nsidc_clim)" ] }, @@ -530,7 +671,7 @@ "nov_ext_sh = nov_nsidc_sh.iloc[:, 4].values\n", "dec_ext_sh = dec_nsidc_sh.iloc[:, 4].values\n", "\n", - "nsidc_clim_sh = [\n", + "nsidc_clim_sh_ext = [\n", " np.nanmean(jan_ext_sh[0:35]),\n", " np.nanmean(feb_ext_sh[0:35]),\n", " np.nanmean(mar_ext_sh[0:35]),\n", @@ -545,6 +686,20 @@ " np.nanmean(dec_ext_sh[0:35]),\n", "]\n", "\n", + "nsidc_clim_sh_area = [\n", + " np.nanmean(jan_area_sh[0:35]),\n", + " np.nanmean(feb_area_sh[0:35]),\n", + " np.nanmean(mar_area_sh[0:35]),\n", + " np.nanmean(apr_area_sh[0:35]),\n", + " np.nanmean(may_area_sh[0:35]),\n", + " np.nanmean(jun_area_sh[0:35]),\n", + " np.nanmean(jul_area_sh[0:35]),\n", + " np.nanmean(aug_area_sh[0:35]),\n", + " np.nanmean(sep_area_sh[0:35]),\n", + " np.nanmean(oct_area_sh[0:35]),\n", + " np.nanmean(nov_area_sh[0:35]),\n", + " np.nanmean(dec_area_sh[0:35]),\n", + "]\n", "# plt.plot(nsidc_clim_sh)" ] }, @@ -581,19 +736,36 @@ "mask_ext1 = xr.DataArray(data=mask_tmp1, dims=[\"month\", \"nj\", \"ni\"])\n", "mask_ext2 = xr.DataArray(data=mask_tmp2, dims=[\"month\", \"nj\", \"ni\"])\n", "\n", + "mask_nh_tmp = np.where(ds1[\"TLAT\"] > 0, tarea, 0.0)\n", + "mask_nh = xr.DataArray(data=mask_nh_tmp, dims=[\"nj\", \"ni\"])\n", "\n", "ext1 = (mask_ext1 * tarea).sum([\"ni\", \"nj\"]) * 1.0e-12\n", "ext2 = (mask_ext2 * tarea).sum([\"ni\", \"nj\"]) * 1.0e-12\n", "\n", - "plt.plot(ext1)\n", - "plt.plot(ext2)\n", - "plt.plot(nsidc_clim)\n", + "area1 = (aice1_month * mask_nh).sum([\"ni\", \"nj\"]) * 1.0e-12\n", + "area2 = (aice2_month * mask_nh).sum([\"ni\", \"nj\"]) * 1.0e-12\n", + "\n", + "p1 = mlines.Line2D([], [], color=\"lightgrey\", label=\"CESM1-LENS\")\n", + "p2 = mlines.Line2D([], [], color=\"lightblue\", label=\"CESM2-LENS\")\n", + "p3 = mlines.Line2D([], [], color=\"red\", label=case_name)\n", + "p4 = mlines.Line2D([], [], color=\"blue\", label=base_case_name)\n", + "p5 = mlines.Line2D([], [], color=\"black\", label=\"NSIDC\")\n", + "\n", + "months = np.linspace(1, 12, 12)\n", + "\n", + "for i in range(0, 38):\n", + " plt.plot(months, cesm1_aicetot_nh_month[i, :] * 1.0e-12, color=\"lightgrey\")\n", + "for i in range(0, 49):\n", + " plt.plot(months, cesm2_aicetot_nh_month[i, :] * 1.0e-12, color=\"lightblue\")\n", + "plt.plot(months, area1, color=\"red\")\n", + "plt.plot(months, area2, color=\"blue\")\n", + "plt.plot(months, nsidc_clim_area, color=\"black\")\n", "\n", "plt.title(\"NH Climatological Seasonal Cycle\")\n", "plt.ylim((0, 25))\n", "plt.xlabel(\"Month\")\n", - "plt.ylabel(\"Climatological Seasonal Cycle Ice Extent $m x 10^{12}$\")\n", - "plt.legend([case_name, base_case_name, \"NSIDC\"])" + "plt.ylabel(\"Climatological Seasonal Cycle Ice Area $m x 10^{12}$\")\n", + "plt.legend(handles=[p1, p2, p3, p4, p5])" ] }, { @@ -617,19 +789,28 @@ "mask_ext1_sh = xr.DataArray(data=mask_tmp1_sh, dims=[\"month\", \"nj\", \"ni\"])\n", "mask_ext2_sh = xr.DataArray(data=mask_tmp2_sh, dims=[\"month\", \"nj\", \"ni\"])\n", "\n", + "mask_sh_tmp = np.where(ds1[\"TLAT\"] < 0, tarea, 0.0)\n", + "mask_sh = xr.DataArray(data=mask_sh_tmp, dims=[\"nj\", \"ni\"])\n", "\n", "ext1_sh = (mask_ext1_sh * tarea).sum([\"ni\", \"nj\"]) * 1.0e-12\n", "ext2_sh = (mask_ext2_sh * tarea).sum([\"ni\", \"nj\"]) * 1.0e-12\n", "\n", - "plt.plot(ext1_sh)\n", - "plt.plot(ext2_sh)\n", - "plt.plot(nsidc_clim_sh)\n", + "area1 = (aice1_month * mask_sh).sum([\"ni\", \"nj\"]) * 1.0e-12\n", + "area2 = (aice2_month * mask_sh).sum([\"ni\", \"nj\"]) * 1.0e-12\n", + "\n", + "for i in range(0, 38):\n", + " plt.plot(months, cesm1_aicetot_sh_month[i, :] * 1.0e-12, color=\"lightgrey\")\n", + "for i in range(0, 49):\n", + " plt.plot(months, cesm2_aicetot_sh_month[i, :] * 1.0e-12, color=\"lightblue\")\n", + "plt.plot(months, area1, color=\"red\")\n", + "plt.plot(months, area2, color=\"blue\")\n", + "plt.plot(months, nsidc_clim_sh_area, color=\"black\")\n", "\n", "plt.title(\"SH Climatological Seasonal Cycle\")\n", - "plt.ylim((0, 25))\n", + "plt.ylim((0, 30))\n", "plt.xlabel(\"Month\")\n", - "plt.ylabel(\"Climatological Seasonal Cycle Ice Extent $m x 10^{12}$\")\n", - "plt.legend([case_name, base_case_name, \"NSIDC\"])" + "plt.ylabel(\"Climatological Seasonal Cycle Ice Area $m x 10^{12}$\")\n", + "plt.legend(handles=[p1, p2, p3, p4, p5])" ] }, { @@ -657,19 +838,39 @@ "ds1_sep = ds1_area.sel(time=(ds1_area.time.dt.month == 9))\n", "ds2_sep = ds2_area.sel(time=(ds2_area.time.dt.month == 9))\n", "\n", - "x1 = np.linspace(1, len(ds1_sep.time), len(ds1_sep.time))\n", - "x2 = np.linspace(1, len(ds2_sep.time), len(ds2_sep.time))\n", - "x3 = np.linspace(1, 36, 36)\n", - "\n", - "plt_plot_len_x_might_be_one(ds1_sep, x1, x2)\n", - "plt_plot_len_x_might_be_one(ds2_sep, x2, x1)\n", - "plt.plot(x3, sep_area[0:36])\n", + "cesm1_aicetot_nh_sep = ds_cesm1_aicetot_nh[\"aice_monthly\"].isel(nmonth=8)\n", + "cesm2_aicetot_nh_sep = ds_cesm2_aicetot_nh[\"aice_monthly\"].isel(nmonth=8)\n", + "\n", + "x1 = np.linspace(end_year - len(ds1_sep.time), end_year, len(ds1_sep.time))\n", + "x2 = np.linspace(end_year - len(ds2_sep.time), base_end_year, len(ds2_sep.time))\n", + "x3 = np.linspace(ds_cesm1_aicetot_nh.year[60], ds_cesm1_aicetot_nh.year[95], 36)\n", + "x4 = np.linspace(ds_cesm2_aicetot_nh.year[110], ds_cesm2_aicetot_nh.year[145], 36)\n", + "obs_first_year = 0\n", + "if first_year > 1:\n", + " obs_first_year = 1979\n", + "x5 = np.linspace(1, 36, 36) + obs_first_year\n", + "\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " x3,\n", + " cesm1_aicetot_nh_sep.isel(n_members=i, nyr=slice(60, 96)) * 1.0e-12,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " x4,\n", + " cesm2_aicetot_nh_sep.isel(n_members=i, nyr=slice(110, 146)) * 1.0e-12,\n", + " color=\"lightblue\",\n", + " )\n", + "plt_plot_len_x_might_be_one(ds1_sep, x1, x2, color=\"red\")\n", + "plt_plot_len_x_might_be_one(ds2_sep, x2, x1, color=\"blue\")\n", + "plt.plot(x5, sep_area[0:36], color=\"black\")\n", "\n", - "plt.title(\"NH Annual Mean Sea Ice Area\")\n", + "plt.title(\"NH September Mean Sea Ice Area\")\n", "plt.ylim((0, 25))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Sea Ice Area $mx10^{12}$\")\n", - "plt.legend([case_name, base_case_name, \"NSIDC\"])" + "plt.legend(handles=[p1, p2, p3, p4, p5])" ] }, { @@ -700,19 +901,41 @@ "ds1_sep_sh = ds1_area_sh.sel(time=(ds1_area_sh.time.dt.month == 9))\n", "ds2_sep_sh = ds2_area_sh.sel(time=(ds2_area_sh.time.dt.month == 9))\n", "\n", - "x1 = np.linspace(1, len(ds1_sep.time), len(ds1_sep.time))\n", - "x2 = np.linspace(1, len(ds2_sep.time), len(ds2_sep.time))\n", - "x3 = np.linspace(1, 36, 36)\n", + "cesm1_aicetot_sh_sep = ds_cesm1_aicetot_sh[\"aice_monthly\"].isel(nmonth=8)\n", + "cesm2_aicetot_sh_sep = ds_cesm2_aicetot_sh[\"aice_monthly\"].isel(nmonth=8)\n", + "\n", + "x1 = np.linspace(end_year - len(ds1_sep.time), end_year, len(ds1_sep.time))\n", + "x2 = np.linspace(end_year - len(ds2_sep.time), base_end_year, len(ds2_sep.time))\n", + "x3 = np.linspace(ds_cesm1_aicetot_nh.year[60], ds_cesm1_aicetot_nh.year[95], 36)\n", + "x4 = np.linspace(ds_cesm2_aicetot_nh.year[110], ds_cesm2_aicetot_nh.year[145], 36)\n", + "\n", "\n", - "plt_plot_len_x_might_be_one(ds1_sep_sh, x1, x2)\n", - "plt_plot_len_x_might_be_one(ds2_sep_sh, x2, x1)\n", - "plt.plot(x3, sep_area_sh[0:36])\n", + "obs_first_year = 0\n", + "if first_year > 1:\n", + " obs_first_year = 1979\n", + "x5 = np.linspace(1, 36, 36) + obs_first_year\n", + "\n", + "for i in range(0, 38):\n", + " plt.plot(\n", + " x3,\n", + " cesm1_aicetot_sh_sep.isel(n_members=i, nyr=slice(60, 96)) * 1.0e-12,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 49):\n", + " plt.plot(\n", + " x4,\n", + " cesm2_aicetot_sh_sep.isel(n_members=i, nyr=slice(110, 146)) * 1.0e-12,\n", + " color=\"lightblue\",\n", + " )\n", + "plt_plot_len_x_might_be_one(ds1_sep_sh, x1, x2, color=\"red\")\n", + "plt_plot_len_x_might_be_one(ds2_sep_sh, x2, x1, color=\"blue\")\n", + "plt.plot(x5, sep_area_sh[0:36], color=\"black\")\n", "\n", - "plt.title(\"SH Annual Mean Sea Ice Area\")\n", + "plt.title(\"SH September Mean Sea Ice Area\")\n", "plt.ylim((0, 25))\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Sea Ice Area $mx10^{12}$\")\n", - "plt.legend([case_name, base_case_name, \"NSIDC\"])" + "plt.legend(handles=[p1, p2, p3, p4, p5])" ] }, { @@ -741,8 +964,8 @@ "ds1_lab = (mask * tarea * ds1.aice).sum(dim=[\"nj\", \"ni\"]) * 1.0e-12\n", "ds2_lab = (mask * tarea * ds2.aice).sum(dim=[\"nj\", \"ni\"]) * 1.0e-12\n", "\n", - "ds1_lab.plot()\n", - "ds2_lab.plot()\n", + "ds1_lab.plot(color=\"red\")\n", + "ds2_lab.plot(color=\"blue\")\n", "\n", "plt.title(\"Labrador Sea Montly Mean Sea Ice Area\")\n", "plt.ylim((0, 10))\n", diff --git a/nblibrary/ice/plot_diff.py b/nblibrary/ice/plot_diff.py index fa5dd87..4fc1701 100644 --- a/nblibrary/ice/plot_diff.py +++ b/nblibrary/ice/plot_diff.py @@ -6,9 +6,73 @@ import matplotlib.path as mpath import matplotlib.pyplot as plt import numpy as np +import xarray as xr from matplotlib.gridspec import GridSpec +def add_cyclic(ds): + + ni = ds.tlon.shape[1] + + xL = int(ni / 2 - 1) + xR = int(xL + ni) + + tlon = ds.tlon.data + tlat = ds.tlat.data + + tlon = np.where(np.greater_equal(tlon, min(tlon[:, 0])), tlon - 360.0, tlon) + lon = np.concatenate((tlon, tlon + 360.0), 1) + lon = lon[:, xL:xR] + + if ni == 320: + lon[367:-3, 0] = lon[367:-3, 0] + 360.0 + lon = lon - 360.0 + + lon = np.hstack((lon, lon[:, 0:1] + 360.0)) + if ni == 320: + lon[367:, -1] = lon[367:, -1] - 360.0 + + # -- trick cartopy into doing the right thing: + # it gets confused when the cyclic coords are identical + lon[:, 0] = lon[:, 0] - 1e-8 + + # -- periodicity + lat = np.concatenate((tlat, tlat), 1) + lat = lat[:, xL:xR] + lat = np.hstack((lat, lat[:, 0:1])) + + TLAT = xr.DataArray(lat, dims=("nlat", "nlon")) + TLONG = xr.DataArray(lon, dims=("nlat", "nlon")) + + dso = xr.Dataset({"TLAT": TLAT, "TLONG": TLONG}) + # copy vars + varlist = [v for v in ds.data_vars if v not in ["TLAT", "TLONG"]] + for v in varlist: + v_dims = ds[v].dims + if not ("nlat" in v_dims and "nlon" in v_dims): + dso[v] = ds[v] + else: + # determine and sort other dimensions + other_dims = set(v_dims) - {"nlat", "nlon"} + other_dims = tuple([d for d in v_dims if d in other_dims]) + lon_dim = ds[v].dims.index("nlon") + field = ds[v].data + field = np.concatenate((field, field), lon_dim) + field = field[..., :, xL:xR] + field = np.concatenate((field, field[..., :, 0:1]), lon_dim) + dso[v] = xr.DataArray( + field, + dims=other_dims + ("nlat", "nlon"), + attrs=ds[v].attrs, + ) + # copy coords + for v, da in ds.coords.items(): + if not ("nlat" in da.dims and "nlon" in da.dims): + dso = dso.assign_coords(**{v: da}) + + return dso + + def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): # make circular boundary for polar stereographic circular plots theta = np.linspace(0, 2 * np.pi, 100) @@ -16,6 +80,16 @@ def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): verts = np.vstack([np.sin(theta), np.cos(theta)]).T circle = mpath.Path(verts * radius + center) + # Read in observed sea ice concentration + + path_nsidc = "/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_obs_data/ice/" + + ds_obs = xr.open_dataset(path_nsidc + "SSMI.ifrac.1981-2005monthlymean.gx1v5.nc") + + ds_pop = add_cyclic(ds_obs) + + aice = title.find("Concentration") + if np.size(levels) > 2: cmap = mpl.colormaps["ocean"] norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N) @@ -28,10 +102,12 @@ def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): ax = fig.add_subplot(gs[0, :2], projection=ccrs.NorthPolarStereo()) # sets the latitude / longitude boundaries of the plot ax.set_extent([0.005, 360, 90, 45], crs=ccrs.PlateCarree()) + ifrac_obs = ds_pop["monthly_ifrac"].isel(month=2) if proj == "S": ax = fig.add_subplot(gs[0, :2], projection=ccrs.SouthPolarStereo()) # sets the latitude / longitude boundaries of the plot ax.set_extent([0.005, 360, -90, -45], crs=ccrs.PlateCarree()) + ifrac_obs = ds_pop["monthly_ifrac"].isel(month=8) ax.set_boundary(circle, transform=ax.transAxes) ax.add_feature(cfeature.LAND, zorder=100, edgecolor="k") @@ -47,6 +123,15 @@ def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): cmap="ocean", transform=ccrs.PlateCarree(), ) + if aice > 0: + plt.contour( + ds_pop.tlon, + ds_pop.tlat, + ifrac_obs, + levels=[0.15], + colors="magenta", + transform=ccrs.PlateCarree(), + ) plt.colorbar(this, orientation="vertical", fraction=0.04, pad=0.01) plt.title(case1, fontsize=10) @@ -70,6 +155,15 @@ def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): cmap="ocean", transform=ccrs.PlateCarree(), ) + if aice > 0: + plt.contour( + ds_pop.tlon, + ds_pop.tlat, + ifrac_obs, + levels=[0.15], + colors="magenta", + transform=ccrs.PlateCarree(), + ) plt.colorbar(this, orientation="vertical", fraction=0.04, pad=0.01) plt.title(case2, fontsize=10)