bug fix (continue)

pcmdi_metrics/variability_mode/lib/calc_stat.py

@@ -1,5 +1,7 @@
 from time import gmtime, strftime
 
+import xarray as xr
+
 from pcmdi_metrics.io import get_grid, region_subset
 from pcmdi_metrics.stats import bias_xy as calcBias
 from pcmdi_metrics.stats import cor_xy as calcSCOR
@@ -10,14 +12,14 @@
 
 
 def calc_stats_save_dict(
- mode,
- dict_head,
- eof,
- eof_lr,
+ mode: str,
+ dict_head: dict,
+ eof: xr.Dataset,
+ eof_lr: xr.Dataset,
  pc,
  stdv_pc,
  frac,
- regions_specs,
+ regions_specs: dict = None,
  eof_obs=None,
  eof_lr_obs=None,
  stdv_pc_obs=None,
@@ -28,6 +30,7 @@ def calc_stats_save_dict(
  """
  NOTE: Calculate statistics and save numbers to dictionary for JSON.
  Input
+ - mode: [str] name of variability mode
  - dict_head: [dict] subset of dictionary
  - eof: [2d field] linear regressed eof pattern (eof domain)
  - eof_lr: [2d field] linear regressed eof pattern (global)
@@ -60,18 +63,17 @@ def calc_stats_save_dict(
  # Note: '_glo' indicates statistics calculated over global domain
  # . . . . . . . . . . . . . . . . . . . . . . . . .
  if obs_compare:
- if method in ["eof", "cbf"]:
- ref_grid_global = get_grid(eof_lr_obs)
- # Regrid (interpolation, model grid to ref grid)
- debug_print("regrid (global) start", debug)
- # eof_model_global = eof_lr.regrid(eof_lr,
- # ref_grid_global, regridTool="regrid2", mkCyclic=True
- # )
- eof_model_global = regrid(eof_lr, ref_grid_global)
- debug_print("regrid end", debug)
- # Extract subdomain
- # eof_model = eof_model_global(region_subdomain)
- eof_model = region_subset(eof_model_global, mode, regions_specs)
+ ref_grid_global = get_grid(eof_lr_obs)
+ # Regrid (interpolation, model grid to ref grid)
+ debug_print("regrid (global) start", debug)
+ # eof_model_global = eof_lr.regrid(eof_lr,
+ # ref_grid_global, regridTool="regrid2", mkCyclic=True
+ # )
+ eof_model_global = regrid(eof_lr, ref_grid_global)
+ debug_print("regrid end", debug)
+ # Extract subdomain
+ # eof_model = eof_model_global(region_subdomain)
+ eof_model = region_subset(eof_model_global, mode, regions_specs=regions_specs)
 
  # Spatial correlation weighted by area ('generate' option for weights)
  cor = calcSCOR(eof_model, eof_obs)

pcmdi_metrics/variability_mode/lib/eof_analysis.py

@@ -174,22 +174,23 @@ def arbitrary_checking(mode, eof_Nth):
 
 
 def linear_regression_on_globe_for_teleconnection(
- pc, model_timeseries, stdv_pc, RmDomainMean=True, EofScaling=False, debug=False
+ pc, ds, data_var, stdv_pc, RmDomainMean=True, EofScaling=False, debug=False
 ):
  """
  - Reconstruct EOF fist mode including teleconnection purpose as well
  - Have confirmed that "eof_lr" is identical to "eof" over EOF domain (i.e., "subdomain")
  - Note that eof_lr has global field
  """
  if debug:
- print("pc.shape, timeseries.shape:", pc.shape, model_timeseries.shape)
+ print("pc.shape, timeseries.shape:", pc.shape, ds[data_var].shape)
 
  # Linear regression to have extended global map; teleconnection purpose
- slope, intercept = linear_regression(pc, model_timeseries)
+ slope, intercept = linear_regression(pc, ds[data_var])
 
- factor = stdv_pc
  if not RmDomainMean and EofScaling:
  factor = 1
+ else:
+ factor = stdv_pc
 
  eof_lr = (slope * factor) + intercept
 

pcmdi_metrics/variability_mode/lib/lib_variability_mode.py

@@ -217,10 +217,14 @@ def diff_month(date1, date2):
  return (date2.year - date1.year) * 12 + date2.month - date1.month + 1
 
 
-def debug_print(string, debug):
+def debug_print(to_check, debug):
  if debug:
  nowtime = strftime("%Y-%m-%d %H:%M:%S", gmtime())
- print("debug: " + nowtime + " " + string)
+ if isinstance(to_check, str):
+ print("debug: " + nowtime + " " + to_check)
+ else:
+ print("debug: " + nowtime)
+ print(to_check)
 
 
 def pick_year_last_day(ds):

pcmdi_metrics/variability_mode/lib/plot_map.py

@@ -19,7 +19,15 @@
 
 
 def plot_map(
- mode, model, syear, eyear, season, eof_Nth, frac_Nth, output_file_name, debug=False
+ mode: str,
+ model: str,
+ syear: int,
+ eyear: int,
+ season: str,
+ eof_pattern: xr.DataArray,
+ eof_variance_fraction: float,
+ output_file_name: str,
+ debug=False,
 ):
  """Plot dive down map and save
 
@@ -35,9 +43,9 @@ def plot_map(
  End year from analysis
  season : str
  season ("DJF", "MAM", "JJA", "SON", "monthly", or "yearly") that was used for analysis and will be shown in figure title
- eof_Nth : cdms2.TransientVariable
+ eof_pattern : cdms2.TransientVariable
  EOF pattern to plot, 2D cdms2 TransientVariable with lat/lon coordinates attached
- frac_Nth : float
+ eof_variance_fraction : float
  Fraction of explained variability (0 to 1), which will be shown in the figure as percentage after multiplying 100
  output_file_name : str
  Name of output image file (e.g., "output_file.png")
@@ -55,26 +63,14 @@ def plot_map(
  sys.exit("Projection for " + mode + "is not defined.")
 
  # title
- if frac_Nth != -999:
+ if eof_variance_fraction != -999:
  percentage = (
- str(round(float(frac_Nth * 100.0), 1)) + "%"
+ str(round(float(eof_variance_fraction * 100.0), 1)) + "%"
  ) # % with one floating number
  else:
  percentage = ""
 
- plot_title = (
- mode
- + ": "
- + model
- + "\n"
- + str(syear)
- + "-"
- + str(eyear)
- + " "
- + season
- + " "
- + percentage
- )
+ plot_title = f"{mode}: {model}\n{syear}-{eyear} {season} {percentage}"
 
  debug_print(
  "plot_map: projection, plot_title:" + projection + ", " + plot_title, debug
@@ -102,13 +98,7 @@ def plot_map(
  central_longitude = 180
 
  # Convert cdms variable to xarray
- lon = get_longitude(eof_Nth)
- lat = get_latitude(eof_Nth)
- data = np.array(eof_Nth)
- lon, lat = np.meshgrid(lon, lat)
- data_array = xr.DataArray(
- data, coords={"lon": lon[0, :], "lat": lat[:, 0]}, dims=("lat", "lon")
- )
+ data_array = eof_pattern
  data_array = data_array.where(data_array != 1e20, np.nan)
 
  plot_map_cartopy(
@@ -164,8 +154,8 @@ def plot_map_cartopy(
 
  debug_print("plot_map_cartopy starts", debug)
 
- lon = data_array.lon
- lat = data_array.lat
+ lon = get_longitude(data_array)
+ lat = get_latitude(data_array)
 
  # Determine the extent based on the longitude range where data exists
  lon_min = lon.min().item()

pcmdi_metrics/variability_mode/variability_modes_driver.py

@@ -387,12 +387,22 @@
  intercept_obs,
  ) = linear_regression_on_globe_for_teleconnection(
  pc_obs[season],
- obs_timeseries_season[obs_var],
+ obs_timeseries_season,
+ obs_var,
  stdv_pc_obs[season],
  RmDomainMean,
  EofScaling,
  debug=debug,
  )
+
+ obs_timeseries_season["eof_lr"] = eof_lr_obs_season
+ obs_timeseries_season["slope"] = slope_obs
+ obs_timeseries_season["intercept"] = intercept_obs
+
+ # Extract subdomain for plot
+ obs_timeseries_season_region = region_subset(
+ obs_timeseries_season, mode, regions_specs=regions_specs
+ )
  eof_lr_obs[season] = eof_lr_obs_season
  # - - - - - - - - - - - - - - - - - - - - - - - - -
  # Record results
@@ -413,16 +423,14 @@
  output_img_file_obs = os.path.join(
  dir_paths["graphics"], output_filename_obs
  )
- eof_lr_obs_season_region = region_subset(
- eof_lr_obs_season, mode, regions_specs
- )
+
  plot_map(
  mode,
  "[REF] " + obs_name,
  osyear,
  oeyear,
  season,
- eof_lr_obs_season_region,
+ obs_timeseries_season_region["eof_lr"],
  frac_obs[season],
  output_img_file_obs,
  debug=debug,
@@ -434,7 +442,7 @@
  oeyear,
  season,
  # eof_lr_obs[season](longitude=(lon1_global, lon2_global)),
- eof_lr_obs[season],
+ eof_lr_obs_season,
  frac_obs[season],
  output_img_file_obs + "_teleconnection",
  debug=debug,
@@ -663,36 +671,39 @@
  intercept_cbf,
  ) = linear_regression_on_globe_for_teleconnection(
  cbf_pc,
- model_timeseries_season[var],
+ model_timeseries_season,
+ var,
  stdv_cbf_pc,
- # cbf_pc, model_timeseries_season_regrid, stdv_cbf_pc,
  RmDomainMean,
  EofScaling,
  debug=debug,
  )
 
+ model_timeseries_season["eof_lr_cbf"] = eof_lr_cbf
+ model_timeseries_season["slope_cbf"] = slope_cbf
+ model_timeseries_season["intercept_cbf"] = intercept_cbf
+
  # Extract subdomain for statistics
  # eof_lr_cbf_subdomain = eof_lr_cbf(region_subdomain)
- eof_lr_cbf_subdomain = region_subset(
- eof_lr_cbf, mode, regions_specs
+ model_timeseries_season_subdomain = region_subset(
+ model_timeseries_season,
+ mode,
+ regions_specs=regions_specs,
  )
 
- # Calculate fraction of variance explained by cbf pc
+ # Calculate fraction of variance explained by cbf pc (native grid)
  frac_cbf = gain_pcs_fraction(
- # model_timeseries_season_regrid_subdomain, # regridded model anomaly space
- model_timeseries_season_subdomain[
- var
- ], # native grid model anomaly space
- eof_lr_cbf_subdomain,
+ model_timeseries_season_subdomain[var],
+ model_timeseries_season_subdomain["eof_lr_cbf"],
  cbf_pc / stdv_cbf_pc,
  debug=debug,
  )
 
  # SENSITIVITY TEST ---
- # Calculate fraction of variance explained by cbf pc (on regrid domain)
+ # Calculate fraction of variance explained by cbf pc (regrid domain)
  frac_cbf_regrid = gain_pcs_fraction(
  model_timeseries_season_regrid_subdomain[var],
- eof_lr_cbf_subdomain,
+ model_timeseries_season_subdomain["eof_lr_cbf"],
  cbf_pc / stdv_cbf_pc,
  debug=debug,
  )
@@ -705,7 +716,7 @@
  dict_head, eof_lr_cbf = calc_stats_save_dict(
  mode,
  dict_head,
- eof_lr_cbf_subdomain,
+ model_timeseries_season_subdomain["eof_lr_cbf"],
  eof_lr_cbf,
  cbf_pc,
  stdv_cbf_pc,
@@ -750,7 +761,7 @@
  msyear,
  meyear,
  season,
- eof_lr_cbf_subdomain,
+ model_timeseries_season_subdomain["eof_lr_cbf"],
  frac_cbf,
  output_img_file + "_cbf",
  debug=debug,
@@ -846,6 +857,7 @@
  # Metrics results -- statistics to JSON
  if obs_compare:
  dict_head, eof_lr = calc_stats_save_dict(
+ mode,
  dict_head,
  eof,
  eof_lr,
@@ -862,6 +874,7 @@
  )
  else:
  dict_head, eof_lr = calc_stats_save_dict(
+ mode,
  dict_head,
  eof,
  eof_lr,