Merge pull request #978 from PCMDI/mean_clim_3d_field_vertical_coord_…

…unit_check

Mean climate 3d field vertical coordinate QC

pcmdi_metrics/io/default_regions_define.py

@@ -5,25 +5,29 @@ def load_regions_specs():
 
  regions_specs = {
  # Mean Climate
+ "global": {},
  "NHEX": {"domain": {"latitude": (30.0, 90)}},
  "SHEX": {"domain": {"latitude": (-90.0, -30)}},
  "TROPICS": {"domain": {"latitude": (-30.0, 30)}},
- "global": {},
  "90S50S": {"domain": {"latitude": (-90.0, -50)}},
  "50S20S": {"domain": {"latitude": (-50.0, -20)}},
  "20S20N": {"domain": {"latitude": (-20.0, 20)}},
  "20N50N": {"domain": {"latitude": (20.0, 50)}},
  "50N90N": {"domain": {"latitude": (50.0, 90)}},
+ "CONUS": {"domain": {"latitude": (24.7, 49.4), "longitude": (-124.78, -66.92)}},
+ "land": {"value": 100},
  "land_NHEX": {"value": 100, "domain": {"latitude": (30.0, 90)}},
  "land_SHEX": {"value": 100, "domain": {"latitude": (-90.0, -30)}},
  "land_TROPICS": {"value": 100, "domain": {"latitude": (-30.0, 30)}},
- "land": {"value": 100},
+ "land_CONUS": {"value": 100, "domain": {"latitude": (24.7, 49.4), "longitude": (-124.78, -66.92)}},
+ "ocean": {"value": 0},
  "ocean_NHEX": {"value": 0, "domain": {"latitude": (30.0, 90)}},
  "ocean_SHEX": {"value": 0, "domain": {"latitude": (-90.0, -30)}},
  "ocean_TROPICS": {"value": 0, "domain": {"latitude": (30.0, 30)}},
- "ocean": {"value": 0},
- "CONUS": {"domain": {"latitude": (24.7, 49.4), "longitude": (-124.78, -66.92)}},
- "land_CONUS": {"value": 100, "domain": {"latitude": (24.7, 49.4), "longitude": (-124.78, -66.92)}},
+ "ocean_50S50N" : {"value":0.,'domain': {"latitude": (-50., 50)}},
+ "ocean_50S20S" : {"value":0.,'domain': {"latitude": (-50., -20)}},
+ "ocean_20S20N": {"value":0.,'domain': {"latitude": (-20., 20)}},
+ "ocean_20N50N" : {"value":0.,'domain': {"latitude": (20., 50)}},
  # Modes of variability
  "NAM": {"domain": {"latitude": (20.0, 90), "longitude": (-180, 180)}},
  "NAO": {"domain": {"latitude": (20.0, 80), "longitude": (-90, 40)}},

pcmdi_metrics/mean_climate/lib/load_and_regrid.py

@@ -49,10 +49,41 @@ def load_and_regrid(data_path, varname, varname_in_file=None, level=None, t_grid
 
  # level - extract a specific level if needed
  if level is not None:
- level = level * 100 # hPa to Pa
- ds = ds.sel(plev=level)
- if debug:
- print('ds:', ds)
+ if isinstance(level, int) or isinstance(level, float):
+ pass
+ else:
+ level = float(level)
+
+ # check vertical coordinate first
+ if 'plev' in list(ds.coords.keys()):
+ if ds.plev.units == 'Pa':
+ level = level * 100 # hPa to Pa
+ try:
+ ds = ds.sel(plev=level)
+ except Exception as ex:
+ print('WARNING: ', ex)
+
+ nearest_level = find_nearest(ds.plev.values, level)
+
+ print(' Given level', level)
+ print(' Selected nearest level from dataset:', nearest_level)
+
+ diff_percentage = abs(nearest_level - level) / level * 100
+ if diff_percentage < 0.1: # acceptable if differance is less than 0.1%
+ ds = ds.sel(plev=level, method='nearest')
+ print(' Difference is in acceptable range.')
+ pass
+ else:
+ print('ERROR: Difference between two levels are too big!')
+ return
+ if debug:
+ print('ds:', ds)
+ print('ds.plev.units:', ds.plev.units)
+ else:
+ print('ERROR: plev is not in the nc file. Check vertical coordinate.')
+ print(' Coordinates keys in the nc file:', list(ds.coords.keys()))
+ print('ERROR: load and regrid can not complete')
+ return
 
  # regrid
  if regrid_tool == 'regrid2':
@@ -78,3 +109,10 @@ def load_and_regrid(data_path, varname, varname_in_file=None, level=None, t_grid
  if debug:
  print('ds_regridded:', ds_regridded)
  return ds_regridded
+
+
+def find_nearest(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return array[idx]
+

pcmdi_metrics/mean_climate/mean_climate_driver.py

@@ -20,6 +20,7 @@
  mean_climate_metrics_to_json,
 )
 from pcmdi_metrics.variability_mode.lib import tree
+from pcmdi_metrics.variability_mode.lib import sort_human
 
 
 parser = create_mean_climate_parser()
@@ -59,15 +60,15 @@
 diagnostics_output_path = diagnostics_output_path.replace('%(case_id)', case_id)
 
 find_all_realizations = False
+first_realization_only = False
 if realization is None:
  realization = ""
- realizations = [realization]
 elif isinstance(realization, str):
  if realization.lower() in ["all", "*"]:
  find_all_realizations = True
-  realizations = "Search for all realizations!!"
- else:
- realizations = [realization]
+ elif realization.lower() in ["first", "first_only"]:
+  first_realization_only = True
+realizations = [realization]
 
 if debug:
  print('regions_specs (before loading internally defined):', regions_specs)
@@ -206,7 +207,7 @@
 
  result_dict["RESULTS"][model][ref]["source"] = ref_dataset_name 
 
- if find_all_realizations:
+ if find_all_realizations or first_realization_only:
  test_data_full_path = os.path.join(
  test_data_path,
  filename_template).replace('%(variable)', varname).replace('%(model)', model).replace('%(model_version)', model).replace('%(realization)', '*')
@@ -215,6 +216,9 @@
  realizations = []
  for ncfile in ncfiles:
  realizations.append(ncfile.split('/')[-1].split('.')[3])
+ realizations = sort_human(realizations)
+ if first_realization_only:
+ realizations = realizations[0:1]
  print('realizations (after search): ', realizations)
 
  for run in realizations: