update

Author: lee1043

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,
+)
+