Merge pull request #1024 from PCMDI/405_sic_ao

Sea ice metrics beta

conda-env/dev.yml

@@ -29,6 +29,7 @@ dependencies:
     - rasterio=1.3.6
     - shapely=2.0.1
     - numdifftools
+    - nc-time-axis
     # ==================
     # Testing
     # ==================

doc/jupyter/Demo/sea_ice_line_plots.py

@@ -0,0 +1,61 @@
+import cftime
+import matplotlib.pyplot as plt
+import xarray as xr
+import xcdat as xc
+
+ds = xc.open_mfdataset(
+    "/p/user_pub/pmp/demo/sea-ice/links_siconc/E3SM-1-0/historical/r1i2p2f1/siconc/siconc_SImon_E3SM-1-0_historical_r1i2p2f1_*_*.nc"
+)
+area = xc.open_dataset(
+    "/p/user_pub/pmp/demo/sea-ice/links_area/E3SM-1-0/areacello_Ofx_E3SM-1-0_historical_r1i1p1f1_gr.nc"
+)
+
+arctic = (ds.where(ds.lat > 0) * 1e-2 * area.areacello * 1e-6).sum(("lat", "lon"))
+
+f_os_n = "/p/user_pub/pmp/demo/sea-ice/EUMETSAT/OSI-SAF-450-a-3-0/v20231201/ice_conc_nh_ease2-250_cdr-v3p0_198801-202012.nc"
+obs = xc.open_dataset(f_os_n)
+obs_area = 625
+obs_arctic = (obs.ice_conc.where(obs.lat > 0) * 1e-2 * obs_area).sum(("xc", "yc"))
+
+# Time series plot
+arctic.siconc.sel({"time": slice("1981-01-01", "2010-12-31")}).plot(label="E3SM-1-0")
+obs_arctic.plot(label="OSI-SAF")
+plt.title("Arctic monthly sea ice extent")
+plt.ylabel("Extent (km${^2}$)")
+plt.xlabel("time")
+plt.xlim(
+    [
+        cftime.DatetimeNoLeap(1981, 1, 16, 12, 0, 0, 0, has_year_zero=True),
+        cftime.DatetimeNoLeap(2010, 12, 16, 12, 0, 0, 0, has_year_zero=True),
+    ]
+)
+plt.legend(loc="upper right", fontsize=9)
+plt.savefig("Arctic_tseries.png")
+plt.close()
+
+# Climatology plot
+arctic_ds = xr.Dataset(
+    data_vars={"siconc": arctic.siconc, "time_bnds": ds.time_bnds},
+    coords={"time": ds.time},
+)
+arctic_clim = arctic_ds.sel(
+    {"time": slice("1981-01-01", "2010-12-31")}
+).temporal.climatology("siconc", freq="month")
+arctic_clim["time"] = [x for x in range(1, 13)]
+
+obs_arc_ds = xr.Dataset(
+    data_vars={"ice_conc": obs_arctic, "time_bnds": obs.time_bnds},
+    coords={"time": obs.time},
+)
+obs_clim = obs_arc_ds.temporal.climatology("ice_conc", freq="month")
+obs_clim["time"] = [x for x in range(1, 13)]
+
+arctic_clim.siconc.plot(label="E3SM-1-0")
+obs_clim.ice_conc.plot(label="OSI-SAF")
+plt.title("Arctic climatological sea ice extent\n1981-2010")
+plt.xlabel("month")
+plt.ylabel("Extent (km${^2}$)")
+plt.xlim([1, 12])
+plt.legend(loc="upper right", fontsize=9)
+plt.savefig("Arctic_clim.png")
+plt.close()

doc/jupyter/Demo/sea_ice_param.py

@@ -0,0 +1,52 @@
+# Sea ice metrics parameter file
+
+# List of models to include in analysis
+test_data_set = [
+    "E3SM-1-0",
+]
+
+# realization can be a single realization, a list of realizations, or "*" for all realizations
+realization = "r1i2p2f1"
+
+# test_data_path is a template for the model data parent directory
+test_data_path = "/p/user_pub/pmp/demo/sea-ice/links_siconc/%(model)/historical/%(realization)/siconc/"
+
+# filename_template is a template for the model data file name
+# combine it with test_data_path to get complete data path
+filename_template = "siconc_SImon_%(model)_historical_%(realization)_*_*.nc"
+
+# The name of the sea ice variable in the model data
+var = "siconc"
+
+# Start and end years for model data
+msyear = 1981
+meyear = 2010
+
+# Factor for adjusting model data to decimal rather than percent units
+ModUnitsAdjust = (True, "multiply", 1e-2)
+
+# Template for the grid area file
+area_template = "/p/user_pub/pmp/demo/sea-ice/links_area/%(model)/*.nc"
+
+# Area variable name; likely 'areacello' or 'areacella' for CMIP6
+area_var = "areacello"
+
+# Factor to convert area units to km-2
+AreaUnitsAdjust = (True, "multiply", 1e-6)
+
+# Directory for writing outputs
+case_id = "ex1"
+metrics_output_path = "sea_ice_demo/%(case_id)/"
+
+# Settings for the observational data
+reference_data_path_nh = "/p/user_pub/pmp/demo/sea-ice/EUMETSAT/OSI-SAF-450-a-3-0/v20231201/ice_conc_nh_ease2-250_cdr-v3p0_198801-202012.nc"
+reference_data_path_sh = "/p/user_pub/pmp/demo/sea-ice/EUMETSAT/OSI-SAF-450-a-3-0/v20231201/ice_conc_sh_ease2-250_cdr-v3p0_198801-202012.nc"
+ObsUnitsAdjust = (True, "multiply", 1e-2)
+reference_data_set = "OSI-SAF"
+osyear = 1988
+oeyear = 2020
+obs_var = "ice_conc"
+ObsAreaUnitsAdjust = (False, 0, 0)
+obs_area_template = None
+obs_area_var = None
+obs_cell_area = 625  # km 2

doc/jupyter/Demo/sea_ice_sector_plots.py

@@ -0,0 +1,156 @@
+import cartopy.crs as ccrs
+import matplotlib.colors as colors
+import matplotlib.pyplot as plt
+import numpy as np
+import regionmask
+import xcdat as xc
+
+from pcmdi_metrics.utils import create_land_sea_mask
+
+# ----------
+# Arctic
+# ----------
+print("Creating Arctic map")
+# Load and process data
+f_os_n = "/p/user_pub/pmp/demo/sea-ice/EUMETSAT/OSI-SAF-450-a-3-0/v20231201/ice_conc_nh_ease2-250_cdr-v3p0_198801-202012.nc"
+obs = xc.open_dataset(f_os_n)
+obs = obs.mean("time")
+mask = create_land_sea_mask(obs, lon_key="lon", lat_key="lat")
+obs["ice_conc"] = obs["ice_conc"].where(mask < 1)
+ds = obs.assign_coords(
+    xc=obs["lon"], yc=obs["lat"]
+)  # Assign these variables to Coordinates, which were originally data variables
+
+# 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()
+ds.ice_conc.plot.pcolormesh(
+    ax=ax,
+    x="xc",
+    y="yc",
+    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])
+plt.title("Arctic regions with mean\nOSI-SAF ice concentration\n1988-2020")
+plt.savefig("Arctic_regions.png")
+plt.close()
+obs.close()
+
+# ----------
+# Antarctic
+# ----------
+print("Creating Antarctic map")
+# Load and process data
+f_os_s = "/p/user_pub/pmp/demo/sea-ice/EUMETSAT/OSI-SAF-450-a-3-0/v20231201/ice_conc_sh_ease2-250_cdr-v3p0_198801-202012.nc"
+obs = xc.open_dataset(f_os_s)
+obs = obs.mean("time")
+mask = create_land_sea_mask(obs, lon_key="lon", lat_key="lat")
+obs["ice_conc"] = obs["ice_conc"].where(mask < 1)
+ds = obs.assign_coords(
+    xc=obs["lon"], yc=obs["lat"]
+)  # Assign these variables to Coordinates, which were originally data variables
+
+# 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.ice_conc.plot.pcolormesh(
+    ax=ax,
+    x="xc",
+    y="yc",
+    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])
+plt.title("Antarctic regions with mean\nOSI-SAF ice concentration\n1988-2020")
+plt.savefig("Antarctic_regions.png")
+plt.close()
+obs.close()

pcmdi_metrics/sea_ice/README.md

@@ -0,0 +1,61 @@
+# Sea Ice
+
+## Summary
+
+## Demo
+
+* Link to notebook
+
+## Inputs
+
+### Sectors
+
+## Run
+
+The PMP sea ice metrics can be controlled via an input parameter file, the command line, or both. With the command line only it is executed via:
+
+```
+sea_ice_driver.py -p parameter_file.py
+```
+
+or as a combination of an input parameter file and the command line, e.g.:
+
+```
+sea_ice_driver.py -p parameter_file.py --msyear 1991 --meyear 2020
+```
+
+## Outputs
+
+The driver produces a JSON file containing mean square error metrics for all input models and realizations relative to the reference data set. It also produces a bar chart displaying these metrics.
+
+## Parameters
+
+* **case_id**: Save JSON and figure files into this subdirectory so that results from multiple tests can be readily organized.
+* **test_data_set**: List of model names.
+* **realization**: List of realizations.
+* **test_data_path**: File path to directory containing model/test data.
+* **filename_template**: File name template for test data, e.g., "CMIP5.historical.%(model_version).r1i1p1.mon.%(variable).19810-200512.AC.v2019022512.nc" where "model_version" and "variable" will be analyzed for each of the entries in test_data_set and vars.
+* **var**: Name of model sea ice variable
+* **msyear**: Start year for test data set.
+* **meyear**: End year for test data set.
+* **ModUnitsAdjust**: Factor to convert model sea ice data to fraction of 1. Uses format (flag (bool), operation (str), value (float)). Operation can be "add", "subtract", "multiply", or "divide". For example, use (True, 'multiply', 1e-2) to convert from percent concentration to decimal concentration.
+* **area_template**: File path of model grid area data.
+* **area_var**: Name of model area variable, e.g. "areacello"
+* **AreaUnitsAdjust**: Factor to convert model area data to units of km2. Uses format (flag (bool), operation (str), value (float)). Operation can be "add", "subtract", "multiply", or "divide". For example, use (True, 'multiply', 1e6) to convert from m2 to km2.
+* **metrics_output_path**: Directory path for metrics output in JSON files, e.g., '~/demo_data/PMP_metrics/'. The %(case_id) variable can be used here. If exists, should be empty before run.
+* **reference_data_path_nh**: The reference data file path for the northern hemisphere. If data is global, provide same path for nh and sh.
+* **reference_data_path_sh**: The reference data file path for the southern hemisphere. If data is global, provide same path for nh and sh.
+* **ObsUnitsAdjust**: Factor to convert reference sea ice data to fraction of 1. Uses format (flag (bool), operation (str), value (float)). Operation can be "add", "subtract", "multiply", or "divide". For example, use (True, 'multiply', 1e-2) to convert from percent concentration to decimal concentration.
+* **reference_data_set**: A short name describing the reference dataset, e.g. "OSI-SAF".
+* **osyear**: Start year for reference data set.
+* **oeyear**: End year for reference data set.
+* **obs_var**: Name of reference sea ice variable.
+* **ObsAreaUnitsAdjust**: Factor to convert model area data to units of km2. Uses format (flag (bool), operation (str), value (float)). Operation can be "add", "subtract", "multiply", or "divide". For example, use (True, 'multiply', 1e6) to convert from m2 to km2.
+* **obs_area_template**: File path of grid area data. If unavailalbe, skip and use "obs_cell_area".
+* **obs_area_var**: Name of reference area variable, if available. If unavailable, skip and use "obs_cell_area".
+* **obs_cell_area**: For equal area grids, the area of a single grid cell in units of km2.
+
+
+## Reference
+
+Ivanova, D. P., P. J. Gleckler, K. E. Taylor, P. J. Durack, and K. D. Marvel, 2016: Moving beyond the Total Sea Ice Extent in Gauging Model Biases. J. Climate, 29, 8965–8987, https://doi.org/10.1175/JCLI-D-16-0026.1. 

pcmdi_metrics/sea_ice/lib/__init__.py

@@ -0,0 +1 @@
+from .sea_ice_parser import create_sea_ice_parser