Merge pull request #90 from rubisco-sfa/ngeea

ADD: changes needed for NGEE-A

Author: nocollier

ilamb3/__init__.py

@@ -20,6 +20,7 @@
 units.define("Tg = 1e12 * g")
 units.define("Pg = 1e15 * g")
 units.define("Sv = 1e6 m**3 / s")
+units.define("gC = g")
 units.define("@alias W = Watt")
 
 # we don't really have data versions for the collection :/

ilamb3/analysis/spatial_distribution.py

@@ -116,6 +116,9 @@ def __call__(
             # Spatial standard deviation
             ref_std = float(rref.std())
             com_std = float(rcom.std())
+            if np.allclose(ref_std, 0):
+                # There is no spatial variance for this region and we should skip
+                continue
             norm_std = com_std / ref_std
 
             # Correlation
@@ -170,6 +173,7 @@ def plots(
     ) -> pd.DataFrame:
         # Add the Taylor diagram
         axs = []
+        df = df[df["analysis"] == "Spatial Distribution"]
         for region in df["region"].unique():
             axs += [
                 {
@@ -180,5 +184,5 @@ def plots(
                     "axis": ilplt.plot_taylor_diagram(df[df["region"] == region]),
                 }
             ]
-        axs = pd.DataFrame(axs).dropna(subset=["axis"])
+        axs = pd.DataFrame(axs)
         return axs

ilamb3/compare.py

@@ -295,8 +295,11 @@ def handle_timescale_mismatch(
             }
         )
         return ref, com
+    if np.allclose(dt_ref, 30, atol=3) and dt_com < 28:
+        com = dset.compute_monthly_mean(com)
+        return ref, com
     raise NotImplementedError(
-        f"We encountered a time scale mismatch that we have no logic to handle. {ref} {com}"
+        f"We encountered a time scale mismatch ({dt_ref=:.1f} vs. {dt_com=:.1f}) that we have no logic to handle. {ref} {com}"
     )
 
 
@@ -327,8 +330,10 @@ def make_comparable(
         com = com.sortby([dset.get_dim_name(com, "lat"), dset.get_dim_name(com, "lon")])
 
     # pick just the sites
-    if dset.is_site(ref[varname]):
+    if dset.is_site(ref[varname]) and dset.is_spatial(com[varname]):
         com = extract_sites(ref, com, varname)
+    if dset.is_site(com[varname]) and dset.is_spatial(ref[varname]):
+        ref = extract_sites(com, ref, varname)
 
     # convert units
     com = dset.convert(com, ref[varname].attrs["units"], varname=varname)

ilamb3/dataset.py

@@ -864,18 +864,17 @@ def convert(
     xr.Dataset or xr.DataArray
         The converted dataset.
     """
-    dset = dset.pint.quantify()
     if isinstance(dset, xr.DataArray):
         da = dset
     else:
         assert varname is not None
         da = dset[varname]
+    da = da.pint.quantify()
     da = scale_by_water_density(da, unit)
     da = da.pint.to(unit)
     if isinstance(dset, xr.DataArray):
         return da.pint.dequantify()
-    dset[varname] = da
-    dset = dset.pint.dequantify()
+    dset[varname] = da.pint.dequantify()
     return dset
 
 
@@ -1007,3 +1006,18 @@ def compute_monthly_mean(ds: xr.Dataset) -> xr.Dataset:
         raise ValueError(f"Input dataset is already coarser than monthly {dt=}.")
     ds = ds.resample(time=xr.groupers.TimeResampler("MS")).mean()
     return ds
+
+
+def shift_time_by_years(ds: xr.Dataset, years: int) -> xr.Dataset:
+    """
+    Return the dataset shifted by the given number of years.
+    """
+    time_dim = get_dim_name(ds, "time")
+    cls = ds[time_dim].values[0].__class__
+    ds[time_dim] = [
+        cls(
+            t.dt.year + years, t.dt.month, t.dt.day, t.dt.hour, t.dt.minute, t.dt.second
+        )
+        for t in ds[time_dim]
+    ]
+    return ds

ilamb3/run.py

@@ -49,6 +49,21 @@ def _fix(units: str) -> str:
     return ds
 
 
+def fix_lndgrid_coords(ds: xr.Dataset) -> xr.Dataset:
+    """
+    Return a dataset with coordinates properly assigned.
+
+    Note
+    ----
+    E3SM/CESM2 raw land model output comes as if it were run over sites in the
+    `lndgrid` dimension. Some of the variables that are listed in these files
+    are only of this dimension and really belong in the coordinates. These tend
+    to be things like `lat`, `lon`, etc. and ilamb3 needs them to be associated
+    with the dataset coordinates to work.
+    """
+    return ds.assign_coords({v: ds[v] for v in ds if ds[v].dims == ("lndgrid",)})
+
+
 def select_analysis_variable(setup: dict[str, Any]) -> str:
     """
     Return the main variable to be used in this analysis.
@@ -358,6 +373,8 @@ def _load_comparison_data(
     # Fix bounds attributes (there is a bounds variable but it isn't in the
     # attributes)
     com = {var: dset.fix_missing_bounds_attrs(ds) for var, ds in com.items()}
+    # Fix lndgrid cooridates if raw E3SM/CESM2 data is given
+    com = {var: fix_lndgrid_coords(ds) for var, ds in com.items()}
     # Merge all the data together
     if len(com) > 1:
         # The grids should be the same, but sometimes models generate output
@@ -547,21 +564,34 @@ def run_single_block(
             "Reference intermediate data was not generated."
         )  # pragma: no cover
 
+    log_file = output_path / "post.log"
+    log_id = logger.add(log_file, backtrace=True, diagnose=True)
+
     # Phase 2: get plots and combine scalars and save
-    plt.rcParams.update({"figure.max_open_warning": 0})
-    df = pd.concat(df_all).drop_duplicates(
-        subset=["source", "region", "analysis", "name"]
-    )
-    df = add_overall_score(df)
-    df_plots = plot_analyses(df, ds_ref, ds_com, analyses, output_path)
+    try:
+        plt.rcParams.update({"figure.max_open_warning": 0})
+        df = pd.concat(df_all).drop_duplicates(
+            subset=["source", "region", "analysis", "name"]
+        )
+        df = add_overall_score(df)
+        df_plots = plot_analyses(df, ds_ref, ds_com, analyses, output_path)
+    except Exception:
+        logger.exception(f"ILAMB analysis '{block_name}' failed in plotting.")
+        return
 
     # Generate an output page
-    if ilamb3.conf["debug_mode"] and (output_path / "index.html").is_file():
+    try:
+        if ilamb3.conf["debug_mode"] and (output_path / "index.html").is_file():
+            logger.remove(log_id)
+            return
+        ds_ref.attrs["header"] = block_name
+        html = generate_html_page(df, ds_ref, ds_com, df_plots)
+        with open(output_path / "index.html", mode="w") as out:
+            out.write(html)
+    except Exception:
+        logger.exception(f"ILAMB analysis '{block_name}' failed in generating html.")
         return
-    ds_ref.attrs["header"] = block_name
-    html = generate_html_page(df, ds_ref, ds_com, df_plots)
-    with open(output_path / "index.html", mode="w") as out:
-        out.write(html)
+    logger.remove(log_id)
 
 
 def run_analyses(

ilamb3/templates/directory.html

@@ -0,0 +1,15 @@
+<!doctype html>
+<html lang="en">
+  <head>
+    <meta charset="utf-8">
+    <meta name="viewport" content="width=device-width, initial-scale=1">
+    <title>Directory of Benchmarking Results</title>
+  </head>
+  <body>
+    <ul>
+{% for name,link in links.items() %}
+<li><a href="{{link}}">{{name}}</a></li>
+{% endfor %}
+    </ul>
+  </body>
+</html>