update

conda-env/dev.yml

@@ -11,7 +11,7 @@ dependencies:
  - python=3.10.10
  - pip=23.1.2
  - numpy=1.23.5
- - cartopy=0.21.1
+ - cartopy=0.22.0
  - matplotlib=3.7.1
  - cdat_info=8.2.1
  - cdms2=3.1.5
@@ -27,7 +27,7 @@ dependencies:
  - netcdf4=1.6.3
  - regionmask=0.9.0
  - rasterio=1.3.6
- - shapely=2.0.1
+ - shapely=2.0.2
  # ==================
  # Testing
  # ==================

pcmdi_metrics/variability_mode/lib/argparse_functions.py

@@ -150,13 +150,25 @@ def AddParserArgument(P):
  "--nc_out",
  type=bool,
  default=True,
- help="Option for generate netCDF file output: True (default) / False",
+ help="Option for generate netCDF file output for models: True (default) / False",
  )
  P.add_argument(
  "--plot",
  type=bool,
  default=True,
- help="Option for generate individual plots: True (default) / False",
+ help="Option for generate individual plots for models: True (default) / False",
+ )
+ P.add_argument(
+ "--nc_out_obs",
+ type=bool,
+ default=True,
+ help="Option for generate netCDF file output for obs: True (default) / False",
+ )
+ P.add_argument(
+ "--plot_obs",
+ type=bool,
+ default=True,
+ help="Option for generate individual plots for obs: True (default) / False",
  )
  P.add_argument(
  "--parallel",

pcmdi_metrics/variability_mode/lib/plot_map.py

@@ -1,16 +1,23 @@
 import sys
 
-import cartopy
+#import cartopy
 import cartopy.crs as ccrs
 import matplotlib.path as mpath
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mticker
 import numpy as np
 from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
 from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
+from pcmdi_metrics.variability_mode.lib import debug_print
+from cartopy.feature import LAND as cartopy_land
+from cartopy.feature import OCEAN as cartopy_ocean
 
+import faulthandler
 
-def plot_map(mode, model, syear, eyear, season, eof_Nth, frac_Nth, output_file_name):
+faulthandler.enable()
+
+
+def plot_map(mode, model, syear, eyear, season, eof_Nth, frac_Nth, output_file_name, debug=False):
  """Plot dive down map and save
 
  Parameters
@@ -67,6 +74,8 @@ def plot_map(mode, model, syear, eyear, season, eof_Nth, frac_Nth, output_file_n
  + percentage
  )
 
+ debug_print('plot_map: projection, plot_title:' + projection +', '+ plot_title, debug)
+
  gridline = True
 
  if mode in [
@@ -97,6 +106,7 @@ def plot_map(mode, model, syear, eyear, season, eof_Nth, frac_Nth, output_file_n
  levels=levels,
  maskout=maskout,
  center_lon_global=center_lon_global,
+ debug=debug,
  )
 
 
@@ -138,6 +148,8 @@ def plot_map_cartopy(
  Switch for debugging print statements (default is False)
  """
 
+ debug_print('plot_map_cartopy starts', debug)
+
  lons = data.getLongitude()
  lats = data.getLatitude()
 
@@ -148,6 +160,8 @@ def plot_map_cartopy(
  if debug:
  print(min_lon, max_lon, min_lat, max_lat)
 
+ debug_print('Central longitude setup starts', debug)
+ debug_print('proj: '+proj, debug)
  # map types example:
  # https://github.com/SciTools/cartopy-tutorial/blob/master/tutorial/projections_crs_and_terms.ipynb
 
@@ -170,22 +184,35 @@ def plot_map_cartopy(
  central_latitude=central_latitude,
  standard_parallels=(20, max_lat),
  )
+ else:
+ print('Error: projection not defined!')
+
+ if debug:
+ debug_print('Central longitude setup completes', debug)
+ print('projection:', projection)
 
  # Generate plot
- fig = plt.figure(figsize=(8, 6))
- ax = plt.axes(projection=projection)
+ debug_print('Generate plot starts', debug)
+ #fig = plt.figure(figsize=(8, 6))
+ debug_print('fig done', debug)
+ #ax = plt.axes(projection=projection)
+ ax = plt.axes(projection=ccrs.NorthPolarStereo())
+ debug_print('ax done', debug)
  im = ax.contourf(
  lons,
  lats,
- data,
+ np.array(data),
  transform=ccrs.PlateCarree(),
  cmap=cmap,
  levels=levels,
  extend="both",
  )
+ debug_print('contourf done', debug)
  ax.coastlines()
+ debug_print('Generate plot completed', debug)
 
  # Grid Lines and tick labels
+ debug_print('projection starts', debug)
  if proj == "PlateCarree":
  if data_area == "global":
  if gridline:
@@ -203,6 +230,7 @@ def plot_map_cartopy(
  if gridline:
  gl = ax.gridlines(alpha=0.5, linestyle="--")
  elif "Stereo" in proj:
+ debug_print(proj + ' start', debug)
  if gridline:
  gl = ax.gridlines(draw_labels=True, alpha=0.5, linestyle="--")
  gl.xlocator = mticker.FixedLocator(
@@ -226,6 +254,7 @@ def plot_map_cartopy(
  verts = np.vstack([np.sin(theta), np.cos(theta)]).T
  circle = mpath.Path(verts * radius + center)
  ax.set_boundary(circle, transform=ax.transAxes)
+ debug_print(proj + ' plotted', debug)
  elif proj == "Lambert":
  # Make a boundary path in PlateCarree projection, I choose to start in
  # the bottom left and go round anticlockwise, creating a boundary point
@@ -250,6 +279,7 @@ def plot_map_cartopy(
  right_label = ea.get_position()[0] > 0
  if right_label:
  ea.set_visible(False)
+ debug_print('projection completed', debug)
 
  # Add title
  plt.title(title, pad=15, fontsize=15)
@@ -268,13 +298,14 @@ def plot_map_cartopy(
  if maskout is not None:
  if maskout == "land":
  ax.add_feature(
- cartopy.feature.LAND, zorder=100, edgecolor="k", facecolor="lightgrey"
+ cartopy_land, zorder=100, edgecolor="k", facecolor="lightgrey"
  )
  if maskout == "ocean":
  ax.add_feature(
- cartopy.feature.OCEAN, zorder=100, edgecolor="k", facecolor="lightgrey"
+ cartopy_ocean, zorder=100, edgecolor="k", facecolor="lightgrey"
  )
 
  # Done, save figure
+ debug_print('plot done, save figure as ' + filename, debug)
  fig.savefig(filename)
  plt.close("all")

pcmdi_metrics/variability_mode/scripts_pcmdi/parallel_submitter_all_modes_models.py

@@ -0,0 +1,131 @@
+import datetime
+import glob
+import os
+
+from pcmdi_metrics.misc.scripts import parallel_submitter
+from pcmdi_metrics.variability_mode.lib import sort_human
+
+
+def find_latest(path):
+ dir_list = [p for p in glob.glob(path + "/v????????")]
+ return sorted(dir_list)[-1]
+
+# ---------------------------------------------------------------
+
+mip = "cmip6"
+exp = "historical"
+
+datadir = '/p/user_pub/pmp/pmp_results/pmp_v1.1.2/additional_xmls/latest'
+param_dir = "../param"
+
+modes = ["NAM", "SAM", "NAO", "PNA", "NPO", "PDO", "NPGO"]
+
+case_id = "{:v%Y%m%d}".format(datetime.datetime.now())
+
+num_workers = 40
+debug = True 
+
+# ---------------------------------------------------------------
+
+cmds_list = list()
+logfilename_list = list()
+
+for mode in modes:
+ if mode in ["PDO", "NPGO", "AMO"]:
+ var = "ts"
+ param_file = "myParam_demo_PDO.py"
+ else:
+ var = "psl"
+ param_file = "myParam_demo_NAM.py"
+
+ if mode in ["SAM"]:
+ osyear = 1955
+ else:
+ osyear = 1900
+
+ if mode in ["NPO", "NPGO"]:
+ eofn_obs = 2
+ eofn_mod = 2
+ else:
+ eofn_obs = 1
+ eofn_mod = 1
+
+ datadir_ver = find_latest(datadir)
+ datadir_final = os.path.join(datadir_ver, mip, exp, 'atmos', 'mon', var)
+
+ xml_list = glob.glob(os.path.join(datadir_final, '*.xml'))
+
+ # get list of models
+ models_list = sort_human(
+ [r.split("/")[-1].split(".")[2] for r in xml_list]
+ )
+ # remove repeat
+ models_list = sort_human(list(dict.fromkeys(models_list)))
+
+ if debug:
+ models_list = models_list[0:3]
+
+ print(models_list)
+ print(len(models_list))
+
+ for model in models_list:
+ file_list_model = glob.glob(os.path.join(datadir_final, '*.' + model + '.*.xml'))
+ runs_list = sort_human(
+ [
+ r.split("/")[-1].split(".")[3]
+ for r in file_list_model
+ ]
+ )
+
+ if debug:
+ runs_list = runs_list[0:1]
+
+ print(model, runs_list)
+
+ for run in runs_list:
+
+ cmd_content = [
+ 'variability_modes_driver.py', 
+ '-p', os.path.join(param_dir, param_file), 
+ '--variability_mode', mode,
+ '--osyear', str(osyear),
+ '--eofn_obs', str(eofn_obs),
+ '--eofn_mod', str(eofn_mod),
+ '--mip', mip, '--exp', exp,
+ '--modnames', model,
+ '--realization', run,
+ ]
+
+ if model != models_list[0] or run != runs_list[0]:
+ cmd_content.extend(["--no_plot_obs", "--no_nc_out_obs"])
+
+ if debug:
+ cmd_content.append('--debug True')
+
+ cmd = ' '.join(cmd_content)
+
+ log_file = '_'.join(['variability_modes', mode, mip, exp, model, run]) + '.txt'
+
+ print(cmd) 
+ cmds_list.append(cmd)
+ logfilename_list.append(log_file)
+
+# =================================================
+# Run subprocesses in parallel
+# -------------------------------------------------
+# log dir
+log_dir = os.path.join(
+ "/p/user_pub/pmp/pmp_results/pmp_v1.1.2",
+ "log",
+ "variability_modes",
+ mip, exp, case_id)
+
+os.makedirs(log_dir, exist_ok=True)
+
+parallel_submitter(
+ cmds_list,
+ log_dir=log_dir,
+ logfilename_list=logfilename_list,
+ num_workers=num_workers,
+)
+