Update main_v12.2-ref after #3410

internal/test_unit/python/comp_dir.py

@@ -14,7 +14,7 @@
     raise e
 sys.path.append(metplus_util_path)
 
-from diff_util import compare_dir
+import diff_util as du
 
 
 def comp_dir(truth_dir, output_dir, debug=True, save_diff=True):
@@ -43,11 +43,15 @@ def comp_dir(truth_dir, output_dir, debug=True, save_diff=True):
     """
 
     start_time = dt.datetime.now()
+
+    # Custom rounding overrides for MET unit test output
+    du.ROUNDING_OVERRIDES['ALAL2010_stat.out'] = 4
+
     print('******************************')
     print("Comparing output to truth data")
-    diff_files = compare_dir(truth_dir, output_dir,
-                             debug=debug,
-                             save_diff=save_diff)
+    diff_files = du.compare_dir(truth_dir, output_dir,
+                                debug=debug,
+                                save_diff=save_diff)
 
     end_time = dt.datetime.now()
     runtime = end_time - start_time

src/libcode/vx_data2d_nc_cf/nc_cf_file.cc

@@ -1912,6 +1912,20 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin
     return;
   }
 
+  // false_easting, optional
+
+  double false_east = get_nc_var_att_double(
+    grid_mapping_var, "false_easting", false);
+
+  if(is_bad_data(false_east)) false_east = 0.0;
+
+  // false_northing, optional
+
+  double false_north = get_nc_var_att_double(
+    grid_mapping_var, "false_northing", false);
+
+  if(is_bad_data(false_north)) false_north = 0.0;
+
   // Look for the x/y dimensions and x/y coordinate variables
   find_xy_vars(method_name);
 
@@ -1987,10 +2001,11 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin
 
   get_nc_data(_yCoordVar, y_values.data());
 
-  // Unit conversion
+  // MET #3374 Support false easting and northing offsets
+  // Unit conversion and false northing and easting offsets
 
-  for (int i = 0; i<x_counts; ++i) x_values[i] *= x_coord_to_m_cf;
-  for (int i = 0; i<y_counts; ++i) y_values[i] *= y_coord_to_m_cf;
+  for (int i = 0; i<x_counts; ++i) x_values[i] = x_values[i] * x_coord_to_m_cf - false_east;
+  for (int i = 0; i<y_counts; ++i) y_values[i] = y_values[i] * y_coord_to_m_cf - false_north;
 
   // Calculate dx and dy assuming they are constant.  MET requires that dx be
   // equal to dy

src/libcode/vx_data2d_nc_cf/var_info_nc_cf.cc

@@ -208,8 +208,23 @@ void VarInfoNcCF::set_magic(const ConcatString &nstr, const ConcatString &lstr)
             as_offset = (*ptr2 != '@');
             if (!as_offset) ptr2++;
 
+            ptr3 = strchr(ptr2, '-');
+            //skip negative sign of the negative value to check the range
+            if (ptr3 != nullptr && ptr3 == ptr2) ptr3 = strchr((ptr2+1), '-');
+
             // Check for a range of levels
-            if ((ptr3 = strchr(ptr2, '-')) != nullptr) {
+            if (ptr3 != nullptr && ptr3 != ptr2) {
+               *ptr3 = 0;
+               ptr3++;
+               if (*ptr3 == '@') {
+                  if (as_offset) {
+                     mlog << Error << "\n" << method_name
+                          << "Can not mix an offset and a value for NetCDF variable \""
+                          << MagicStr << "\".\n\n";
+                     exit(1);
+                  }
+                  ptr3++;    // to support @vlevel_lower-@vlevel_upper
+               }
 
                // Check if a range has already been supplied
                if (Dimension.has(range_flag)) {

src/libcode/vx_data2d_nc_wrf/data2d_nc_wrf.cc

@@ -121,33 +121,75 @@ void MetNcWrfDataFile::dump(ostream & out, int depth) const {
 
 ////////////////////////////////////////////////////////////////////////
 
+bool MetNcWrfDataFile::get_real_dimension(const VarInfoNcWrf *vinfo_nc,
+                                          const NcVarInfo *info, LongArray &dimension) const {
+   static const string method_name
+         = "MetNcWrfDataFile::get_real_dimension() ->";
+   int dim_count = dimension.n_elements();
+   for (int k=0; k<dim_count; k++) {
+      if (dimension[k] == range_flag && vinfo_nc->level().is_offset()) {
+         dimension[k] = nint(vinfo_nc->level().lower());
+         continue;
+      }
+
+      if (dimension[k] != vx_data2d_dim_by_value && dimension[k] != range_flag) continue;
+
+      string dim_name = GET_NC_NAME(get_nc_dim(info->var, k));
+      NcVarInfo *var_info = find_var_info_by_dim_name(WrfNc->Var, dim_name,
+                                                         WrfNc->Nvars);
+      if (var_info == nullptr) continue;
+
+      if (dimension[k] == vx_data2d_dim_by_value) {
+         long new_offset = get_index_at_nc_data(var_info->var,
+                                                vinfo_nc->dim_value(k),
+                                                dim_name, (k == info->t_slot));
+         if (new_offset != bad_data_int) dimension[k] = new_offset;
+         else {
+            mlog << Warning << "\n" << method_name
+                 << "for \"" << vinfo_nc->req_name()
+                 << "\" variable, the dimension value "
+                 << vinfo_nc->dim_value(k) << " for " << dim_name
+                 << " does not exist\n\n";
+            return false;
+         }
+      }
+      else if (dimension[k] == range_flag) {
+         double lower = vinfo_nc->level().lower();
+         double upper = vinfo_nc->level().upper();
+         long new_offset = get_index_at_nc_data(var_info->var, lower, upper,
+                                                dim_name, (k == info->t_slot));
+         if (new_offset != bad_data_int) dimension[k] = new_offset;
+         else {
+            mlog << Warning << "\n" << method_name
+                 << "for \"" << vinfo_nc->req_name()
+                 << "\" variable, the dimension value between "
+                 << lower <<  " and " << upper << " for " << dim_name
+                 << " does not exist\n\n";
+            return false;
+         }
+      }
+   }
+   return true;
+ }
+
+////////////////////////////////////////////////////////////////////////
+
 bool MetNcWrfDataFile::data_plane(VarInfo &vinfo, DataPlane &plane) {
    bool status = false;
    double pressure;
    ConcatString level_str;
-   VarInfoNcWrf * vinfo_nc = (VarInfoNcWrf *) &vinfo;
-   NcVarInfo *info = (NcVarInfo *) nullptr;
+   auto vinfo_nc = (VarInfoNcWrf *) &vinfo;
+   auto info = (NcVarInfo *) nullptr;
 
    // Initialize the data plane
    plane.clear();
 
    // Read the data
    WrfNc->get_nc_var_info(vinfo_nc->req_name().c_str(), info);
+
    LongArray dimension = vinfo_nc->dimension();
-   int dim_count = dimension.n_elements();
-   for (int k=0; k<dim_count; k++) {
-      if (dimension[k] == vx_data2d_dim_by_value) {
-         string dim_name = GET_NC_NAME(get_nc_dim(info->var, k));
-         NcVarInfo *var_info = find_var_info_by_dim_name(WrfNc->Var, dim_name,
-                                                         WrfNc->Nvars);
-         if (var_info) {
-            long new_offset = get_index_at_nc_data(var_info->var,
-                                                   vinfo_nc->dim_value(k),
-                                                   dim_name, (k == info->t_slot));
-            if (new_offset != bad_data_int) dimension[k] = new_offset;
-         }
-      }
-   }
+   if (! get_real_dimension(vinfo_nc, info, dimension))
+      return false;
 
    status = WrfNc->data(vinfo_nc->req_name().c_str(),
                         dimension, plane, pressure, info);
@@ -198,7 +240,7 @@ bool MetNcWrfDataFile::data_plane(VarInfo &vinfo, DataPlane &plane) {
 ////////////////////////////////////////////////////////////////////////
 
 int MetNcWrfDataFile::data_plane_array(VarInfo &vinfo,
-                                           DataPlaneArray &plane_array) {
+                                       DataPlaneArray &plane_array) {
    int i, i_dim, n_level, status, lower, upper;
    ConcatString level_str;
    double pressure, min_level, max_level;

src/libcode/vx_data2d_nc_wrf/data2d_nc_wrf.h

@@ -34,6 +34,8 @@ class MetNcWrfDataFile : public Met2dDataFile {
    private:
 
       void nc_wrf_init_from_scratch();
+      bool get_real_dimension(const VarInfoNcWrf *vinfo_nc, const NcVarInfo *info,
+                              LongArray &dimension) const;
 
       MetNcWrfDataFile(const MetNcWrfDataFile &);
       MetNcWrfDataFile & operator=(const MetNcWrfDataFile &);

src/libcode/vx_data2d_nc_wrf/var_info_nc_wrf.cc

@@ -209,8 +209,23 @@ void VarInfoNcWrf::set_magic(const ConcatString &nstr, const ConcatString &lstr)
             as_offset = (*ptr2 != '@');
             if (!as_offset) ptr2++;
 
+            ptr3 = strchr(ptr2, '-');
+            //skip negative sign of the negative value to check the range
+            if (ptr3 != nullptr && ptr3 == ptr2) ptr3 = strchr((ptr2+1), '-');
+
             // Check for a range of levels
-            if((ptr3 = strchr(ptr2, '-')) != nullptr) {
+            if (ptr3 != nullptr && ptr3 != ptr2) {
+               *ptr3 = 0;
+               ptr3++;
+               if (*ptr3 == '@') {
+                  if (as_offset) {
+                     mlog << Error << "\n" << method_name
+                          << "Can not mix an offset and a value for NetCDF variable \""
+                          << MagicStr << "\".\n\n";
+                     exit(1);
+                  }
+                  ptr3++;    // to support @vlevel_lower-@vlevel_upper
+               }
 
                // Check if a range has already been supplied
                if(Dimension.has(range_flag)) {
@@ -223,10 +238,11 @@ void VarInfoNcWrf::set_magic(const ConcatString &nstr, const ConcatString &lstr)
                else {
                   add_dimension(range_flag);
                   Level.set_lower(atoi(ptr2));
-                  Level.set_upper(atoi(++ptr3));
+                  Level.set_upper(atoi(ptr3));
 
                   // Assume pressure level type for a range of levels
                   Level.set_type(LevelType_Pres);
+                  Level.set_is_offset(as_offset);
                }
             }
             // Single level

src/libcode/vx_data2d_ugrid/ugrid_file.cc

@@ -38,6 +38,7 @@ using namespace netCDF;
 static const char *def_user_config = "UGridConfig_user";
 static const char *def_config_prefix = "UGridConfig_";
 static const char *def_config_prefix2 = "MET_BASE/config/UGridConfig_";
+constexpr double lat_epsilon = 0.00001;
 
 array<string, UG_DIM_COUNT> DIM_KEYS = {
       "dim_face", "dim_node", "dim_edge", "dim_time", "dim_vert"
@@ -890,6 +891,40 @@ void UGridFile::read_config(ConcatString config_filename) {
 
 }
 
+////////////////////////////////////////////////////////////////////////
+
+void UGridFile::radian_to_degree(vector<double> &lat_values, const int lat_count) const {
+  const char *method_name = "UGridFile::radian_to_degree() -> ";
+  int lat_adjusted = 0;
+  int lat_adjusted_total = 0;
+  for (int idx=0; idx<lat_count; idx++) {
+    lat_values[idx] /= rad_per_deg;
+    if (lat_values[idx] > 90.0) {
+      if (!is_eq(lat_values[idx], 90.0, lat_epsilon)) {
+        mlog << Warning << "\n" << method_name << "adjusted " << lat_values[idx]
+             << " (delta: " << (lat_values[idx] - 90.0) << ") to 90.0\n\n";
+        lat_adjusted++;
+      }
+      lat_values[idx] = 90.0;
+      lat_adjusted_total++;
+    }
+    else if (lat_values[idx] < -90.0) {
+      if (!is_eq(lat_values[idx], -90.0, lat_epsilon)) {
+        mlog << Warning << "\n" << method_name << "adjusted " << lat_values[idx]
+             << " (delta: " << (lat_values[idx] + 90.0) << ") to -90.0\n\n";
+        lat_adjusted++;
+      }
+      lat_values[idx] = -90.0;
+      lat_adjusted_total++;
+    }
+  }
+  if (lat_adjusted_total > 0) {
+    mlog << Debug(4) << method_name << "adjusted " << lat_adjusted << " ("
+         << lat_adjusted_total << ") latitudes\n";
+  }
+}
+
+
 ////////////////////////////////////////////////////////////////////////
 
 
@@ -927,13 +962,13 @@ void UGridFile::read_netcdf_grid()
 
   if (get_var_units(_latVar, units_value)) {
     if (units_value == "rad" || units_value == "radian") {
-      mlog << Debug(6) << method_name << "convert  " << units_value << " to degree for lat\n";
-      for (int idx=0; idx<face_count; idx++) _lat[idx] /= rad_per_deg;
+      mlog << Debug(6) << method_name << "convert " << units_value << " to degree for lat\n";
+      radian_to_degree(_lat, face_count);
     }
   }
   if (get_var_units(_lonVar, units_value)) {
     if (units_value == "rad" || units_value == "radian") {
-      mlog << Debug(6) << method_name << "  convert " << units_value << " to degree for lon\n";
+      mlog << Debug(6) << method_name << "convert " << units_value << " to degree for lon\n";
       for (int idx=0; idx<face_count; idx++) _lon[idx] /= rad_per_deg;
     }
   }

src/libcode/vx_data2d_ugrid/ugrid_file.h

@@ -165,6 +165,7 @@ class UGridFile {
 
       std::string find_metadata_name(std::string &key, StringArray &available_names);
       StringArray get_metadata_names(std::string &key);
+      void radian_to_degree(std::vector<double> &lat_values, const int lat_count) const;
       void read_config(ConcatString config_filename);
       void read_netcdf_grid();