From a4bd9ca44ecac23cd02afbf96a25825c085ed562 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Mon, 3 Feb 2025 14:42:50 -0700 Subject: [PATCH 1/3] Add in CESM2 and CESM1 LENS data and fix axes in sea ice plots. Also add contours to sea ice concentration maps. Still a work in progress. --- .../Hemis_seaice_visual_compare_contour.ipynb | 13 +- ...Hemis_seaice_visual_compare_obs_lens.ipynb | 323 +++++++++++++----- nblibrary/ice/plot_diff.py | 95 ++++++ 3 files changed, 345 insertions(+), 86 deletions(-) diff --git a/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb b/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb index a7a9f610..a3d924f4 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 = {}" ] }, { @@ -86,7 +87,7 @@ }, "outputs": [], "source": [ - "cupid_run = 1\n", + "cupid_run = 0\n", "\n", "if cupid_run == 1:\n", "\n", @@ -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 33dedb23..ef22bc00 100644 --- a/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb +++ b/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb @@ -38,6 +38,7 @@ "import yaml\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", + "import matplotlib.lines as mlines\n", "import nc_time_axis" ] }, @@ -72,9 +73,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,42 +179,72 @@ "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", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [ - "hide-input" - ] - }, + "id": "31c8d5bf-b76f-46b2-bb5b-a654fa2150ed", + "metadata": {}, "outputs": [], "source": [ - "### Read in the CESM1 and CESM2 LE sea ice area\n", + "### Read in the CESM LENS historical data\n", "\n", - "# cesm1_path = \"/glade/campaign/cesm/collections/cesmLE/CESM-CAM5-BGC-LE/ice/proc/tseries/monthly/aice/\"\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_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_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", - "# 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", + "cesm1_years = xr.cftime_range(start=\"1920\", end=\"2100\", freq=\"YS\", calendar=\"noleap\")\n", "\n", - "# nj_nh = ds_cesm1_nh['nj']\n", - "# nj_sh = ds_cesm2_sh['nh']\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", + "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_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_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", + "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\")" ] }, { @@ -225,21 +258,26 @@ "# 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 +304,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", + " ds_cesm1_hitot_nh[\"hi_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", + " ds_cesm1_hstot_nh[\"hs_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", + " ds_cesm1_aicetot_nh[\"aice_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 +397,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", + " ds_cesm1_hitot_sh[\"hi_ann\"][i, 0:95] * 1.0e-13,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 39):\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", + " ds_cesm1_hstot_sh[\"hs_ann\"][i, 0:95] * 1.0e-13,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 39):\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", + " ds_cesm1_aicetot_sh[\"aice_ann\"][i, 0:95] * 1.0e-12,\n", + " color=\"lightgrey\",\n", + " )\n", + "for i in range(0, 39):\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 +535,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 +550,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 +648,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 +663,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 +713,34 @@ "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", + "for i in range(0, 38):\n", + " plt.plot(cesm1_aicetot_nh_month[i, :] * 1.0e-12, color=\"lightgrey\")\n", + "for i in range(0, 49):\n", + " plt.plot(cesm2_aicetot_nh_month[i, :] * 1.0e-12, color=\"lightblue\")\n", + "plt.plot(area1, color=\"red\")\n", + "plt.plot(area2, color=\"blue\")\n", + "plt.plot(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 +764,29 @@ "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(cesm1_aicetot_sh_month[i, :] * 1.0e-12, color=\"lightgrey\")\n", + "for i in range(0, 39):\n", + " plt.plot(cesm2_aicetot_sh_month[i, :] * 1.0e-12, color=\"lightblue\")\n", + "plt.plot(area1, color=\"red\")\n", + "plt.plot(area1, color=\"red\")\n", + "plt.plot(area2, color=\"blue\")\n", + "plt.plot(nsidc_clim_sh_area, color=\"black\")\n", "\n", "plt.title(\"SH 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])" ] }, { @@ -657,19 +814,22 @@ "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", + "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", + "obs_first_year = 0\n", + "if first_year > 1:\n", + " obs_first_year = 1979\n", + "x3 = np.linspace(1, 36, 36) + obs_first_year\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", + "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(x3, 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 +860,22 @@ "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", + "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", + "obs_first_year = 0\n", + "if first_year > 1:\n", + " obs_first_year = 1979\n", + "x3 = np.linspace(1, 36, 36) + obs_first_year\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", + "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(x3, 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])" ] }, { diff --git a/nblibrary/ice/plot_diff.py b/nblibrary/ice/plot_diff.py index fa5dd87e..30a8f0da 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,18 @@ 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) + + ifrac_obs = ds_pop.monthly_ifrac.mean(dim="month") + + aice = title.find("Concentration") + if np.size(levels) > 2: cmap = mpl.colormaps["ocean"] norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N) @@ -25,6 +101,7 @@ def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): gs = GridSpec(2, 4) if proj == "N": + 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()) @@ -47,6 +124,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 +156,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) From ca3af951f0c8336f06f4d30bb447e59d9f24aefa Mon Sep 17 00:00:00 2001 From: David Bailey Date: Mon, 3 Feb 2025 18:15:34 -0700 Subject: [PATCH 2/3] Adding LENS to more plots and fixing some axes --- ...Hemis_seaice_visual_compare_obs_lens.ipynb | 101 +++++++++++++----- 1 file changed, 73 insertions(+), 28 deletions(-) diff --git a/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb b/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb index ef22bc00..bc202c1a 100644 --- a/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb +++ b/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb @@ -204,6 +204,14 @@ "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", @@ -315,7 +323,7 @@ "for i in range(0, 38):\n", " plt.plot(\n", " cesm1_years[0:95],\n", - " ds_cesm1_hitot_nh[\"hi_ann\"][i, 0:95] * 1.0e-13,\n", + " cesm1_hitot_nh_ann[i, 0:95] * 1.0e-13,\n", " color=\"lightgrey\",\n", " )\n", "for i in range(0, 49):\n", @@ -335,7 +343,7 @@ "for i in range(0, 38):\n", " plt.plot(\n", " cesm1_years[0:95],\n", - " ds_cesm1_hstot_nh[\"hs_ann\"][i, 0:95] * 1.0e-13,\n", + " cesm1_hstot_nh_ann[i, 0:95] * 1.0e-13,\n", " color=\"lightgrey\",\n", " )\n", "for i in range(0, 49):\n", @@ -354,7 +362,7 @@ "for i in range(0, 38):\n", " plt.plot(\n", " cesm1_years[0:95],\n", - " ds_cesm1_aicetot_nh[\"aice_ann\"][i, 0:95] * 1.0e-12,\n", + " cesm1_aicetot_nh_ann[i, 0:95] * 1.0e-12,\n", " color=\"lightgrey\",\n", " )\n", "for i in range(0, 49):\n", @@ -400,10 +408,10 @@ "for i in range(0, 38):\n", " plt.plot(\n", " cesm1_years[0:95],\n", - " ds_cesm1_hitot_sh[\"hi_ann\"][i, 0:95] * 1.0e-13,\n", + " cesm1_hitot_sh_ann[i, 0:95] * 1.0e-13,\n", " color=\"lightgrey\",\n", " )\n", - "for i in range(0, 39):\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", @@ -420,10 +428,10 @@ "for i in range(0, 38):\n", " plt.plot(\n", " cesm1_years[0:95],\n", - " ds_cesm1_hstot_sh[\"hs_ann\"][i, 0:95] * 1.0e-13,\n", + " cesm1_hstot_sh_ann[i, 0:95] * 1.0e-13,\n", " color=\"lightgrey\",\n", " )\n", - "for i in range(0, 39):\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", @@ -439,10 +447,10 @@ "for i in range(0, 38):\n", " plt.plot(\n", " cesm1_years[0:95],\n", - " ds_cesm1_aicetot_sh[\"aice_ann\"][i, 0:95] * 1.0e-12,\n", + " cesm1_aicetot_sh_ann[i, 0:95] * 1.0e-12,\n", " color=\"lightgrey\",\n", " )\n", - "for i in range(0, 39):\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", @@ -728,13 +736,15 @@ "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(cesm1_aicetot_nh_month[i, :] * 1.0e-12, color=\"lightgrey\")\n", + " plt.plot(months, cesm1_aicetot_nh_month[i, :] * 1.0e-12, color=\"lightgrey\")\n", "for i in range(0, 49):\n", - " plt.plot(cesm2_aicetot_nh_month[i, :] * 1.0e-12, color=\"lightblue\")\n", - "plt.plot(area1, color=\"red\")\n", - "plt.plot(area2, color=\"blue\")\n", - "plt.plot(nsidc_clim_area, color=\"black\")\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", @@ -774,16 +784,15 @@ "area2 = (aice2_month * mask_sh).sum([\"ni\", \"nj\"]) * 1.0e-12\n", "\n", "for i in range(0, 38):\n", - " plt.plot(cesm1_aicetot_sh_month[i, :] * 1.0e-12, color=\"lightgrey\")\n", - "for i in range(0, 39):\n", - " plt.plot(cesm2_aicetot_sh_month[i, :] * 1.0e-12, color=\"lightblue\")\n", - "plt.plot(area1, color=\"red\")\n", - "plt.plot(area1, color=\"red\")\n", - "plt.plot(area2, color=\"blue\")\n", - "plt.plot(nsidc_clim_sh_area, color=\"black\")\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 Area $m x 10^{12}$\")\n", "plt.legend(handles=[p1, p2, p3, p4, p5])" @@ -814,16 +823,33 @@ "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", + "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", - "x3 = np.linspace(1, 36, 36) + obs_first_year\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(x3, sep_area[0:36], color=\"black\")\n", + "plt.plot(x5, sep_area[0:36], color=\"black\")\n", "\n", "plt.title(\"NH September Mean Sea Ice Area\")\n", "plt.ylim((0, 25))\n", @@ -860,16 +886,35 @@ "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", + "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", "obs_first_year = 0\n", "if first_year > 1:\n", " obs_first_year = 1979\n", - "x3 = np.linspace(1, 36, 36) + obs_first_year\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(x3, sep_area_sh[0:36], color=\"black\")\n", + "plt.plot(x5, sep_area_sh[0:36], color=\"black\")\n", "\n", "plt.title(\"SH September Mean Sea Ice Area\")\n", "plt.ylim((0, 25))\n", @@ -904,8 +949,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", From 5345b7b31d747591d833c42d1191a8e2bbf23890 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Tue, 4 Feb 2025 09:38:24 -0700 Subject: [PATCH 3/3] Remove nc_time_axis and fix hide-input for cells --- .../Hemis_seaice_visual_compare_contour.ipynb | 2 +- ...Hemis_seaice_visual_compare_obs_lens.ipynb | 23 +++++++++++++++---- nblibrary/ice/plot_diff.py | 5 ++-- 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb b/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb index a3d924f4..f7ca899f 100644 --- a/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb +++ b/nblibrary/ice/Hemis_seaice_visual_compare_contour.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "cupid_run = 0\n", + "cupid_run = 1\n", "\n", "if cupid_run == 1:\n", "\n", diff --git a/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb b/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb index bc202c1a..7a9e0e24 100644 --- a/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb +++ b/nblibrary/ice/Hemis_seaice_visual_compare_obs_lens.ipynb @@ -38,8 +38,7 @@ "import yaml\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", - "import matplotlib.lines as mlines\n", - "import nc_time_axis" + "import matplotlib.lines as mlines" ] }, { @@ -191,7 +190,15 @@ "cell_type": "code", "execution_count": null, "id": "31c8d5bf-b76f-46b2-bb5b-a654fa2150ed", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, "outputs": [], "source": [ "### Read in the CESM LENS historical data\n", @@ -259,7 +266,15 @@ "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", diff --git a/nblibrary/ice/plot_diff.py b/nblibrary/ice/plot_diff.py index 30a8f0da..4fc17011 100644 --- a/nblibrary/ice/plot_diff.py +++ b/nblibrary/ice/plot_diff.py @@ -88,8 +88,6 @@ def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): ds_pop = add_cyclic(ds_obs) - ifrac_obs = ds_pop.monthly_ifrac.mean(dim="month") - aice = title.find("Concentration") if np.size(levels) > 2: @@ -101,14 +99,15 @@ def plot_diff(field1, field2, levels, case1, case2, title, proj, TLAT, TLON): gs = GridSpec(2, 4) if proj == "N": - 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")