Add IODA comparisons

src/swell/suites/compare/flow.cylc

@@ -41,9 +41,11 @@
             {% for path in comparison_experiment_paths %}
             JediOopsLogParser-{{model_component}}-{{ loop.index0 }}
             {% endfor %}
-            JediLogComparison-{{model_component}}[^]:fail? => EvaComparisonIncrement-{{model_component}} => PublishComparisons => comparison_fail
-            JediLogComparison-{{model_component}}[^]:fail? => EvaComparisonJediLog-{{model_component}} => PublishComparisons => comparison_fail
-	    JediLogComparison-{{model_component}}[^]:fail? => EvaComparisonObservations-{{model_component}} => PublishComparisons => comparison_fail
+            CompareIodaObservations-{{model_component}}?
+            CompareIodaObservations-{{model_component}}:fail? | JediLogComparison-{{model_component}}[^]:fail? => EvaComparisonObservations-{{model_component}}? => PublishComparisons?
+            JediLogComparison-{{model_component}}[^]:fail? => EvaComparisonIncrement-{{model_component}} => PublishComparisons?
+            JediLogComparison-{{model_component}}[^]:fail? => EvaComparisonJediLog-{{model_component}} => PublishComparisons?
+            PublishComparisons? => comparison_fail?
         {% endif %}
         {% endfor %}
         """
@@ -75,6 +77,12 @@
     [[EvaComparisonJediLog-{{model_component}}]]
         script = "swell task EvaComparisonJediLog $config -d $datetime -m {{model_component}}"
 
+    [[PublishComparisons]]
+         script = "swell task PublishComparisons $config -d $datetime -m {{model_component}}"
+
+    [[CompareIodaObservations-{{model_component}}]]
+        script = "swell task CompareIodaObservations $config -d $datetime -m {{model_component}}"
+
     [[EvaComparisonObservations-{{model_component}}]]
         script = "swell task EvaComparisonObservations $config -d $datetime -m {{model_component}}"
         platform = {{platform}}
@@ -84,9 +92,6 @@
             --{{key}} = {{value}}
         {%- endfor %}
 
-    [[PublishComparisons]]
-         script = "swell task PublishComparisons $config -d $datetime -m {{model_component}}"
-
     {% if comparison_experiment_paths is mapping %}
     {% for path in comparison_experiment_paths.values() %}
     [[JediOopsLogParser-{{model_component}}-{{ loop.index0 }}]]

src/swell/suites/compare/suite_config.py

@@ -30,6 +30,9 @@ class SuiteConfig(QuestionContainer, Enum):
             qd.cycle_times(default_value=[None], widget_type=WidgetType.STRING_CHECK_LIST),
             qd.model_components(),
             qd.runahead_limit(),
+            qd.ioda_fields_for_comparison(['hofx0/{variable}',
+                                           'hofx1/{variable}']),
+            qd.observations([]),
         ]
     )
 
@@ -41,6 +44,21 @@ class SuiteConfig(QuestionContainer, Enum):
             compare,
             qd.comparison_log_type('variational'),
             qd.model_components(['geos_marine']),
+            qd.observations([
+                "adt_cryosat2n",
+                "adt_jason3",
+                "adt_saral",
+                "adt_sentinel3a",
+                "adt_sentinel3b",
+                "insitu_profile_argo",
+                "sst_ostia",
+                "sss_smos",
+                "sss_smapv5",
+                "sst_abi_g16_l3c",
+                "sst_gmi_l3u",
+                "sst_viirs_n20_l3u",
+                "temp_profile_xbt"
+            ])
         ]
     )
 
@@ -52,6 +70,42 @@ class SuiteConfig(QuestionContainer, Enum):
             compare,
             qd.comparison_log_type('variational'),
             qd.model_components(['geos_atmosphere']),
+            qd.observations([
+                "aircraft_temperature",
+                "aircraft_wind",
+                "airs_aqua",
+                "amsr2_gcom-w1",
+                "amsua_aqua",
+                "amsua_metop-b",
+                "amsua_metop-c",
+                "amsua_n15",
+                "amsua_n18",
+                "amsua_n19",
+                "atms_n20",
+                "atms_npp",
+                "avhrr3_metop-b",
+                "avhrr3_n18",
+                "avhrr3_n19",
+                "cris-fsr_n20",
+                "cris-fsr_npp",
+                "gmi_gpm",
+                "gps",
+                "iasi_metop-b",
+                "iasi_metop-c",
+                "mhs_metop-b",
+                "mhs_metop-c",
+                "mhs_n19",
+                "mls55_aura",
+                "omi_aura",
+                "ompsnm_npp",
+                "pibal",
+                "satwind",
+                "scatwind",
+                "sfcship",
+                "sfc",
+                "sondes",
+                "ssmis_f17"
+            ])
         ]
     )
 
@@ -63,6 +117,10 @@ class SuiteConfig(QuestionContainer, Enum):
             compare,
             qd.comparison_log_type('variational'),
             qd.model_components(['geos_cf']),
+            qd.observations([
+                "tempo_no2_tropo",
+                "tropomi_s5p_no2_tropo",
+            ])
         ]
     )
 
@@ -74,6 +132,95 @@ class SuiteConfig(QuestionContainer, Enum):
             compare,
             qd.comparison_log_type('fgat'),
             qd.model_components(['geos_marine']),
+            qd.observations([
+                "adt_cryosat2n",
+                "adt_jason3",
+                "adt_jason3n",
+                "adt_saral",
+                "adt_sentinel3a",
+                "adt_sentinel3b",
+                "adt_sentinel6a",
+                "adt_swot_nadir",
+                "insitu_profile_argo",
+                "insitu_profile_ctd",
+                "insitu_profile_pirata",
+                "insitu_profile_rama",
+                "insitu_profile_tao",
+                "icec_amsr2_north",
+                "icec_amsr2_south",
+                "icec_nsidc_nh",
+                "icec_nsidc_sh",
+                "sst_ostia",
+                "sss_smos",
+                "sss_smapv5",
+                "sst_abi_g16_l3c",
+                "sst_avhrrf_mb_l3u",
+                "sst_avhrrf_mc_l3u",
+                "sst_viirs_n20_l3u",
+                "sst_viirs_npp_l3u",
+                "temp_profile_xbt"
+            ])
+        ]
+    )
+
+    # --------------------------------------------------------------------------------------------------
+
+    compare_hofx = QuestionList(
+        list_name="compare_hofx",
+        questions=[
+            compare,
+            qd.comparison_log_type('hofx'),
+            qd.model_components(['geos_atmosphere']),
+            qd.ioda_fields_for_comparison(['hofx/{variable}']),
+            qd.observations([
+                "aircraft_temperature",
+                "aircraft_wind",
+                "airs_aqua",
+                "amsr2_gcom-w1",
+                "amsua_aqua",
+                "amsua_metop-b",
+                "amsua_metop-c",
+                "amsua_n15",
+                "amsua_n18",
+                "amsua_n19",
+                "atms_n20",
+                "atms_npp",
+                "avhrr3_metop-b",
+                "avhrr3_n18",
+                "avhrr3_n19",
+                "cris-fsr_n20",
+                "cris-fsr_npp",
+                "gmi_gpm",
+                "gps",
+                "iasi_metop-b",
+                "iasi_metop-c",
+                "mhs_metop-b",
+                "mhs_metop-c",
+                "mhs_n19",
+                "mls55_aura",
+                "omi_aura",
+                "ompsnm_npp",
+                "pibal",
+                "satwind",
+                "scatwind",
+                "sfcship",
+                "sfc",
+                "sondes",
+                "ssmis_f17"
+            ]),
+            qd.mock_experiment(True)
+        ]
+    )
+
+    # --------------------------------------------------------------------------------------------------
+
+    compare_hofx_cf = QuestionList(
+        list_name="compare_hofx_cf",
+        questions=[
+            compare,
+            qd.comparison_log_type('hofx'),
+            qd.ioda_fields_for_comparison(['hofx/{variable}']),
+            qd.model_components(['geos_cf']),
         ]
     )
 

src/swell/tasks/compare_ioda_observations.py

@@ -0,0 +1,146 @@
+# (C) Copyright 2021- United States Government as represented by the Administrator of the
+# National Aeronautics and Space Administration. All Rights Reserved.
+
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+
+# --------------------------------------------------------------------------------------------------
+
+import os
+import numpy as np
+from ruamel.yaml import YAML
+from netCDF4 import Dataset
+from pathlib import Path
+from swell.tasks.base.task_base import taskBase
+from swell.utilities.comparisons import comparison_tags
+
+# --------------------------------------------------------------------------------------------------
+
+
+class CompareIodaObservations(taskBase):
+
+    def ioda_means(self, exp_path: str, observation: str, field: str, cutoff: int | None = None):
+        ''' Reads observation files for ioda variables, returning total length of
+        the array and mean of all data points.
+
+        Arguments:
+        exp_path: experiment.yaml location for experiment to compare
+        observation: ioda name of observation
+        field: name of data field
+        cutoff: int array index to cut off mean of array
+
+        Returns:
+        Dictionary of simulated variable with length and mean
+        '''
+
+        exp_path = Path(exp_path)
+
+        window_length = self.config.window_length()
+        window_begin = self.da_window_params.window_begin(window_length)
+
+        obs_path = (Path(os.path.dirname(exp_path)) / '..' / 'run' /
+                    self.__datetime__.string_directory() / self.get_model())
+
+        # Get list of obs files
+        obs_files = list(obs_path.glob(f'*.{observation}.{window_begin}.nc4'))
+
+        # If empty obs return black dictionary
+        if len(obs_files) > 0:
+            obs_file = obs_files[0]
+        else:
+            return {}
+
+        # Read obs config yaml to get names of simulated variables
+        yaml = YAML(typ='safe')
+        with open(obs_path / 'obs.yaml', 'r') as f:
+            obs_config = yaml.load(f)
+
+        for ob in obs_config:
+            if ob['observation_name'] == observation:
+                simulated_variables = ob['obs space']['simulated variables']
+
+        field_means = {}
+        with Dataset(obs_file, 'r') as ds:
+            for sim_var in simulated_variables:
+                var_name = field.format(variable=sim_var)
+                obs_var = ds[var_name]
+                field_means[var_name] = {}
+                field_means[var_name]['length'] = len(obs_var)
+                field_means[var_name]['mean'] = np.mean(obs_var[0:cutoff])
+
+        return field_means
+
+    # --------------------------------------------------------------------------------------------------
+
+    def execute(self) -> None:
+
+        '''
+        Reads observation files for ioda variables, compare total length of array and mean
+        of data points to evaluate diff. Output is sent to
+        <cycle_dir>/ioda_<observation>_comparison.txt
+        '''
+
+        comparison_experiment_paths = self.config.comparison_experiment_paths()
+        observations = self.config.observations()
+        ioda_fields = self.config.ioda_fields_for_comparison()
+
+        experiment_tag_paths = comparison_tags(comparison_experiment_paths, self.logger)
+
+        tag_1 = list(experiment_tag_paths.keys())[0]
+        tag_2 = list(experiment_tag_paths.keys())[1]
+
+        path_1 = list(experiment_tag_paths.values())[0]
+        path_2 = list(experiment_tag_paths.values())[1]
+
+        for observation in observations:
+
+            output_str = f'{observation} Comparison Results\n'
+            output_str += f'{tag_1}: {path_1}\n'
+            output_str += f'{tag_2}: {path_2}\n'
+            output_str += '\n'
+            passed = True
+
+            for field in ioda_fields:
+
+                field_means_1 = self.ioda_means(path_1, observation, field, cutoff=int(1e4))
+                field_means_2 = self.ioda_means(path_2, observation, field, cutoff=int(1e4))
+
+                if len(field_means_1) != len(field_means_2):
+                    raise Exception(f'Length of {field} fields does not '
+                                    'match between experiments.')
+
+                output_str += f'{observation}\n'
+                for sim_var in field_means_1.keys():
+                    len_1 = field_means_1[sim_var]['length']
+                    len_2 = field_means_2[sim_var]['length']
+
+                    mean_1 = field_means_1[sim_var]['mean']
+                    mean_2 = field_means_2[sim_var]['mean']
+
+                    output_str += f'{sim_var}\n'
+                    if len_1 != len_2 or mean_1 != mean_2:
+                        tag_length = max(len(tag_1), len(tag_2)) + 2
+                        len_length = max(len(str(len_1)), len(str(len_2))) + 2
+                        mean_length = max(len(str(mean_1)), len(str(mean_2))) + 2
+                        output_str += (f'{"":<{tag_length}} {"Length":<{len_length}} '
+                                       f'{"Mean":<{mean_length}}\n')
+                        output_str += (f'{tag_1:<{tag_length}} {len_1:<{len_length}} '
+                                       f'{mean_1:<{mean_length}}\n')
+                        output_str += (f'{tag_2:<{tag_length}} {len_2:<{len_length}} '
+                                       f'{mean_2:<{mean_length}}\n\n')
+                        passed = False
+                    else:
+                        output_str += f'Passed\n\n'
+
+                # Fail suite if not passed
+                if not passed:
+                    output_file = Path(self.cycle_dir()) / f'ioda_{observation}_comparison.txt'
+
+                    # Output to file
+                    with open(output_file, 'w') as f:
+                        f.write(output_str)
+                    raise Exception(f'Mismatch in IODA field length or average, '
+                                    f'check {output_file}')
+
+# --------------------------------------------------------------------------------------------------