From 22d45dafa88d6140cbff623dc0a9349c70c2bfbc Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 17 May 2024 17:01:00 -0700 Subject: [PATCH 01/14] add nc code --- pcmdi_metrics/sea_ice/sea_ice_driver.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/pcmdi_metrics/sea_ice/sea_ice_driver.py b/pcmdi_metrics/sea_ice/sea_ice_driver.py index 64471a6c7..4f33c3ce8 100644 --- a/pcmdi_metrics/sea_ice/sea_ice_driver.py +++ b/pcmdi_metrics/sea_ice/sea_ice_driver.py @@ -49,6 +49,7 @@ oeyear = parameter.oeyear plot = parameter.plot pole = parameter.pole + to_nc = parameter.netcdf print("Model list:", model_list) model_list.sort() @@ -358,6 +359,22 @@ mask = create_land_sea_mask(ds, lon_key=xvar, lat_key=yvar) ds[var] = ds[var].where(mask < 1) + if to_nc: + # Generate netcdf files of climatologies + # TODO: this requires a refactoring of the get_clim function + # to accept valid datasets as well as data arrays + nc_dir = os.path.join(metrics_output_path, "netcdf") + if not os.path.exists(nc_dir): + os.mkdir(nc_dir) + nc_climo = lib.get_clim(ds, var, ds=None) + fname = ( + "sic_climatology_" + + "_".join(model, run, yr_range[0], yr_range[1]) + + ".nc" + ) + fname = os.path.join(nc_dir, fname) + nc_climo.to_netcdf(fname) + # Get regions clims, means = lib.process_by_region(ds, var, area[area_var].data, pole) From ac2ef421ec1b36d0a3878481b03406f4df338bcc Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 20 May 2024 16:10:04 -0700 Subject: [PATCH 02/14] add nc opt --- pcmdi_metrics/sea_ice/lib/sea_ice_parser.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_parser.py b/pcmdi_metrics/sea_ice/lib/sea_ice_parser.py index 4d1fe9806..7613a4036 100644 --- a/pcmdi_metrics/sea_ice/lib/sea_ice_parser.py +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_parser.py @@ -230,4 +230,11 @@ def create_sea_ice_parser(): default=90.1, help="Set to a latitude value to exclude sea ice at North pole. Must be > 80.", ) + + parser.add_argument( + "--netcdf", + action="store_true", + default=False, + help="Use to save netcdf intermediate results", + ) return parser From cd71c71821895fe7218c78c0b10c9e4460d05321 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 20 May 2024 16:10:26 -0700 Subject: [PATCH 03/14] refactor clim --- pcmdi_metrics/sea_ice/lib/sea_ice_lib.py | 29 ++++++++++++++---------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py b/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py index 759b39e3b..9dcfc7259 100644 --- a/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py @@ -222,17 +222,22 @@ def get_total_extent(data, ds_area): return total_extent, te_mean -def get_clim(total_extent, ds, ds_var): - try: - clim = to_ice_con_ds(total_extent, ds, ds_var).temporal.climatology( - ds_var, freq="month" - ) - except IndexError: # Issue with time bounds - tmp = to_ice_con_ds(total_extent, ds, ds_var) - tbkey = xcdat_dataset_io.get_time_bounds_key(tmp) - tmp = tmp.drop_vars(tbkey) - tmp = tmp.bounds.add_missing_bounds() - clim = tmp.temporal.climatology(ds_var, freq="month") +def get_clim(total_extent, ds_var, ds=None): + # ds is a dataset that contains the dimensions + # needed to turn total_extent into a dataset + if ds is None: + clim = total_extent.temporal.climatology(ds_var, freq="month") + else: + try: + clim = to_ice_con_ds(total_extent, ds, ds_var).temporal.climatology( + ds_var, freq="month" + ) + except IndexError: # Issue with time bounds + tmp = to_ice_con_ds(total_extent, ds, ds_var) + tbkey = xcdat_dataset_io.get_time_bounds_key(tmp) + tmp = tmp.drop_vars(tbkey) + tmp = tmp.bounds.add_missing_bounds() + clim = tmp.temporal.climatology(ds_var, freq="month") return clim @@ -245,7 +250,7 @@ def process_by_region(ds, ds_var, ds_area, pole): yvar = find_lat(ds) data = choose_region(region, ds, ds_var, xvar, yvar, pole) total_extent, te_mean = get_total_extent(data, ds_area) - clim = get_clim(total_extent, ds, ds_var) + clim = get_clim(total_extent, ds_var, ds) clims[region] = clim means[region] = te_mean del data From 90c682c2e7a77abc0a189ea15775e59ff0c3b0b0 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 20 May 2024 16:10:46 -0700 Subject: [PATCH 04/14] add print --- pcmdi_metrics/sea_ice/sea_ice_driver.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pcmdi_metrics/sea_ice/sea_ice_driver.py b/pcmdi_metrics/sea_ice/sea_ice_driver.py index 4f33c3ce8..6c8609def 100644 --- a/pcmdi_metrics/sea_ice/sea_ice_driver.py +++ b/pcmdi_metrics/sea_ice/sea_ice_driver.py @@ -363,13 +363,15 @@ # Generate netcdf files of climatologies # TODO: this requires a refactoring of the get_clim function # to accept valid datasets as well as data arrays + print("Generating climatology for netcdf") nc_dir = os.path.join(metrics_output_path, "netcdf") if not os.path.exists(nc_dir): os.mkdir(nc_dir) nc_climo = lib.get_clim(ds, var, ds=None) + print("Writing climatology netcdf") fname = ( - "sic_climatology_" - + "_".join(model, run, yr_range[0], yr_range[1]) + "sic_clim_" + + "_".join([model, run, yr_range[0], yr_range[1]]) + ".nc" ) fname = os.path.join(nc_dir, fname) From 0a7feff32c290e235f44bd9a862050680dc2a7dd Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 24 May 2024 13:33:03 -0700 Subject: [PATCH 05/14] add figs --- pcmdi_metrics/sea_ice/sea_ice_driver.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/pcmdi_metrics/sea_ice/sea_ice_driver.py b/pcmdi_metrics/sea_ice/sea_ice_driver.py index 6c8609def..29a3c9e1c 100644 --- a/pcmdi_metrics/sea_ice/sea_ice_driver.py +++ b/pcmdi_metrics/sea_ice/sea_ice_driver.py @@ -11,6 +11,7 @@ from pcmdi_metrics.io.base import Base from pcmdi_metrics.sea_ice.lib import create_sea_ice_parser +from pcmdi_metrics.sea_ice.lib import sea_ice_figures as fig from pcmdi_metrics.sea_ice.lib import sea_ice_lib as lib from pcmdi_metrics.utils import create_land_sea_mask @@ -361,8 +362,6 @@ if to_nc: # Generate netcdf files of climatologies - # TODO: this requires a refactoring of the get_clim function - # to accept valid datasets as well as data arrays print("Generating climatology for netcdf") nc_dir = os.path.join(metrics_output_path, "netcdf") if not os.path.exists(nc_dir): @@ -377,6 +376,27 @@ fname = os.path.join(nc_dir, fname) nc_climo.to_netcdf(fname) + # Filling in NaNs with zeros for prettier plot + nc_climo = nc_climo.fillna(0) + nc_climo[var] = nc_climo[var].where(mask < 1) + fig_dir = os.path.join(metrics_output_path, "plot") + if not os.path.exists(fig_dir): + os.mkdir(fig_dir) + for mo in range(0, 12): + tmp_title = "_".join([model, str(mo + 1)]) + fig.create_arctic_map( + nc_climo.isel({"time": mo}), + fig_dir, + varname="siconc", + title=tmp_title, + ) + fig.create_antarctic_map( + nc_climo.isel({"time": mo}), + fig_dir, + varname="siconc", + title=tmp_title, + ) + # Get regions clims, means = lib.process_by_region(ds, var, area[area_var].data, pole) From b09fa33cc88e6ae2ad7d65a465b41259385146aa Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 24 May 2024 13:34:15 -0700 Subject: [PATCH 06/14] add figs --- pcmdi_metrics/sea_ice/lib/sea_ice_figures.py | 165 +++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 pcmdi_metrics/sea_ice/lib/sea_ice_figures.py diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py new file mode 100644 index 000000000..a64d94a6d --- /dev/null +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py @@ -0,0 +1,165 @@ +import os + +import cartopy.crs as ccrs +import matplotlib.colors as colors +import matplotlib.pyplot as plt +import numpy as np +import regionmask + +from pcmdi_metrics.sea_ice.lib import sea_ice_lib as lib + + +def create_arctic_map(ds, metrics_output_path, varname="siconc", title=None): + print("Creating Arctic map") + # Load and process data + xvar = lib.find_lon(ds) + yvar = lib.find_lat(ds) + + # Set up regions + region_NA = np.array([[-120, 45], [-120, 80], [90, 80], [90, 45]]) + region_NP = np.array([[90, 45], [90, 65], [240, 65], [240, 45]]) + names = ["North_Atlantic", "North_Pacific"] + abbrevs = ["NA", "NP"] + arctic_regions = regionmask.Regions( + [region_NA, region_NP], names=names, abbrevs=abbrevs, name="arctic" + ) + + # Do plotting + cmap = colors.LinearSegmentedColormap.from_list( + "", [[0, 85 / 255, 182 / 255], "white"] + ) + proj = ccrs.NorthPolarStereo() + ax = plt.subplot(111, projection=proj) + ax.set_global() + # TODO: get xvar and yvar + ds[varname].plot.pcolormesh( + ax=ax, + x=xvar, + y=yvar, + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice concentration (%)"}, + ) + arctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + plt.annotate( + "North Atlantic", + (0.5, 0.2), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="white", + ) + plt.annotate( + "North Pacific", + (0.65, 0.88), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="white", + ) + plt.annotate( + "Central\nArctic ", + (0.56, 0.56), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + ) + ax.set_facecolor([0.55, 0.55, 0.6]) + if title is not None: + plt.title(title) + fig_path = os.path.join( + metrics_output_path, title.replace(" ", "_") + "_arctic_regions.png" + ) + else: + fig_path = os.path.join(metrics_output_path, "arctic_regions.png") + plt.savefig(fig_path) + plt.close() + + +# ---------- +# Antarctic +# ---------- +def create_antarctic_map(ds, metrics_output_path, varname="siconc", title=None): + print("Creating Antarctic map") + # Load and process data + xvar = lib.find_lon(ds) + yvar = lib.find_lat(ds) + + # Set up regions + region_IO = np.array([[20, -90], [90, -90], [90, -55], [20, -55]]) + region_SA = np.array([[20, -90], [-60, -90], [-60, -55], [20, -55]]) + region_SP = np.array([[90, -90], [300, -90], [300, -55], [90, -55]]) + + names = ["Indian Ocean", "South Atlantic", "South Pacific"] + abbrevs = ["IO", "SA", "SP"] + arctic_regions = regionmask.Regions( + [region_IO, region_SA, region_SP], + names=names, + abbrevs=abbrevs, + name="antarctic", + ) + + # Do plotting + cmap = colors.LinearSegmentedColormap.from_list( + "", [[0, 85 / 255, 182 / 255], "white"] + ) + proj = ccrs.SouthPolarStereo() + ax = plt.subplot(111, projection=proj) + ax.set_global() + ds[varname].plot.pcolormesh( + ax=ax, + x=xvar, + y=yvar, + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice concentration (%)"}, + ) + arctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + plt.annotate( + "South Pacific", + (0.50, 0.2), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="black", + ) + plt.annotate( + "Indian\nOcean", + (0.89, 0.66), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="black", + ) + plt.annotate( + "South Atlantic", + (0.54, 0.82), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="black", + ) + ax.set_facecolor([0.55, 0.55, 0.6]) + if title is not None: + plt.title(title) + fig_path = os.path.join( + metrics_output_path, title.replace(" ", "_") + "_antarctic_regions.png" + ) + else: + fig_path = os.path.join(metrics_output_path, "antarctic_regions.png") + plt.savefig(fig_path) + plt.close() From a60fe93f4e911ba72210de193f5f36887020c1ed Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Thu, 6 Jun 2024 10:37:35 -0700 Subject: [PATCH 07/14] work in progress --- pcmdi_metrics/sea_ice/lib/sea_ice_figures.py | 25 ++++++++++++++++++-- pcmdi_metrics/sea_ice/sea_ice_driver.py | 2 +- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py index a64d94a6d..9238580f9 100644 --- a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py @@ -5,16 +5,32 @@ import matplotlib.pyplot as plt import numpy as np import regionmask +import xarray as xr from pcmdi_metrics.sea_ice.lib import sea_ice_lib as lib +def replace_nan_zero(da): + da_new = xr.where(np.isnan(da), 0, da) + return da_new + + +def replace_fill_zero(da, val): + da_new = xr.where(da > val, 0, da) + return da_new + + def create_arctic_map(ds, metrics_output_path, varname="siconc", title=None): print("Creating Arctic map") # Load and process data xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) + # Some models have NaN values in coordinates + # that can't be plotted by pcolormesh + ds[xvar] = replace_nan_zero(ds[xvar]) + ds[yvar] = replace_nan_zero(ds[yvar]) + # Set up regions region_NA = np.array([[-120, 45], [-120, 80], [90, 80], [90, 45]]) region_NP = np.array([[90, 45], [90, 65], [240, 65], [240, 45]]) @@ -38,7 +54,7 @@ def create_arctic_map(ds, metrics_output_path, varname="siconc", title=None): y=yvar, transform=ccrs.PlateCarree(), cmap=cmap, - cbar_kwargs={"label": "ice concentration (%)"}, + cbar_kwargs={"label": "ice fraction"}, ) arctic_regions.plot_regions( ax=ax, @@ -92,6 +108,11 @@ def create_antarctic_map(ds, metrics_output_path, varname="siconc", title=None): xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) + # Some models have NaN values in coordinates + # that can't be plotted by pcolormesh + ds[xvar] = replace_nan_zero(ds[xvar]) + ds[yvar] = replace_nan_zero(ds[yvar]) + # Set up regions region_IO = np.array([[20, -90], [90, -90], [90, -55], [20, -55]]) region_SA = np.array([[20, -90], [-60, -90], [-60, -55], [20, -55]]) @@ -119,7 +140,7 @@ def create_antarctic_map(ds, metrics_output_path, varname="siconc", title=None): y=yvar, transform=ccrs.PlateCarree(), cmap=cmap, - cbar_kwargs={"label": "ice concentration (%)"}, + cbar_kwargs={"label": "ice fraction"}, ) arctic_regions.plot_regions( ax=ax, diff --git a/pcmdi_metrics/sea_ice/sea_ice_driver.py b/pcmdi_metrics/sea_ice/sea_ice_driver.py index 204a589ea..4e4df4887 100644 --- a/pcmdi_metrics/sea_ice/sea_ice_driver.py +++ b/pcmdi_metrics/sea_ice/sea_ice_driver.py @@ -391,7 +391,7 @@ if not os.path.exists(fig_dir): os.mkdir(fig_dir) for mo in range(0, 12): - tmp_title = "_".join([model, str(mo + 1)]) + tmp_title = "_".join([model, run, str(mo + 1)]) fig.create_arctic_map( nc_climo.isel({"time": mo}), fig_dir, From 8ca0c7dc7e059fd574c442b6adeecdc43d16000f Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 7 Jun 2024 14:27:50 -0700 Subject: [PATCH 08/14] refactor clim --- pcmdi_metrics/sea_ice/lib/sea_ice_lib.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py b/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py index 9dcfc7259..85054da97 100644 --- a/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_lib.py @@ -226,18 +226,16 @@ def get_clim(total_extent, ds_var, ds=None): # ds is a dataset that contains the dimensions # needed to turn total_extent into a dataset if ds is None: - clim = total_extent.temporal.climatology(ds_var, freq="month") + ds_new = total_extent else: - try: - clim = to_ice_con_ds(total_extent, ds, ds_var).temporal.climatology( - ds_var, freq="month" - ) - except IndexError: # Issue with time bounds - tmp = to_ice_con_ds(total_extent, ds, ds_var) - tbkey = xcdat_dataset_io.get_time_bounds_key(tmp) - tmp = tmp.drop_vars(tbkey) - tmp = tmp.bounds.add_missing_bounds() - clim = tmp.temporal.climatology(ds_var, freq="month") + ds_new = to_ice_con_ds(total_extent, ds, ds_var) + try: + clim = ds_new.temporal.climatology(ds_var, freq="month") + except IndexError: # Issue with time bounds + tbkey = xcdat_dataset_io.get_time_bounds_key(ds_new) + ds_new = ds_new.drop_vars(tbkey) + ds_new = ds_new.bounds.add_missing_bounds() + clim = ds_new.temporal.climatology(ds_var, freq="month") return clim From 33be6bf44e34ccd071eb312b84c14fe57ed30cb5 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 18 Jun 2024 15:17:12 -0700 Subject: [PATCH 09/14] update driver --- pcmdi_metrics/sea_ice/sea_ice_driver.py | 259 ++++++++++-------------- 1 file changed, 111 insertions(+), 148 deletions(-) diff --git a/pcmdi_metrics/sea_ice/sea_ice_driver.py b/pcmdi_metrics/sea_ice/sea_ice_driver.py index 4e4df4887..f2c619a18 100644 --- a/pcmdi_metrics/sea_ice/sea_ice_driver.py +++ b/pcmdi_metrics/sea_ice/sea_ice_driver.py @@ -4,7 +4,6 @@ import logging import os -import matplotlib.pyplot as plt import numpy as np import xarray import xcdat as xc @@ -124,6 +123,22 @@ "np": means["np"], "na": means["na"], } + if to_nc: + # Generate netcdf files of climatologies + print("Generating climatology for netcdf") + nc_dir = os.path.join(metrics_output_path, "netcdf") + if not os.path.exists(nc_dir): + os.mkdir(nc_dir) + nc_climo = lib.get_clim(obs, obs_var, ds=None) + print("Writing climatology netcdf") + fname_nh = ( + "sic_clim_" + + "_".join([reference_data_set, "nh", str(osyear), str(oeyear)]) + + ".nc" + ) + fname_nh = os.path.join(nc_dir, fname_nh) + nc_climo.to_netcdf(fname_nh) + del nc_climo obs.close() antarctic_clims = {} @@ -166,6 +181,22 @@ "sp": means["sp"], "sa": means["sa"], } + if to_nc: + # Generate netcdf files of climatologies + print("Generating climatology for netcdf") + nc_dir = os.path.join(metrics_output_path, "netcdf") + if not os.path.exists(nc_dir): + os.mkdir(nc_dir) + nc_climo = lib.get_clim(obs, obs_var, ds=None) + print("Writing climatology netcdf") + fname_sh = ( + "sic_clim_" + + "_".join([reference_data_set, "sh", str(osyear), str(oeyear)]) + + ".nc" + ) + fname_sh = os.path.join(nc_dir, fname_sh) + nc_climo.to_netcdf(fname_sh) + del nc_climo obs.close() obs_clims = {reference_data_set: {}} @@ -383,27 +414,7 @@ ) fname = os.path.join(nc_dir, fname) nc_climo.to_netcdf(fname) - - # Filling in NaNs with zeros for prettier plot - nc_climo = nc_climo.fillna(0) - nc_climo[var] = nc_climo[var].where(mask < 1) - fig_dir = os.path.join(metrics_output_path, "plot") - if not os.path.exists(fig_dir): - os.mkdir(fig_dir) - for mo in range(0, 12): - tmp_title = "_".join([model, run, str(mo + 1)]) - fig.create_arctic_map( - nc_climo.isel({"time": mo}), - fig_dir, - varname="siconc", - title=tmp_title, - ) - fig.create_antarctic_map( - nc_climo.isel({"time": mo}), - fig_dir, - varname="siconc", - title=tmp_title, - ) + del nc_climo # Get regions clims, means = lib.process_by_region(ds, var, area[area_var].data, pole) @@ -557,137 +568,88 @@ # Make figure # ---------------- if plot: - sector_list = [ - "Central Arctic Sector", - "North Atlantic Sector", - "North Pacific Sector", - "Indian Ocean Sector", - "South Atlantic Sector", - "South Pacific Sector", - ] - sector_short = ["ca", "na", "np", "io", "sa", "sp"] - fig7, ax7 = plt.subplots(6, 1, figsize=(5, 9)) - mlabels = model_list - ind = np.arange(len(mlabels)) # the x locations for the groups - width = 0.3 - n = len(ind) - for inds, sector in enumerate(sector_list): - # Assemble data - mse_clim = [] - mse_ext = [] - clim_range = [] - ext_range = [] - clim_err_x = [] - clim_err_y = [] - ext_err_y = [] - rgn = sector_short[inds] - for nmod, model in enumerate(model_list): - mse_clim.append( - float( - metrics["RESULTS"][model][rgn]["model_mean"][ - reference_data_set - ]["monthly_clim"]["mse"] - ) - ) - mse_ext.append( - float( - metrics["RESULTS"][model][rgn]["model_mean"][ - reference_data_set - ]["total_extent"]["mse"] - ) - ) - # Get spread, only if there are multiple realizations - if len(metrics["RESULTS"][model][rgn].keys()) > 2: - for r in metrics["RESULTS"][model][rgn]: - if r != "model_mean": - clim_err_x.append(ind[nmod]) - clim_err_y.append( - float( - metrics["RESULTS"][model][rgn][r][ - reference_data_set - ]["monthly_clim"]["mse"] - ) - ) - ext_err_y.append( - float( - metrics["RESULTS"][model][rgn][r][ - reference_data_set - ]["total_extent"]["mse"] - ) - ) + fig_dir = os.path.join(metrics_output_path, "plot") + if not os.path.exists(fig_dir): + os.mkdir(fig_dir) + + # Make metrics bar chart for all models + print("Creating metrics bar chart.") + meta = fig.metrics_bar_chart( + model_list, metrics, reference_data_set, fig_dir, meta + ) - # plot data - if len(model_list) < 4: - mark_size = 9 - elif len(model_list) < 12: - mark_size = 3 - else: - mark_size = 1 - ax7[inds].bar( - ind - width / 2.0, mse_clim, width, color="b", label="Ann. Cycle" - ) - ax7[inds].bar(ind, mse_ext, width, color="r", label="Ann. Mean") - if len(clim_err_x) > 0: - ax7[inds].scatter( - [x - width / 2.0 for x in clim_err_x], - clim_err_y, - marker="D", - s=mark_size, - color="k", + # Make annual cycle line plots for all models + print("Creating annual cycle figures.") + meta = fig.annual_cycle_plots(df, fig_dir, reference_data_set, meta) + + # Make maps + print("Creating maps.") + if os.path.exists(fname_nh) and os.path.exists(fname_sh): + obs_nh = xc.open_dataset(fname_nh) + obs_sh = xc.open_dataset(fname_sh) + else: + obs_nh = None + obs_sh = None + nc_dir = os.path.join(metrics_output_path, "netcdf/*") + count = 0 + for file in glob.glob(nc_dir): + model = os.path.basename(file).split("_")[2] + run = os.path.basename(file).split("_")[3] + nc_climo = xc.open_dataset(file) + if model == reference_data_set: + continue + try: + tmp_title = "_".join([model, run]) + meta = fig.create_arctic_map( + nc_climo, + obs_nh, + fig_dir, + meta, + varname=var, + title=tmp_title, ) - ax7[inds].scatter( - clim_err_x, ext_err_y, marker="D", s=mark_size, color="k" + meta = fig.create_antarctic_map( + nc_climo, + obs_sh, + fig_dir, + meta, + varname=var, + title=tmp_title, ) - # xticks - if inds == len(sector_list) - 1: - ax7[inds].set_xticks(ind + width / 2.0, mlabels, rotation=90, size=7) - else: - ax7[inds].set_xticks(ind + width / 2.0, labels="") - # yticks - if len(clim_err_y) > 0: - datamax = np.max( - np.concatenate((np.array(clim_err_y), np.array(ext_err_y))) + except Exception as e: + print("Error making figures for model", model, "realization", run) + print(" ", e) + + # Model Mean maps + for model in model_list: + count = 0 + nc_dir = os.path.join(metrics_output_path, "netcdf/*" + model + "_*") + for file in glob.glob(nc_dir): + nc_climo = xc.open_dataset(file) + if count == 0: + nc_climo_mean = nc_climo.copy(deep=True) + else: + nc_climo_mean[var] = nc_climo_mean[var] + nc_climo[var] + nc_climo_mean[var] = nc_climo_mean[var] / (count + 1) + try: + tmp_title = "_".join([model, "model_mean"]) + meta = fig.create_arctic_map( + nc_climo_mean, + fig_dir, + meta, + varname="siconc", + title=tmp_title, ) - else: - datamax = np.max( - np.concatenate((np.array(mse_clim), np.array(mse_ext))) + meta = fig.create_antarctic_map( + nc_climo_mean, + fig_dir, + meta, + varname="siconc", + title=tmp_title, ) - ymax = (datamax) * 1.3 - ax7[inds].set_ylim(0.0, ymax) - if ymax < 0.1: - ticks = np.linspace(0, 0.1, 6) - labels = [str(round(x, 3)) for x in ticks] - elif ymax < 1: - ticks = np.linspace(0, 1, 5) - labels = [str(round(x, 1)) for x in ticks] - elif ymax < 4: - ticks = np.linspace(0, round(ymax), num=round(ymax / 2) * 2 + 1) - labels = [str(round(x, 1)) for x in ticks] - elif ymax > 10: - ticks = range(0, round(ymax), 5) - labels = [str(round(x, 0)) for x in ticks] - else: - ticks = range(0, round(ymax)) - labels = [str(round(x, 0)) for x in ticks] - - ax7[inds].set_yticks(ticks, labels, fontsize=8) - # labels etc - ax7[inds].set_ylabel("10${^1}{^2}$km${^4}$", size=8) - ax7[inds].grid(True, linestyle=":") - ax7[inds].annotate( - sector, - (0.35, 0.8), - xycoords="axes fraction", - size=9, - ) - # Add legend, save figure - ax7[0].legend(loc="upper right", fontsize=6) - plt.suptitle("Mean Square Error relative to " + reference_data_set, y=0.91) - figfile = os.path.join(metrics_output_path, "MSE_bar_chart.png") - plt.savefig(figfile) - meta.update_plots( - "bar_chart", figfile, "regional_bar_chart", "Bar chart of regional MSE" - ) + except Exception as e: + print("Error making figures for model", model, "mean.") + print(" ", e) # ----------------- # Update and write @@ -701,6 +663,7 @@ # Skip provenance if there's an issue print("Error: Could not get provenance from metrics json for output.json.") + print("Writing metadata file.") meta.update_provenance("modeldata", test_data_path) if reference_data_path_nh is not None: meta.update_provenance("obsdata_nh", reference_data_path_nh) From 521b654bcbe27991d4ba540ab1ceb56f46fac424 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 18 Jun 2024 15:18:19 -0700 Subject: [PATCH 10/14] update figures --- pcmdi_metrics/sea_ice/lib/sea_ice_figures.py | 537 +++++++++++++++---- 1 file changed, 440 insertions(+), 97 deletions(-) diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py index 9238580f9..e899cf68b 100644 --- a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py @@ -20,16 +20,15 @@ def replace_fill_zero(da, val): return da_new -def create_arctic_map(ds, metrics_output_path, varname="siconc", title=None): - print("Creating Arctic map") +def create_arctic_map(ds, obs, metrics_output_path, meta, varname="siconc", title=None): # Load and process data xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) # Some models have NaN values in coordinates # that can't be plotted by pcolormesh - ds[xvar] = replace_nan_zero(ds[xvar]) - ds[yvar] = replace_nan_zero(ds[yvar]) + # ds[xvar] = replace_nan_zero(ds[xvar]) + # ds[yvar] = replace_nan_zero(ds[yvar]) # Set up regions region_NA = np.array([[-120, 45], [-120, 80], [90, 80], [90, 45]]) @@ -45,73 +44,146 @@ def create_arctic_map(ds, metrics_output_path, varname="siconc", title=None): "", [[0, 85 / 255, 182 / 255], "white"] ) proj = ccrs.NorthPolarStereo() - ax = plt.subplot(111, projection=proj) - ax.set_global() - # TODO: get xvar and yvar - ds[varname].plot.pcolormesh( - ax=ax, - x=xvar, - y=yvar, - transform=ccrs.PlateCarree(), - cmap=cmap, - cbar_kwargs={"label": "ice fraction"}, - ) - arctic_regions.plot_regions( - ax=ax, - add_label=False, - label="abbrev", - line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, - ) - ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) - ax.coastlines(color=[0.3, 0.3, 0.3]) - plt.annotate( - "North Atlantic", - (0.5, 0.2), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="white", - ) - plt.annotate( - "North Pacific", - (0.65, 0.88), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="white", - ) - plt.annotate( - "Central\nArctic ", - (0.56, 0.56), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - ) - ax.set_facecolor([0.55, 0.55, 0.6]) + f1, axs = plt.subplots(12, 2, figsize=(10, 60), subplot_kw={"projection": proj}) + pos1 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] + pos2 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + months = [ + "January", + "February", + "March", + "April", + "May", + "June", + "July", + "August", + "September", + "October", + "November", + "December", + ] + for mo in range(0, 12): + ax = axs[pos1[mo], pos2[mo]] + ax.set_global() + + ds[varname].isel({"time": mo}).plot.pcolormesh( + ax=ax, + x=xvar, + y=yvar, + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) + arctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + """plt.annotate( + "North Atlantic", + (0.5, 0.2), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="white", + ) + plt.annotate( + "North Pacific", + (0.65, 0.88), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="white", + ) + plt.annotate( + "Central\nArctic ", + (0.56, 0.56), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + )""" + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title(months[mo]) + for mo in range(0, 12): + ax = axs[pos1[mo], pos2[mo]] + ax.set_global() + + obs[varname].isel({"time": mo}).plot.pcolormesh( + ax=ax, + x=xvar, + y=yvar, + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) + arctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + """plt.annotate( + "North Atlantic", + (0.5, 0.2), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="white", + ) + plt.annotate( + "North Pacific", + (0.65, 0.88), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="white", + ) + plt.annotate( + "Central\nArctic ", + (0.56, 0.56), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + )""" + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title("Obs", months[mo]) if title is not None: - plt.title(title) + plt.suptitle(title.replace("_", " "), fontsize=20) fig_path = os.path.join( metrics_output_path, title.replace(" ", "_") + "_arctic_regions.png" ) else: fig_path = os.path.join(metrics_output_path, "arctic_regions.png") + plt.tight_layout(rect=[0, 0.03, 1, 0.97]) plt.savefig(fig_path) plt.close() + meta.update_plots( + "Arctic" + title.replace(" ", "_"), + fig_path, + "Arctic_map", + "Maps of monthly Antarctic ice fraction", + ) + return meta # ---------- # Antarctic # ---------- -def create_antarctic_map(ds, metrics_output_path, varname="siconc", title=None): - print("Creating Antarctic map") +def create_antarctic_map( + ds, obs, metrics_output_path, meta, varname="siconc", title=None +): # Load and process data xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) # Some models have NaN values in coordinates # that can't be plotted by pcolormesh - ds[xvar] = replace_nan_zero(ds[xvar]) - ds[yvar] = replace_nan_zero(ds[yvar]) + # ds[xvar] = replace_nan_zero(ds[xvar]) + # ds[yvar] = replace_nan_zero(ds[yvar]) # Set up regions region_IO = np.array([[20, -90], [90, -90], [90, -55], [20, -55]]) @@ -120,7 +192,7 @@ def create_antarctic_map(ds, metrics_output_path, varname="siconc", title=None): names = ["Indian Ocean", "South Atlantic", "South Pacific"] abbrevs = ["IO", "SA", "SP"] - arctic_regions = regionmask.Regions( + antarctic_regions = regionmask.Regions( [region_IO, region_SA, region_SP], names=names, abbrevs=abbrevs, @@ -128,59 +200,330 @@ def create_antarctic_map(ds, metrics_output_path, varname="siconc", title=None): ) # Do plotting + cmap = colors.LinearSegmentedColormap.from_list( "", [[0, 85 / 255, 182 / 255], "white"] ) proj = ccrs.SouthPolarStereo() - ax = plt.subplot(111, projection=proj) - ax.set_global() - ds[varname].plot.pcolormesh( - ax=ax, - x=xvar, - y=yvar, - transform=ccrs.PlateCarree(), - cmap=cmap, - cbar_kwargs={"label": "ice fraction"}, - ) - arctic_regions.plot_regions( - ax=ax, - add_label=False, - label="abbrev", - line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, - ) - ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) - ax.coastlines(color=[0.3, 0.3, 0.3]) - plt.annotate( - "South Pacific", - (0.50, 0.2), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="black", - ) - plt.annotate( - "Indian\nOcean", - (0.89, 0.66), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="black", - ) - plt.annotate( - "South Atlantic", - (0.54, 0.82), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="black", - ) - ax.set_facecolor([0.55, 0.55, 0.6]) + f1, axs = plt.subplots(12, 2, figsize=(10, 60), subplot_kw={"projection": proj}) + pos1 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] + pos2 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + months = [ + "January", + "February", + "March", + "April", + "May", + "June", + "July", + "August", + "September", + "October", + "November", + "December", + ] + for mo in range(0, 12): + ax = axs[pos1[mo], pos2[mo]] + ax.set_global() + ds[varname].isel({"time": mo}).plot.pcolormesh( + ax=ax, + x=xvar, + y=yvar, + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) + antarctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + """ax.annotate( + "South Pacific", + (0.50, 0.2), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="black", + ) + ax.annotate( + "Indian\nOcean", + (0.89, 0.66), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="black", + ) + ax.annotate( + "South Atlantic", + (0.54, 0.82), + xycoords="axes fraction", + horizontalalignment="right", + verticalalignment="bottom", + color="black", + )""" + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title(months[mo]) if title is not None: - plt.title(title) + plt.suptitle(title.replace("_", " "), fontsize=20) fig_path = os.path.join( metrics_output_path, title.replace(" ", "_") + "_antarctic_regions.png" ) else: fig_path = os.path.join(metrics_output_path, "antarctic_regions.png") + plt.tight_layout(rect=[0, 0.03, 1, 0.97]) plt.savefig(fig_path) plt.close() + meta.update_plots( + "Antarctic" + title.replace(" ", "_"), + fig_path, + "Antarctic_map", + "Maps of monthly Antarctic ice fraction", + ) + return meta + + +def annual_cycle_plots(df, fig_dir, reference_data_set, meta): + # Annual cycle line figure + # Model mean climatology + labels = { + "antarctic": "Antarctic", + "arctic": "Arctic", + "ca": "Central Arctic", + "na": "North Atlantic", + "np": "North Pacific", + "io": "Indian Ocean", + "sa": "South Atlantic", + "sp": "South Pacific", + } + pos = { + "antarctic": (0, 1), + "arctic": (0, 0), + "ca": (1, 0), + "na": (2, 0), + "np": (3, 0), + "io": (1, 1), + "sa": (2, 1), + "sp": (3, 1), + } + wts = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) / 365 + for mod in df: + try: + f1, axs = plt.subplots(4, 2, figsize=(8, 8)) + if mod != "Reference": + for rgn in ["arctic", "antarctic", "ca", "sa", "np", "sp", "na", "io"]: + # Get realization spread + climr = [] + for real in df[mod][rgn]: + if real != "model_mean": + tmp = df[mod][rgn][real]["monthly_climatology"] + tmp = [float(x) for x in tmp] + climr = climr + tmp + climr = np.array(climr) + climr = np.reshape(climr, (int(len(climr) / 12), 12)) + mins = np.min(climr, axis=0) + maxs = np.max(climr, axis=0) + + # Get means + clim = df[mod][rgn]["model_mean"]["monthly_climatology"] + climobs = df["Reference"][rgn][reference_data_set][ + "monthly_climatology" + ] + clim = np.array([float(x) for x in clim]) + climobs = np.array([float(x) for x in climobs]) + time = np.array([x for x in range(1, 13)]) + climmean = np.ones(12) * np.average(clim, weights=wts) + climmeanobs = np.ones(12) * np.average(climobs, weights=wts) + + # Make figure + ind1 = pos[rgn][0] + ind2 = pos[rgn][1] + axs[ind1, ind2].plot( + time, + climmean, + color="darkorange", + linewidth=1.5, + linestyle="dashed", + zorder=1, + ) + axs[ind1, ind2].plot( + time, + climmeanobs, + color=[0.2, 0.2, 0.2], + linewidth=1.5, + linestyle="dashed", + zorder=1, + ) + axs[ind1, ind2].fill_between( + time, mins, maxs, color="darkorange", alpha=0.3, zorder=1 + ) + axs[ind1, ind2].plot(time, clim, color="darkorange", label=mod) + axs[ind1, ind2].plot( + time, climobs, color=[0.2, 0.2, 0.2], label=reference_data_set + ) + axs[ind1, ind2].set_title(labels[rgn], fontsize=10) + axs[ind1, ind2].xaxis.set_tick_params(which="minor", bottom=False) + axs[ind1, ind2].set_xlim([1, 12]) + axs[ind1, ind2].set_xticks([]) + axs[ind1, ind2].minorticks_on() + if ind2 == 0: + axs[ind1, ind2].set_ylabel("$10^6$ $km^2$") + axs[3, 0].set_xticks(time) + axs[3, 1].set_xticks(time) + axs[3, 0].set_xlabel("month") + axs[3, 1].set_xlabel("month") + axs[0, 0].legend(loc="upper right", fontsize=9) + plt.suptitle(mod) + figfile = os.path.join(fig_dir, mod + "_annual_cycle.png") + plt.savefig(figfile) + plt.close() + meta.update_plots( + "ann_cycle_" + mod, + figfile, + "annual_cycle_plot", + "Plot of model mean regional annual cycles in sea ice area", + ) + except Exception as e: + print("Could not create annual cycle plots for model", mod) + print(e) + return meta + + +def metrics_bar_chart(model_list, metrics, reference_data_set, fig_dir, meta): + try: + sector_list = [ + "Central Arctic Sector", + "North Atlantic Sector", + "North Pacific Sector", + "Indian Ocean Sector", + "South Atlantic Sector", + "South Pacific Sector", + ] + sector_short = ["ca", "na", "np", "io", "sa", "sp"] + fig7, ax7 = plt.subplots(6, 1, figsize=(5, 9)) + mlabels = model_list + ind = np.arange(len(mlabels)) # the x locations for the groups + width = 0.3 + for inds, sector in enumerate(sector_list): + # Assemble data + mse_clim = [] + mse_ext = [] + clim_err_x = [] + clim_err_y = [] + ext_err_y = [] + rgn = sector_short[inds] + for nmod, model in enumerate(model_list): + mse_clim.append( + float( + metrics["RESULTS"][model][rgn]["model_mean"][ + reference_data_set + ]["monthly_clim"]["mse"] + ) + ) + mse_ext.append( + float( + metrics["RESULTS"][model][rgn]["model_mean"][ + reference_data_set + ]["total_extent"]["mse"] + ) + ) + # Get spread, only if there are multiple realizations + if len(metrics["RESULTS"][model][rgn].keys()) > 2: + for r in metrics["RESULTS"][model][rgn]: + if r != "model_mean": + clim_err_x.append(ind[nmod]) + clim_err_y.append( + float( + metrics["RESULTS"][model][rgn][r][ + reference_data_set + ]["monthly_clim"]["mse"] + ) + ) + ext_err_y.append( + float( + metrics["RESULTS"][model][rgn][r][ + reference_data_set + ]["total_extent"]["mse"] + ) + ) + + # plot data + if len(model_list) < 4: + mark_size = 9 + elif len(model_list) < 12: + mark_size = 3 + else: + mark_size = 1 + ax7[inds].bar( + ind - width / 2.0, mse_clim, width, color="b", label="Ann. Cycle" + ) + ax7[inds].bar(ind, mse_ext, width, color="r", label="Ann. Mean") + if len(clim_err_x) > 0: + ax7[inds].scatter( + [x - width / 2.0 for x in clim_err_x], + clim_err_y, + marker="D", + s=mark_size, + color="k", + ) + ax7[inds].scatter( + clim_err_x, ext_err_y, marker="D", s=mark_size, color="k" + ) + # xticks + if inds == len(sector_list) - 1: + ax7[inds].set_xticks(ind + width / 2.0, mlabels, rotation=90, size=7) + else: + ax7[inds].set_xticks(ind + width / 2.0, labels="") + # yticks + if len(clim_err_y) > 0: + datamax = np.max( + np.concatenate((np.array(clim_err_y), np.array(ext_err_y))) + ) + else: + datamax = np.max( + np.concatenate((np.array(mse_clim), np.array(mse_ext))) + ) + ymax = (datamax) * 1.3 + ax7[inds].set_ylim(0.0, ymax) + if ymax < 0.1: + ticks = np.linspace(0, 0.1, 6) + labels = [str(round(x, 3)) for x in ticks] + elif ymax < 1: + ticks = np.linspace(0, 1, 5) + labels = [str(round(x, 1)) for x in ticks] + elif ymax < 4: + ticks = np.linspace(0, round(ymax), num=round(ymax / 2) * 2 + 1) + labels = [str(round(x, 1)) for x in ticks] + elif ymax > 10: + ticks = range(0, round(ymax), 5) + labels = [str(round(x, 0)) for x in ticks] + else: + ticks = range(0, round(ymax)) + labels = [str(round(x, 0)) for x in ticks] + + ax7[inds].set_yticks(ticks, labels, fontsize=8) + # labels etc + ax7[inds].set_ylabel("10${^1}{^2}$km${^4}$", size=8) + ax7[inds].grid(True, linestyle=":") + ax7[inds].annotate( + sector, + (0.35, 0.8), + xycoords="axes fraction", + size=9, + ) + # Add legend, save figure + ax7[0].legend(loc="upper right", fontsize=6) + plt.suptitle("Mean Square Error relative to " + reference_data_set, y=0.91) + figfile = os.path.join(fig_dir, "MSE_bar_chart.png") + plt.savefig(figfile) + plt.close() + meta.update_plots( + "bar_chart", figfile, "regional_bar_chart", "Bar chart of regional MSE" + ) + except Exception as e: + print("Could not create metrics plot.") + print(e) + return meta From 0d8002eab554cb2d5f044b8ae184bedd602a5117 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 21 Jun 2024 15:12:42 -0700 Subject: [PATCH 11/14] update figures --- pcmdi_metrics/sea_ice/sea_ice_driver.py | 55 +++++++++++++------------ 1 file changed, 29 insertions(+), 26 deletions(-) diff --git a/pcmdi_metrics/sea_ice/sea_ice_driver.py b/pcmdi_metrics/sea_ice/sea_ice_driver.py index f2c619a18..ddb727152 100644 --- a/pcmdi_metrics/sea_ice/sea_ice_driver.py +++ b/pcmdi_metrics/sea_ice/sea_ice_driver.py @@ -125,7 +125,6 @@ } if to_nc: # Generate netcdf files of climatologies - print("Generating climatology for netcdf") nc_dir = os.path.join(metrics_output_path, "netcdf") if not os.path.exists(nc_dir): os.mkdir(nc_dir) @@ -183,7 +182,6 @@ } if to_nc: # Generate netcdf files of climatologies - print("Generating climatology for netcdf") nc_dir = os.path.join(metrics_output_path, "netcdf") if not os.path.exists(nc_dir): os.mkdir(nc_dir) @@ -401,7 +399,6 @@ if to_nc: # Generate netcdf files of climatologies - print("Generating climatology for netcdf") nc_dir = os.path.join(metrics_output_path, "netcdf") if not os.path.exists(nc_dir): os.mkdir(nc_dir) @@ -603,18 +600,20 @@ meta = fig.create_arctic_map( nc_climo, obs_nh, + var, + obs_var, fig_dir, meta, - varname=var, - title=tmp_title, + tmp_title, ) meta = fig.create_antarctic_map( nc_climo, obs_sh, + var, + obs_var, fig_dir, meta, - varname=var, - title=tmp_title, + tmp_title, ) except Exception as e: print("Error making figures for model", model, "realization", run) @@ -631,25 +630,29 @@ else: nc_climo_mean[var] = nc_climo_mean[var] + nc_climo[var] nc_climo_mean[var] = nc_climo_mean[var] / (count + 1) - try: - tmp_title = "_".join([model, "model_mean"]) - meta = fig.create_arctic_map( - nc_climo_mean, - fig_dir, - meta, - varname="siconc", - title=tmp_title, - ) - meta = fig.create_antarctic_map( - nc_climo_mean, - fig_dir, - meta, - varname="siconc", - title=tmp_title, - ) - except Exception as e: - print("Error making figures for model", model, "mean.") - print(" ", e) + # try: + tmp_title = "_".join([model, "model_mean"]) + meta = fig.create_arctic_map( + nc_climo_mean, + obs_nh, + var, + obs_var, + fig_dir, + meta, + tmp_title, + ) + meta = fig.create_antarctic_map( + nc_climo_mean, + obs_sh, + var, + obs_var, + fig_dir, + meta, + tmp_title, + ) + # except Exception as e: + # print("Error making figures for model",model,"mean.") + # print(" ",e) # ----------------- # Update and write From 02e35c8d299ea4faba02f63603f7e06053d844dd Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 28 Jun 2024 13:47:52 -0700 Subject: [PATCH 12/14] add summary --- pcmdi_metrics/sea_ice/lib/sea_ice_figures.py | 515 +++++++++++++------ 1 file changed, 350 insertions(+), 165 deletions(-) diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py index e899cf68b..17342b3a2 100644 --- a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py @@ -20,8 +20,7 @@ def replace_fill_zero(da, val): return da_new -def create_arctic_map(ds, obs, metrics_output_path, meta, varname="siconc", title=None): - # Load and process data +def create_summary_maps_arctic(ds, var_ice, metrics_output_path, meta, model): xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) @@ -39,39 +38,29 @@ def create_arctic_map(ds, obs, metrics_output_path, meta, varname="siconc", titl [region_NA, region_NP], names=names, abbrevs=abbrevs, name="arctic" ) - # Do plotting cmap = colors.LinearSegmentedColormap.from_list( "", [[0, 85 / 255, 182 / 255], "white"] ) - proj = ccrs.NorthPolarStereo() - f1, axs = plt.subplots(12, 2, figsize=(10, 60), subplot_kw={"projection": proj}) - pos1 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] - pos2 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] - months = [ - "January", - "February", - "March", - "April", - "May", - "June", - "July", - "August", - "September", - "October", - "November", - "December", - ] - for mo in range(0, 12): - ax = axs[pos1[mo], pos2[mo]] - ax.set_global() - ds[varname].isel({"time": mo}).plot.pcolormesh( + # Do plotting + try: + proj = ccrs.NorthPolarStereo() + f1, axs = plt.subplots( + 1, 2, figsize=(10, 5), subplot_kw={"projection": proj}, layout="compressed" + ) + + # Model arctic Feb + ax = axs[0] + ax.set_global() # to use cartopy polar proj + + ds[var_ice].isel({"time": 1}).plot( ax=ax, x=xvar, y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], transform=ccrs.PlateCarree(), cmap=cmap, - cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + add_colorbar=False, ) arctic_regions.plot_regions( ax=ax, @@ -81,42 +70,26 @@ def create_arctic_map(ds, obs, metrics_output_path, meta, varname="siconc", titl ) ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) ax.coastlines(color=[0.3, 0.3, 0.3]) - """plt.annotate( - "North Atlantic", - (0.5, 0.2), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="white", - ) - plt.annotate( - "North Pacific", - (0.65, 0.88), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="white", - ) - plt.annotate( - "Central\nArctic ", - (0.56, 0.56), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - )""" ax.set_facecolor([0.55, 0.55, 0.6]) - ax.set_title(months[mo]) - for mo in range(0, 12): - ax = axs[pos1[mo], pos2[mo]] - ax.set_global() + ax.set_title("Feb\n" + model.replace("_", " "), fontsize=12) - obs[varname].isel({"time": mo}).plot.pcolormesh( - ax=ax, - x=xvar, - y=yvar, - transform=ccrs.PlateCarree(), - cmap=cmap, - cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + # Model Arctic Sept + ax = axs[1] + ax.set_global() # to use cartopy polar proj + + fds = ( + ds[var_ice] + .isel({"time": 8}) + .plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + add_colorbar=False, + # cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) ) arctic_regions.plot_regions( ax=ax, @@ -126,57 +99,32 @@ def create_arctic_map(ds, obs, metrics_output_path, meta, varname="siconc", titl ) ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) ax.coastlines(color=[0.3, 0.3, 0.3]) - """plt.annotate( - "North Atlantic", - (0.5, 0.2), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="white", - ) - plt.annotate( - "North Pacific", - (0.65, 0.88), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="white", - ) - plt.annotate( - "Central\nArctic ", - (0.56, 0.56), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - )""" ax.set_facecolor([0.55, 0.55, 0.6]) - ax.set_title("Obs", months[mo]) - if title is not None: - plt.suptitle(title.replace("_", " "), fontsize=20) + ax.set_title("Sep\n" + model.replace("_", " "), fontsize=12) + + # plt.suptitle("Arctic", fontsize=30) fig_path = os.path.join( - metrics_output_path, title.replace(" ", "_") + "_arctic_regions.png" + metrics_output_path, model.replace(" ", "_") + "_Feb_Sep_NH.png" ) - else: - fig_path = os.path.join(metrics_output_path, "arctic_regions.png") - plt.tight_layout(rect=[0, 0.03, 1, 0.97]) - plt.savefig(fig_path) - plt.close() - meta.update_plots( - "Arctic" + title.replace(" ", "_"), - fig_path, - "Arctic_map", - "Maps of monthly Antarctic ice fraction", - ) + # plt.tight_layout(rect=[0, 0.03, 1, 0.97]) + plt.colorbar(fds, label="ice fraction", ax=axs) + plt.savefig(fig_path) + plt.close() + meta.update_plots( + "Summary_NH_" + model.replace(" ", "_"), + fig_path, + "Feb_Sep_maps", + "Summary map of Feb and Sep ice areas for Northern hemispheres", + ) + except Exception as e: + print("Could not create summary maps.") + print(e) + if plt.get_fignums(): + plt.close() return meta -# ---------- -# Antarctic -# ---------- -def create_antarctic_map( - ds, obs, metrics_output_path, meta, varname="siconc", title=None -): - # Load and process data +def create_summary_maps_antarctic(ds, var_ice, metrics_output_path, meta, model): xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) @@ -199,39 +147,32 @@ def create_antarctic_map( name="antarctic", ) - # Do plotting - cmap = colors.LinearSegmentedColormap.from_list( "", [[0, 85 / 255, 182 / 255], "white"] ) - proj = ccrs.SouthPolarStereo() - f1, axs = plt.subplots(12, 2, figsize=(10, 60), subplot_kw={"projection": proj}) - pos1 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] - pos2 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] - months = [ - "January", - "February", - "March", - "April", - "May", - "June", - "July", - "August", - "September", - "October", - "November", - "December", - ] - for mo in range(0, 12): - ax = axs[pos1[mo], pos2[mo]] + + # Do plotting + try: + # proj = ccrs.NorthPolarStereo() + f1, axs = plt.subplots( + 1, + 2, + figsize=(10, 5), + subplot_kw={"projection": ccrs.SouthPolarStereo()}, + layout="compressed", + ) + + # Model Antarctic September + ax = axs[0] ax.set_global() - ds[varname].isel({"time": mo}).plot.pcolormesh( + ds[var_ice].isel({"time": 8}).plot( ax=ax, x=xvar, y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], transform=ccrs.PlateCarree(), cmap=cmap, - cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + add_colorbar=False, ) antarctic_regions.plot_regions( ax=ax, @@ -241,48 +182,288 @@ def create_antarctic_map( ) ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) ax.coastlines(color=[0.3, 0.3, 0.3]) - """ax.annotate( - "South Pacific", - (0.50, 0.2), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="black", + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title("Sep\n" + model.replace("_", " "), fontsize=12) + + # Model Antarctic Feb + ax = axs[1] + ax.set_global() + fds = ( + ds[var_ice] + .isel({"time": 1}) + .plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + add_colorbar=False, + # cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) ) - ax.annotate( - "Indian\nOcean", - (0.89, 0.66), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="black", + antarctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, ) - ax.annotate( - "South Atlantic", - (0.54, 0.82), - xycoords="axes fraction", - horizontalalignment="right", - verticalalignment="bottom", - color="black", - )""" + ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) ax.set_facecolor([0.55, 0.55, 0.6]) - ax.set_title(months[mo]) - if title is not None: - plt.suptitle(title.replace("_", " "), fontsize=20) + ax.set_title("Feb\n" + model.replace("_", " "), fontsize=12) + + # plt.suptitle("Arctic", fontsize=30) fig_path = os.path.join( - metrics_output_path, title.replace(" ", "_") + "_antarctic_regions.png" + metrics_output_path, model.replace(" ", "_") + "_Feb_Sep_SH.png" + ) + # plt.tight_layout(rect=[0, 0.03, 1, 0.97]) + plt.colorbar(fds, label="ice fraction", ax=axs) + plt.savefig(fig_path) + plt.close() + meta.update_plots( + "Summary_SH_" + model.replace(" ", "_"), + fig_path, + "Feb_Sep_maps", + "Summary map of Feb and Sep ice areas for Southern hemispheres", ) - else: - fig_path = os.path.join(metrics_output_path, "antarctic_regions.png") - plt.tight_layout(rect=[0, 0.03, 1, 0.97]) - plt.savefig(fig_path) - plt.close() - meta.update_plots( - "Antarctic" + title.replace(" ", "_"), - fig_path, - "Antarctic_map", - "Maps of monthly Antarctic ice fraction", + except Exception as e: + print("Could not create summary maps.") + print(e) + if plt.get_fignums(): + plt.close() + return meta + + +def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, model): + # Load and process data + xvar = lib.find_lon(ds) + yvar = lib.find_lat(ds) + + # Some models have NaN values in coordinates + # that can't be plotted by pcolormesh + # ds[xvar] = replace_nan_zero(ds[xvar]) + # ds[yvar] = replace_nan_zero(ds[yvar]) + + # Set up regions + region_NA = np.array([[-120, 45], [-120, 80], [90, 80], [90, 45]]) + region_NP = np.array([[90, 45], [90, 65], [240, 65], [240, 45]]) + names = ["North_Atlantic", "North_Pacific"] + abbrevs = ["NA", "NP"] + arctic_regions = regionmask.Regions( + [region_NA, region_NP], names=names, abbrevs=abbrevs, name="arctic" + ) + + # Do plotting + try: + cmap = colors.LinearSegmentedColormap.from_list( + "", [[0, 85 / 255, 182 / 255], "white"] + ) + proj = ccrs.NorthPolarStereo() + f1, axs = plt.subplots(3, 8, figsize=(30, 12), subplot_kw={"projection": proj}) + pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] + pos2 = [1, 1, 3, 3, 3, 5, 5, 5, 7, 7, 7, 1] + months = [ + "January", + "February", + "March", + "April", + "May", + "June", + "July", + "August", + "September", + "October", + "November", + "December", + ] + for mo in range(0, 12): + ax = axs[pos1[mo], pos2[mo]] + ax.set_global() # to use cartopy polar proj + + ds[var_ice].isel({"time": mo}).plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) + arctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title(months[mo] + "\n" + model.replace("_", " "), fontsize=12) + pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] + pos2 = [0, 0, 2, 2, 2, 4, 4, 4, 6, 6, 6, 0] + xvar = lib.find_lon(obs) + yvar = lib.find_lat(obs) + for mo in range(0, 12): + ax = axs[pos1[mo], pos2[mo]] + # ax.set_global() + obs[var_obs].isel({"time": mo}).plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) + arctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title(months[mo] + "\nReference", fontsize=12) + + plt.suptitle("Arctic", fontsize=30) + fig_path = os.path.join( + metrics_output_path, model.replace(" ", "_") + "_arctic_regions.png" + ) + plt.tight_layout(rect=[0, 0.03, 1, 0.97]) + plt.savefig(fig_path) + plt.close() + meta.update_plots( + "Arctic" + model.replace(" ", "_"), + fig_path, + "Arctic_map", + "Maps of monthly Antarctic ice fraction", + ) + except Exception as e: + print("Could not create Arctic maps.") + print(e) + if plt.get_fignums(): + plt.close() + return meta + + +# ---------- +# Antarctic +# ---------- +def create_antarctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, model): + # Load and process data + xvar = lib.find_lon(ds) + yvar = lib.find_lat(ds) + + # Some models have NaN values in coordinates + # that can't be plotted by pcolormesh + # ds[xvar] = replace_nan_zero(ds[xvar]) + # ds[yvar] = replace_nan_zero(ds[yvar]) + + # Set up regions + region_IO = np.array([[20, -90], [90, -90], [90, -55], [20, -55]]) + region_SA = np.array([[20, -90], [-60, -90], [-60, -55], [20, -55]]) + region_SP = np.array([[90, -90], [300, -90], [300, -55], [90, -55]]) + + names = ["Indian Ocean", "South Atlantic", "South Pacific"] + abbrevs = ["IO", "SA", "SP"] + antarctic_regions = regionmask.Regions( + [region_IO, region_SA, region_SP], + names=names, + abbrevs=abbrevs, + name="antarctic", ) + + # Do plotting + try: + cmap = colors.LinearSegmentedColormap.from_list( + "", [[0, 85 / 255, 182 / 255], "white"] + ) + proj = ccrs.SouthPolarStereo() + f1, axs = plt.subplots(3, 8, figsize=(30, 12), subplot_kw={"projection": proj}) + pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] + pos2 = [1, 1, 3, 3, 3, 5, 5, 5, 7, 7, 7, 1] + months = [ + "January", + "February", + "March", + "April", + "May", + "June", + "July", + "August", + "September", + "October", + "November", + "December", + ] + for mo in range(0, 12): + ax = axs[pos1[mo], pos2[mo]] + ax.set_global() + ds[var_ice].isel({"time": mo}).plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) + antarctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title(months[mo] + "\n" + model.replace("_", " "), fontsize=12) + pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] + pos2 = [0, 0, 2, 2, 2, 4, 4, 4, 6, 6, 6, 0] + xvar = lib.find_lon(obs) + yvar = lib.find_lat(obs) + for mo in range(0, 12): + ax = axs[pos1[mo], pos2[mo]] + # ax.set_global() + obs[var_obs].isel({"time": mo}).plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + ) + antarctic_regions.plot_regions( + ax=ax, + add_label=False, + label="abbrev", + line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, + ) + ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) + ax.coastlines(color=[0.3, 0.3, 0.3]) + ax.set_facecolor([0.55, 0.55, 0.6]) + ax.set_title(months[mo] + "\nReference", fontsize=12) + plt.suptitle(model.replace("_", " "), fontsize=20) + fig_path = os.path.join( + metrics_output_path, model.replace(" ", "_") + "_antarctic_regions.png" + ) + plt.tight_layout(rect=[0, 0.03, 1, 0.97]) + plt.savefig(fig_path) + plt.close() + meta.update_plots( + "Antarctic" + model.replace(" ", "_"), + fig_path, + "Antarctic_map", + "Maps of monthly Antarctic ice fraction", + ) + except Exception as e: + print("Could not create Antarctic maps.") + print(e) + if plt.get_fignums(): + plt.close() return meta @@ -389,6 +570,8 @@ def annual_cycle_plots(df, fig_dir, reference_data_set, meta): except Exception as e: print("Could not create annual cycle plots for model", mod) print(e) + if plt.get_fignums(): + plt.close() return meta @@ -526,4 +709,6 @@ def metrics_bar_chart(model_list, metrics, reference_data_set, fig_dir, meta): except Exception as e: print("Could not create metrics plot.") print(e) + if plt.get_fignums(): + plt.close() return meta From 5e9212f09f54430a449ffc443c9d7dab211b4d57 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 5 Jul 2024 16:18:05 -0700 Subject: [PATCH 13/14] update figures --- pcmdi_metrics/sea_ice/lib/sea_ice_figures.py | 171 ++++++++++++++----- 1 file changed, 128 insertions(+), 43 deletions(-) diff --git a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py index 17342b3a2..f8485783e 100644 --- a/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py +++ b/pcmdi_metrics/sea_ice/lib/sea_ice_figures.py @@ -235,7 +235,9 @@ def create_summary_maps_antarctic(ds, var_ice, metrics_output_path, meta, model) return meta -def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, model): +def create_arctic_map( + ds, obs, var_ice, var_obs, metrics_output_path, meta, model, title +): # Load and process data xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) @@ -260,9 +262,11 @@ def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, mode "", [[0, 85 / 255, 182 / 255], "white"] ) proj = ccrs.NorthPolarStereo() - f1, axs = plt.subplots(3, 8, figsize=(30, 12), subplot_kw={"projection": proj}) - pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] - pos2 = [1, 1, 3, 3, 3, 5, 5, 5, 7, 7, 7, 1] + f1, axs = plt.subplots( + 4, 6, figsize=(21, 13), subplot_kw={"projection": proj}, layout="compressed" + ) + pos1 = [1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3] # row + pos2 = [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5] # col months = [ "January", "February", @@ -279,16 +283,18 @@ def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, mode ] for mo in range(0, 12): ax = axs[pos1[mo], pos2[mo]] - ax.set_global() # to use cartopy polar proj - - ds[var_ice].isel({"time": mo}).plot( - ax=ax, - x=xvar, - y=yvar, - levels=[0.15, 0.4, 0.6, 0.8, 1], - transform=ccrs.PlateCarree(), - cmap=cmap, - cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + fds = ( + ds[var_ice] + .isel({"time": mo}) + .plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + add_colorbar=False, + ) ) arctic_regions.plot_regions( ax=ax, @@ -296,17 +302,27 @@ def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, mode label="abbrev", line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, ) + if pos2[mo] == 0: + ax.text( + 0.03, + 0.04, + model.replace("_", " "), + horizontalalignment="left", + verticalalignment="center", + transform=ax.transAxes, + backgroundcolor="white", + fontsize=14, + ) ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) ax.coastlines(color=[0.3, 0.3, 0.3]) ax.set_facecolor([0.55, 0.55, 0.6]) - ax.set_title(months[mo] + "\n" + model.replace("_", " "), fontsize=12) - pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] - pos2 = [0, 0, 2, 2, 2, 4, 4, 4, 6, 6, 6, 0] + ax.set_title("") + pos1 = [0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2] # row + pos2 = [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5] # col xvar = lib.find_lon(obs) yvar = lib.find_lat(obs) for mo in range(0, 12): ax = axs[pos1[mo], pos2[mo]] - # ax.set_global() obs[var_obs].isel({"time": mo}).plot( ax=ax, x=xvar, @@ -314,7 +330,7 @@ def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, mode levels=[0.15, 0.4, 0.6, 0.8, 1], transform=ccrs.PlateCarree(), cmap=cmap, - cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + add_colorbar=False, ) arctic_regions.plot_regions( ax=ax, @@ -322,16 +338,42 @@ def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, mode label="abbrev", line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, ) + if pos2[mo] == 0: + ax.text( + 0.03, + 0.04, + "Reference", + horizontalalignment="left", + verticalalignment="center", + transform=ax.transAxes, + backgroundcolor="white", + fontsize=14, + ) ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) ax.coastlines(color=[0.3, 0.3, 0.3]) ax.set_facecolor([0.55, 0.55, 0.6]) - ax.set_title(months[mo] + "\nReference", fontsize=12) + ax.set_title(months[mo], fontsize=18) - plt.suptitle("Arctic", fontsize=30) + plt.suptitle(title, fontsize=26) fig_path = os.path.join( metrics_output_path, model.replace(" ", "_") + "_arctic_regions.png" ) - plt.tight_layout(rect=[0, 0.03, 1, 0.97]) + cbar = plt.colorbar( + fds, + ax=axs[0:2, :], + pad=0.02, + aspect=25, + ) + cbar.ax.tick_params(labelsize=18) + cbar.set_label(label="ice fraction", size=18) + cbar = plt.colorbar( + fds, + ax=axs[2:, :], + pad=0.02, + aspect=25, + ) + cbar.ax.tick_params(labelsize=18) + cbar.set_label(label="ice fraction", size=18) plt.savefig(fig_path) plt.close() meta.update_plots( @@ -351,7 +393,9 @@ def create_arctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, mode # ---------- # Antarctic # ---------- -def create_antarctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, model): +def create_antarctic_map( + ds, obs, var_ice, var_obs, metrics_output_path, meta, model, title +): # Load and process data xvar = lib.find_lon(ds) yvar = lib.find_lat(ds) @@ -381,9 +425,11 @@ def create_antarctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, m "", [[0, 85 / 255, 182 / 255], "white"] ) proj = ccrs.SouthPolarStereo() - f1, axs = plt.subplots(3, 8, figsize=(30, 12), subplot_kw={"projection": proj}) - pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] - pos2 = [1, 1, 3, 3, 3, 5, 5, 5, 7, 7, 7, 1] + f1, axs = plt.subplots( + 4, 6, figsize=(21, 13), subplot_kw={"projection": proj}, layout="compressed" + ) + pos1 = [1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3] # row + pos2 = [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5] # col months = [ "January", "February", @@ -400,15 +446,18 @@ def create_antarctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, m ] for mo in range(0, 12): ax = axs[pos1[mo], pos2[mo]] - ax.set_global() - ds[var_ice].isel({"time": mo}).plot( - ax=ax, - x=xvar, - y=yvar, - levels=[0.15, 0.4, 0.6, 0.8, 1], - transform=ccrs.PlateCarree(), - cmap=cmap, - cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + fds = ( + ds[var_ice] + .isel({"time": mo}) + .plot( + ax=ax, + x=xvar, + y=yvar, + levels=[0.15, 0.4, 0.6, 0.8, 1], + transform=ccrs.PlateCarree(), + cmap=cmap, + add_colorbar=False, + ) ) antarctic_regions.plot_regions( ax=ax, @@ -416,17 +465,27 @@ def create_antarctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, m label="abbrev", line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, ) + if pos2[mo] == 0: + ax.text( + 0.03, + 0.04, + model.replace("_", " "), + horizontalalignment="left", + verticalalignment="center", + transform=ax.transAxes, + backgroundcolor="white", + fontsize=14, + ) ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) ax.coastlines(color=[0.3, 0.3, 0.3]) ax.set_facecolor([0.55, 0.55, 0.6]) - ax.set_title(months[mo] + "\n" + model.replace("_", " "), fontsize=12) - pos1 = [1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0] - pos2 = [0, 0, 2, 2, 2, 4, 4, 4, 6, 6, 6, 0] + ax.set_title("") + pos1 = [0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2] # row + pos2 = [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5] # col xvar = lib.find_lon(obs) yvar = lib.find_lat(obs) for mo in range(0, 12): ax = axs[pos1[mo], pos2[mo]] - # ax.set_global() obs[var_obs].isel({"time": mo}).plot( ax=ax, x=xvar, @@ -434,7 +493,7 @@ def create_antarctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, m levels=[0.15, 0.4, 0.6, 0.8, 1], transform=ccrs.PlateCarree(), cmap=cmap, - cbar_kwargs={"label": "ice fraction", "fraction": 0.046, "pad": 0.04}, + add_colorbar=False, ) antarctic_regions.plot_regions( ax=ax, @@ -442,15 +501,41 @@ def create_antarctic_map(ds, obs, var_ice, var_obs, metrics_output_path, meta, m label="abbrev", line_kws={"color": [0.2, 0.2, 0.25], "linewidth": 3}, ) + if pos2[mo] == 0: + ax.text( + 0.03, + 0.04, + "Reference", + horizontalalignment="left", + verticalalignment="center", + transform=ax.transAxes, + backgroundcolor="white", + fontsize=14, + ) ax.set_extent([-180, 180, -53, -90], ccrs.PlateCarree()) ax.coastlines(color=[0.3, 0.3, 0.3]) ax.set_facecolor([0.55, 0.55, 0.6]) - ax.set_title(months[mo] + "\nReference", fontsize=12) - plt.suptitle(model.replace("_", " "), fontsize=20) + ax.set_title(months[mo], fontsize=18) + plt.suptitle(title, fontsize=26) + cbar = plt.colorbar( + fds, + ax=axs[0:2, :], + pad=0.02, + aspect=25, + ) + cbar.ax.tick_params(labelsize=18) + cbar.set_label(label="ice fraction", size=18) + cbar = plt.colorbar( + fds, + ax=axs[2:, :], + pad=0.02, + aspect=25, + ) + cbar.ax.tick_params(labelsize=18) + cbar.set_label(label="ice fraction", size=18) fig_path = os.path.join( metrics_output_path, model.replace(" ", "_") + "_antarctic_regions.png" ) - plt.tight_layout(rect=[0, 0.03, 1, 0.97]) plt.savefig(fig_path) plt.close() meta.update_plots( From 1dbcfc0b57ada64c66d2da4b266c22ae209425f1 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 5 Jul 2024 16:18:50 -0700 Subject: [PATCH 14/14] update figures --- pcmdi_metrics/sea_ice/sea_ice_driver.py | 37 ++++++++++++------------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/pcmdi_metrics/sea_ice/sea_ice_driver.py b/pcmdi_metrics/sea_ice/sea_ice_driver.py index ddb727152..d400e4c63 100644 --- a/pcmdi_metrics/sea_ice/sea_ice_driver.py +++ b/pcmdi_metrics/sea_ice/sea_ice_driver.py @@ -596,7 +596,8 @@ if model == reference_data_set: continue try: - tmp_title = "_".join([model, run]) + tmp_model = "_".join([model, run]) + tmp_title = "{0}-{1} Arctic sea ice".format(yr_range[0], yr_range[1]) meta = fig.create_arctic_map( nc_climo, obs_nh, @@ -604,8 +605,10 @@ obs_var, fig_dir, meta, + tmp_model, tmp_title, ) + tmp_title = "{0}-{1} Antarctic sea ice".format(yr_range[0], yr_range[1]) meta = fig.create_antarctic_map( nc_climo, obs_sh, @@ -613,6 +616,7 @@ obs_var, fig_dir, meta, + tmp_model, tmp_title, ) except Exception as e: @@ -630,29 +634,22 @@ else: nc_climo_mean[var] = nc_climo_mean[var] + nc_climo[var] nc_climo_mean[var] = nc_climo_mean[var] / (count + 1) - # try: - tmp_title = "_".join([model, "model_mean"]) + tmp_model = "_".join([model, "model_mean"]) + tmp_title = "{0}-{1} Arctic sea ice".format(yr_range[0], yr_range[1]) meta = fig.create_arctic_map( - nc_climo_mean, - obs_nh, - var, - obs_var, - fig_dir, - meta, - tmp_title, + nc_climo_mean, obs_nh, var, obs_var, fig_dir, meta, tmp_model, tmp_title ) + tmp_title = "{0}-{1} Antarctic sea ice".format(yr_range[0], yr_range[1]) meta = fig.create_antarctic_map( - nc_climo_mean, - obs_sh, - var, - obs_var, - fig_dir, - meta, - tmp_title, + nc_climo_mean, obs_sh, var, obs_var, fig_dir, meta, tmp_model, tmp_title + ) + + meta = fig.create_summary_maps_arctic( + nc_climo_mean, var, fig_dir, meta, tmp_model + ) + meta = fig.create_summary_maps_antarctic( + nc_climo_mean, var, fig_dir, meta, tmp_model ) - # except Exception as e: - # print("Error making figures for model",model,"mean.") - # print(" ",e) # ----------------- # Update and write