Merge pull request #1112 from PCMDI/zhan391

Code changes to implement PSA1 and PSA2 mode varibility analysis

Author: lee1043

doc/jupyter/Demo/basic_mov_param.py.in

@@ -23,7 +23,6 @@ varModel = 'psl'
 ModUnitsAdjust = (True, 'divide', 100.0)  # Pa to hPa
 msyear = 1900
 meyear = 2005
-eofn_mod = 1
 
 # OBSERVATIONS SETTINGS
 reference_data_path = '$INPUT_DIR$/obs4MIPs_PCMDI_monthly/NOAA-ESRL-PSD/20CR/mon/psl/gn/v20210727/psl_mon_20CR_PCMDI_gn_187101-201212.nc'
@@ -32,7 +31,6 @@ varOBS = 'psl'
 ObsUnitsAdjust = (True, 'divide', 100.0)  # Pa to hPa; or (False, 0, 0)
 osyear = 1900
 oeyear = 2005
-eofn_obs = 1
 
 # DIRECTORY WHERE TO PUT RESULTS
 results_dir = os.path.join(

doc/jupyter/Demo/basic_mov_param_sst.py.in

@@ -25,7 +25,6 @@ varModel = 'ts'
 ModUnitsAdjust = (True, "subtract", 273.15)  # degK to degC
 msyear = 1900
 meyear = 2005
-eofn_mod = 1
 
 # OBSERVATIONS SETTINGS
 reference_data_path = '$INPUT_DIR$/obs4MIPs_PCMDI_monthly/MOHC/HadISST-1-1/mon/ts/gn/v20210727/ts_mon_HadISST-1-1_PCMDI_gn_187001-201907.nc'
@@ -34,7 +33,6 @@ varOBS = 'ts'
 ObsUnitsAdjust = (True, "subtract", 273.15)  # degK to degC
 osyear = 1900
 oeyear = 2005
-eofn_obs = 1
 
 # DIRECTORY WHERE TO PUT RESULTS
 results_dir = os.path.join(

pcmdi_metrics/io/regions.py

@@ -39,6 +39,8 @@ def load_regions_specs() -> dict:
         "NAM": {"domain": {"latitude": (20.0, 90), "longitude": (-180, 180)}},
         "NAO": {"domain": {"latitude": (20.0, 80), "longitude": (-90, 40)}},
         "SAM": {"domain": {"latitude": (-20.0, -90), "longitude": (0, 360)}},
+        "PSA1": {"domain": {"latitude": (-20.0, -90), "longitude": (0, 360)}},
+        "PSA2": {"domain": {"latitude": (-20.0, -90), "longitude": (0, 360)}},
         "PNA": {"domain": {"latitude": (20.0, 85), "longitude": (120, 240)}},
         "NPO": {"domain": {"latitude": (20.0, 85), "longitude": (120, 240)}},
         "PDO": {"domain": {"latitude": (20.0, 70), "longitude": (110, 260)}},

pcmdi_metrics/variability_mode/lib/argparse_functions.py

@@ -33,6 +33,8 @@ def AddParserArgument(P):
         "- NAM: Northern Annular Mode\n"
         "- NAO: Northern Atlantic Oscillation\n"
         "- SAM: Southern Annular Mode\n"
+        "- PSA1: Pacific–South American pattern 1\n"
+        "- PSA2: Pacific–South American pattern 2\n"
         "- PNA: Pacific North American Pattern\n"
         "- PDO: Pacific Decadal Oscillation\n"
         "- NPO: North Pacific Oscillation\n"
@@ -224,6 +226,8 @@ def VariabilityModeCheck(mode, P):
             "NAM",
             "NAO",
             "SAM",
+            "PSA1",
+            "PSA2",
             "PNA",
             "PDO",
             "NPO",

pcmdi_metrics/variability_mode/lib/eof_analysis.py

@@ -155,6 +155,22 @@ def arbitrary_checking(mode, eof_Nth):
     elif mode == "SAM":
         if eof_Nth.sel({lat_key: slice(-60, -90)}).mean().item() >= 0:
             reverse_sign = True
+    elif mode == "PSA1":
+        if (
+            eof_Nth.sel({lat_key: slice(-59.5, -64.5), lon_key: slice(207.5, 212.5)})
+            .mean()
+            .item()
+            >= 0
+        ):
+            reverse_sign = True
+    elif mode == "PSA2":
+        if (
+            eof_Nth.sel({lat_key: slice(-57.5, -62.5), lon_key: slice(277.5, 282.5)})
+            .mean()
+            .item()
+            >= 0
+        ):
+            reverse_sign = True
     else:  # Minimum sign control part was left behind for any future usage..
         if not np.isnan(eof_Nth[-1, -1].item()):
             if eof_Nth[-1, -1].item() >= 0:

pcmdi_metrics/variability_mode/lib/lib_variability_mode.py

@@ -13,7 +13,7 @@
 
 import pcmdi_metrics
 from pcmdi_metrics.io import get_time, select_subset, xcdat_open
-from pcmdi_metrics.utils import apply_landmask
+from pcmdi_metrics.utils import apply_landmask, check_monthly_time_axis
 
 
 def tree():
@@ -67,6 +67,9 @@ def read_data_in(
     # Open data file
     ds = xcdat_open(path)
 
+    # Data QC check -- time axis check
+    check_monthly_time_axis(ds)
+
     # Time subset
     ds_time_subsetted = subset_time(ds, syear, eyear, debug=debug)
     data_timeseries = ds_time_subsetted[var_in_data]

pcmdi_metrics/variability_mode/lib/plot_map.py

@@ -57,7 +57,7 @@ def plot_map(
         projection = "Lambert"
     elif mode in ["NAM"]:
         projection = "Stereo_north"
-    elif mode in ["SAM"]:
+    elif mode in ["SAM", "PSA1", "PSA2"]:
         projection = "Stereo_south"
     else:
         sys.exit("Projection for " + mode + "is not defined.")

pcmdi_metrics/variability_mode/variability_modes_driver.py

@@ -15,6 +15,10 @@
 ## EOF2 based variability modes
 - NPO: North Pacific Oscillation (2nd EOFs of PNA domain)
 - NPGO: North Pacific Gyre Oscillation (2nd EOFs of PDO domain)
+- PSA2: Pacific South America Mode (2nd EOFs of SAM domain)
+
+## EOF3 based variability modes
+- PSA3: Pacific South America Mode (3rd EOFs of SAM domain)
 
 ## Reference:
 Lee, J., K. Sperber, P. Gleckler, C. Bonfils, and K. Taylor, 2019:
@@ -175,25 +179,34 @@
 eofn_obs = param.eofn_obs
 eofn_mod = param.eofn_mod
 
+if mode in ["NAM", "NAO", "SAM", "PNA", "PDO", "AMO"]:
+    eofn_expected = 1
+elif mode in ["NPGO", "NPO", "PSA1"]:
+    eofn_expected = 2
+elif mode in ["PSA2"]:
+    eofn_expected = 3
+else:
+    raise ValueError(
+        f"Mode '{mode}' is not defiend with associated expected EOF number"
+    )
+
 if eofn_obs is None:
-    if mode in ["NAM", "NAO", "SAM", "PNA", "PDO", "AMO"]:
-        eofn_obs = 1
-    elif mode in ["NPGO", "NPO"]:
-        eofn_obs = 2
-    else:
-        raise ValueError(f"{eofn_obs} is not given for {mode}")
+    eofn_obs = eofn_expected
 else:
     eofn_obs = int(eofn_obs)
+    if eofn_obs != eofn_expected:
+        raise ValueError(
+            f"Observation EOF number ({eofn_obs}) does not match expected EOF number ({eofn_expected}) for mode {mode}"
+        )
 
 if eofn_mod is None:
-    if mode in ["NAM", "NAO", "SAM", "PNA", "PDO", "AMO"]:
-        eofn_mod = 1
-    elif mode in ["NPGO", "NPO"]:
-        eofn_mod = 2
-    else:
-        raise ValueError(f"{eofn_mod} is not given for {mode}")
+    eofn_mod = eofn_expected
 else:
     eofn_mod = int(eofn_mod)
+    if eofn_mod != eofn_expected:
+        raise ValueError(
+            f"Model EOF number ({eofn_mod}) does not match expected EOF number ({eofn_expected}) for mode {mode}"
+        )
 
 print("eofn_obs:", eofn_obs)
 print("eofn_mod:", eofn_mod)
@@ -830,7 +843,7 @@
                 # Conventional EOF approach as supplementary
                 # - - - - - - - - - - - - - - - - - - - - - - - - -
                 if ConvEOF:
-                    eofn_mod_max = 3
+                    eofn_mod_max = max(3, eofn_mod)
 
                     # EOF analysis
                     debug_print("conventional EOF analysis start", debug)
@@ -976,15 +989,15 @@
                                 debug=debug,
                             )
                             plot_map(
-                                mode + "_teleconnection",
+                                f"{mode}_teleconnection",
                                 f"{mip.upper()} {model} ({run}) - EOF{n + 1}",
                                 msyear,
                                 meyear,
                                 season,
                                 # eof_lr(longitude=(lon1_global, lon2_global)),
                                 eof_lr,
                                 frac,
-                                output_img_file + "_teleconnection",
+                                f"{output_img_file}_teleconnection",
                                 debug=debug,
                             )