Add stratosphere and troposphere ozone (#481)

* Add_column_Ozone

* address reviews

* address review 2

* bring arrows on colorbars back

acme_diags/derivations/acme.py

@@ -1351,6 +1351,8 @@ def cosp_histogram_standardize(cld: "FileVariable"):
             (("tasmax",), lambda t: convert_units(t, target_units="DegC")),
         ]
     ),
+    "TCO": OrderedDict([(("TCO",), rename)]),
+    "SCO": OrderedDict([(("SCO",), rename)]),
     # Land variables
     "SOILWATER_10CM": OrderedDict([(("mrsos",), rename)]),
     "SOILWATER_SUM": OrderedDict([(("mrso",), rename)]),

acme_diags/derivations/default_regions.py

@@ -14,6 +14,7 @@
     "20N50N": {"domain": cdutil.region.domain(latitude=(20.0, 50, "ccb"))},
     "50N90N": {"domain": cdutil.region.domain(latitude=(50.0, 90, "ccb"))},
     "60S90N": {"domain": cdutil.region.domain(latitude=(-60.0, 90, "ccb"))},
+    "60S60N": {"domain": cdutil.region.domain(latitude=(-60.0, 60, "ccb"))},
     "ocean": {
         "value": 0.65,
     },

acme_diags/driver/annual_cycle_zonal_mean_driver.py

@@ -11,6 +11,7 @@
     get_name_and_yrs,
     regrid_to_lower_res,
     save_ncfiles,
+    select_region_lat_lon,
 )
 from acme_diags.plot import plot
 
@@ -50,15 +51,51 @@ def run_diag(parameter: "CoreParameter"):
         )
 
         test_ac_zonal_mean = cdutil.averager(test_ac, axis="x", weights="generate")
-        ref_ac_zonal_mean = cdutil.averager(ref_ac, axis="x", weights="generate")
         test_ac_reg_zonal_mean = cdutil.averager(
             test_ac_reg, axis="x", weights="generate"
         )
-        ref_ac_reg_zonal_mean = cdutil.averager(
-            ref_ac_reg, axis="x", weights="generate"
-        )
+
+        if (
+            parameter.ref_name == "OMI-MLS"
+        ):  # SCO from OMI-MLS only available as (time, lat)
+            test_ac_reg_zonal_mean = select_region_lat_lon(
+                "60S60N", test_ac_reg_zonal_mean, parameter
+            )
+            test_ac_zonal_mean = select_region_lat_lon(
+                "60S60N", test_ac_zonal_mean, parameter
+            )
+            if var == "SCO":
+                ref_ac_zonal_mean = ref_ac
+                ref_ac_reg_zonal_mean = ref_ac_reg
+            else:
+                ref_ac_zonal_mean = cdutil.averager(
+                    ref_ac, axis="x", weights="generate"
+                )
+                ref_ac_reg_zonal_mean = cdutil.averager(
+                    ref_ac_reg, axis="x", weights="generate"
+                )
+
+        else:
+            ref_ac_zonal_mean = cdutil.averager(ref_ac, axis="x", weights="generate")
+            ref_ac_reg_zonal_mean = cdutil.averager(
+                ref_ac_reg, axis="x", weights="generate"
+            )
+
+        # if var == 'SCO' and parameter.ref_name=='OMI-MLS':  # SCO from OMI-MLS only available as (time, lat)
+        #    ref_ac_zonal_mean = ref_ac
+        #    ref_ac_reg_zonal_mean = ref_ac_reg
+
+        #    test_ac_reg_zonal_mean = select_region_lat_lon("60S60N", test_ac_reg_zonal_mean, parameter)
+        #    test_ac_zonal_mean = select_region_lat_lon("60S60N", test_ac_zonal_mean, parameter)
+        # else:
+        #    ref_ac_zonal_mean = cdutil.averager(ref_ac, axis="x", weights="generate")
+        #    ref_ac_reg_zonal_mean = cdutil.averager(
+        #        ref_ac_reg, axis="x", weights="generate"
+        #    )
 
         diff_ac = test_ac_reg_zonal_mean - ref_ac_reg_zonal_mean
+        diff_ac.setAxis(1, test_ac_reg_zonal_mean.getAxis(1))
+        diff_ac.setAxis(0, test_ac_reg_zonal_mean.getAxis(0))
 
         parameter.var_id = var
         parameter.output_file = "-".join([ref_name, var, "Annual-Cycle"])
@@ -117,8 +154,10 @@ def _create_annual_cycle(dataset: Dataset, variable: str) -> "TransientVariable"
 
             var_ann_cycle.setAxis(0, time)
             time.designateTime()
-            var_ann_cycle.setAxis(1, var.getAxis(0))
-            var_ann_cycle.setAxis(2, var.getAxis(1))
+
+            for iax in list(range(len(var.shape))):
+                var_ann_cycle.setAxis(1 + iax, var.getAxis(iax))
+            # var_ann_cycle.setAxis(2, var.getAxis(1))
 
             var_ann_cycle[0] = var
         else:

acme_diags/driver/default_diags/annual_cycle_zonal_mean_model_vs_model.cfg

@@ -269,3 +269,24 @@ reference_colormap = "Purples"
 diff_colormap = "RdBu"
 contour_levels = [0., 0.025, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5]
 diff_levels = [-0.1, -0.08, -0.06, -0.05, -0.04, -0.02, -0.01,  0.01, 0.02, 0.04, 0.05, 0.06, 0.08, 0.1]
+
+[#]
+sets = ["annual_cycle_zonal_mean"]
+case_id = "model_vs_model"
+variables = ["SCO"]
+test_colormap = "WhiteBlueGreenYellowRed.rgb"
+reference_colormap = "WhiteBlueGreenYellowRed.rgb"
+diff_colormap = "diverging_bwr.rgb"
+contour_levels = [180,200,220,240,260,280,300,320,340,360,380]
+diff_levels = [-50, -40, -30, -20, -10, -5,-2,2, 5, 10, 20, 30,40,50]
+
+
+[#]
+sets = ["annual_cycle_zonal_mean"]
+case_id = "model_vs_model"
+variables = ["TCO"]
+test_colormap = "WhiteBlueGreenYellowRed.rgb"
+reference_colormap = "WhiteBlueGreenYellowRed.rgb"
+diff_colormap = "diverging_bwr.rgb"
+contour_levels = [12,16,20,24,28,32,36,40,44]
+diff_levels = [-20,-15,-10,-5,-2,2,5,10,15,20]

acme_diags/driver/default_diags/annual_cycle_zonal_mean_model_vs_obs.cfg

@@ -437,3 +437,27 @@ reference_colormap = "RdBu"
 diff_colormap = "RdBu"
 contour_levels = [-16, -11, -7, -4,-2,-1,1,2,4,7,11,16]
 diff_levels = [-10, -7, -4, -2, -1,-0.5,0.5, 1, 2, 4, 7, 10]
+
+[#]
+sets = ["annual_cycle_zonal_mean"]
+case_id = "OMI-MLS"
+variables = ["SCO"]
+ref_name = "OMI-MLS"
+reference_name = "OMI-MLS"
+test_colormap = "WhiteBlueGreenYellowRed.rgb"
+reference_colormap = "WhiteBlueGreenYellowRed.rgb"
+diff_colormap = "diverging_bwr.rgb"
+contour_levels = [180,200,220,240,260,280,300,320,340,360,380]
+diff_levels = [ -40, -30, -20, -15, -10, -5,-2,2, 5, 10, 15, 20, 30,40]
+
+[#]
+sets = ["annual_cycle_zonal_mean"]
+case_id = "OMI-MLS"
+variables = ["TCO"]
+ref_name = "OMI-MLS"
+reference_name = "OMI-MLS"
+test_colormap = "WhiteBlueGreenYellowRed.rgb"
+reference_colormap = "WhiteBlueGreenYellowRed.rgb"
+diff_colormap = "diverging_bwr.rgb"
+contour_levels = [12,16,20,24,28,32,36,40,44]
+diff_levels = [-20,-15,-10,-5,-2,2,5,10,15,20]

acme_diags/driver/default_diags/lat_lon_model_vs_model.cfg

@@ -633,3 +633,14 @@ reference_colormap = "Purples"
 diff_colormap = "RdBu"
 contour_levels = [0., 0.025, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5]
 diff_levels = [-0.1, -0.08, -0.06, -0.05, -0.04, -0.02, -0.01,  0.01, 0.02, 0.04, 0.05, 0.06, 0.08, 0.1]
+
+[#]
+sets = ["lat_lon"]
+case_id = "model_vs_model"
+variables = ["TCO"]
+seasons = ["ANN", "DJF", "MAM", "JJA", "SON"]
+test_colormap = "WhiteBlueGreenYellowRed.rgb"
+reference_colormap = "WhiteBlueGreenYellowRed.rgb"
+diff_colormap = "diverging_bwr.rgb"
+contour_levels = [12,16,20,24,28,32,36,40,44]
+diff_levels = [-20,-15,-10,-5,-2,2,5,10,15,20]

acme_diags/driver/default_diags/lat_lon_model_vs_obs.cfg

@@ -1316,3 +1316,16 @@ reference_colormap = "RdBu"
 diff_colormap = "RdBu"
 contour_levels = [-16, -11, -7, -4,-2,-1,1,2,4,7,11,16]
 diff_levels = [-10, -7, -4, -2, -1,-0.5,0.5, 1, 2, 4, 7, 10]
+
+[#]
+sets = ["lat_lon"]
+case_id = "OMI-MLS"
+variables = ["TCO"]
+ref_name = "OMI-MLS"
+reference_name = "OMI-MLS"
+seasons = ["ANN", "DJF", "MAM", "JJA", "SON"]
+test_colormap = "WhiteBlueGreenYellowRed.rgb"
+reference_colormap = "WhiteBlueGreenYellowRed.rgb"
+diff_colormap = "diverging_bwr.rgb"
+contour_levels = [12,16,20,24,28,32,36,40,44]
+diff_levels = [-20,-15,-10,-5,-2,2,5,10,15,20]

acme_diags/driver/utils/dataset.py

@@ -313,10 +313,11 @@ def _find_climo_file(self, path_name, data_name, season):
         for filename in dir_files:
             if filename.startswith(data_name + "_" + season):
                 return os.path.join(path_name, filename)
-        # The below is only ran on model data, because a shorter name is passed into this software.
+        # The below is only ran on model data, because a shorter name is passed into this software. Won't work when use month name such as '01' as season.
         for filename in dir_files:
-            if filename.startswith(data_name) and season in filename:
-                return os.path.join(path_name, filename)
+            if season in ["ANN", "DJF", "MAM", "JJA", "SON"]:
+                if filename.startswith(data_name) and season in filename:
+                    return os.path.join(path_name, filename)
         # No file found.
         return ""
 

acme_diags/driver/utils/general.py

@@ -170,6 +170,22 @@ def pressure_to_plevs(var, plev):
     return var_p
 
 
+def select_region_lat_lon(region, var, parameter):
+    """Select desired regions from transient variables (no mask)."""
+    try:
+        # if region.find('global') == -1:
+        domain = regions_specs[region]["domain"]  # type: ignore
+        # print('Domain: ', domain)
+    except Exception:
+        pass
+        # print("No domain selector.")
+
+    var_selected = var(domain)
+    var_selected.units = var.units
+
+    return var_selected
+
+
 def select_region(region, var, land_frac, ocean_frac, parameter):
     """Select desired regions from transient variables."""
     domain = None

acme_diags/plot/cartopy/annual_cycle_zonal_mean_plot.py

@@ -53,14 +53,15 @@ def plot_panel(n, fig, var, clevels, cmap, title, parameters, stats=None):
     # Contour plot
     ax = fig.add_axes(panel[n])
     cmap = get_colormap(cmap, parameters)
-    p1 = ax.contourf(
+    # p1 = ax.contourf(
+    p1 = ax.pcolor(
         mon,
         lat,
         var,
         norm=norm,
-        levels=levels,
         cmap=cmap,
-        extend="both",
+        edgecolors="face",
+        shading="auto",
     )
 
     if title[0] is not None:
@@ -80,7 +81,7 @@ def plot_panel(n, fig, var, clevels, cmap, title, parameters, stats=None):
 
     # Color bar
     cbax = fig.add_axes((panel[n][0] + 0.6635, panel[n][1] + 0.0215, 0.0326, 0.1792))
-    cbar = fig.colorbar(p1, cax=cbax)
+    cbar = fig.colorbar(p1, cax=cbax, extend="both")
     w, h = get_ax_size(fig, cbax)
 
     if levels is None:
@@ -99,6 +100,9 @@ def plot_panel(n, fig, var, clevels, cmap, title, parameters, stats=None):
             pad = 30
         cbar.set_ticks(levels[1:-1])
         labels = [fmt % level for level in levels[1:-1]]
+        if all(x[-2:] == ".0" for x in labels):
+            labels = [x[:-2] for x in labels]
+            pad = pad - 5
         cbar.ax.set_yticklabels(labels, ha="right")
         cbar.ax.tick_params(labelsize=9.0, pad=pad, length=0)
 

analysis_data_preprocess/create_OMI-MLS_climo.sh

@@ -0,0 +1,47 @@
+#! /usr/bin/env bash
+#Tropospheric column ozone (TCO) data is downloaded from https://acd-ext.gsfc.nasa.gov/Data_services/cloud_slice/new_data.html
+#stratospheric column ozone (SCO) data is provided by J. R. Ziemke based on https://doi.org/10.5194/acp-19-3257-2019
+#Both data are originally in ASCII format and converted into netCDF files by Qi Tang.
+
+path='/Users/zhang40/Documents/ACME_simulations/O3/'
+
+original_data_path=$path'original_data/'
+time_series_output_path=$path'time_series/'
+climo_output_path=$path'climatology/'
+tmp=$path'tmp/'
+
+#start_yr=1997
+#end_yr=2016
+
+mkdir $time_series_output_path
+mkdir $climo_output_path
+mkdir $tmp
+
+
+#Add lon dimention to zonel mean
+ncap2 -s 'SCO=O3strat' ${original_data_path}O3strat_ZMK.nc ${time_series_output_path}SCO_200410_201712.nc
+#cp ${original_data_path}O3strat_ZMK.nc ${time_series_output_path}SCO_200410_201712.nc
+cdo splityear ${time_series_output_path}SCO_200410_201712.nc ${tmp}sco
+cp ${original_data_path}tropO3clmn.nc ${time_series_output_path}TCO_200410_201712.nc
+cdo splityear ${original_data_path}tropO3clmn.nc ${tmp}tco
+
+
+cd $original_data_path
+for yr in {2005..2017}; do      # Loop over years
+    yyyy=`printf "%04d" $yr`
+    echo $yyyy
+
+    for mth in {1..12}; do
+        mm=`printf "%02d" $mth`
+        ncks -O -F -d time,${mth} ${tmp}sco${yyyy}.nc ${tmp}sco_${yyyy}${mm}.nc
+        ncks -O -F -d time,${mth} ${tmp}tco${yyyy}.nc ${tmp}OMI-MLS_${yyyy}${mm}.nc
+        ncks -A ${tmp}sco_${yyyy}${mm}.nc ${tmp}OMI-MLS_${yyyy}${mm}.nc
+        done
+done
+ncclimo -a sdd -c OMI-MLS_200501.nc -s 2005 -e 2017 -i ${tmp} -o ${climo_output_path}
+
+
+exit
+
+# Concatenate monthly files together
+#ncrcat -O in_??????.nc out.nc

tests/integration/test_run.py

@@ -84,14 +84,15 @@ def test_all_sets_and_all_seasons(self):
                 msg = "Count for {} is invalid: {}"
                 self.fail(msg.format(set_name, count))
 
-        all_season_counts = list(season_counter.values())
         # enso_diags and streamflow only run ANN, no seasons
         # So, reduce the ANN count by the number of times these appear
-        all_season_counts[0] -= set_counter["enso_diags"]
-        all_season_counts[0] -= set_counter["streamflow"]
-        if not all(all_season_counts[0] == count for count in all_season_counts):
+        season_counter["ANN"] -= set_counter["enso_diags"]
+        season_counter["ANN"] -= set_counter["streamflow"]
+        if not all(season_counter["ANN"] == count for count in season_counter.values()):
             self.fail(
-                "Counts for the seasons don't match: {}".format(all_season_counts)
+                "In .cfg files, at least one season does not match the count for ANN: {}".format(
+                    season_counter
+                )
             )
 
     def test_zonal_mean_2d(self):