Per #2486, update grid_proj_comparison utility to read any named grid…

…. Still do need to implement ProjGrid support for non-lambert conformal grids.

internal/test_util/libcode/vx_grid/grid_proj_comparison.cc

@@ -1,7 +1,25 @@
-
+// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
+// ** Copyright UCAR (c) 1992 - 2023
+// ** University Corporation for Atmospheric Research (UCAR)
+// ** National Center for Atmospheric Research (NCAR)
+// ** Research Applications Lab (RAL)
+// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
+// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 
 ////////////////////////////////////////////////////////////////////////
-
+//
+//   Filename:   grid_proj_comparison.cc
+//
+//   Description:
+//      Locate the requested grid by name, instantiate it
+//      using the legacy grid code and the newer Proj-based
+//      grid code, and track the maximum location difference
+//      for each point in the grid.
+//
+//   Mod#   Date      Name            Description
+//   ----   ----      ----            -----------
+//
+////////////////////////////////////////////////////////////////////////
 
 using namespace std;
 
@@ -15,13 +33,23 @@ using namespace std;
 #include "vx_math.h"
 #include "vx_grid.h"
 
+////////////////////////////////////////////////////////////////////////
+
+static ConcatString program_name;
+
+static ConcatString GridName;
 
 ////////////////////////////////////////////////////////////////////////
 
+static void process_command_line(int, char **);
+static void usage();
+static void set_grid(const StringArray &);
+
+////////////////////////////////////////////////////////////////////////
 
-static const LambertData data = {
+static const LambertData default_data = {
 
-   "adam_clark",        // name
+   "LambertGrid",       // name
 
    'N',                 // hemisphere
 
@@ -30,7 +58,6 @@ static const LambertData data = {
     19.834, 120.051,    //  lat_pin, lot_pin
 
     0.0, 0.0,           //  x_pin, y_pin
-    // 5.0, 10.0,           //  x_pin, y_pin
 
     97.0,               //  lon_orient
 
@@ -42,136 +69,154 @@ static const LambertData data = {
 
 };
 
-
 ////////////////////////////////////////////////////////////////////////
 
+int met_main(int argc, char * argv[]) {
 
-int main()
-
-{
-
-const Grid g(data);
-ProjGrid p;
-Affine aff;
-ConcatString proj_params;
-double xx, yy;
-char junk[256];
-
-
-cout << "\n\n";
+   program_name = get_short_name(argv[0]);
 
    //
-   //  create the PROJ part of the transformation
+   // process the command line arguments
    //
+   process_command_line(argc, argv);
+
+   Grid g;
+   ProjGrid p;
 
    //
-   //  Note: PROJ lambert params explained at:
-   //
-   //     https://proj.org/en/9.3/operations/projections/lcc.html#lambert-conformal-conic
+   // instantiate the grids
    //
+   if(GridName.nonempty()) {
 
-   //   specify lambert conformal projection
+      mlog << Debug(1) << "Reading grid named \"" << GridName << "\".\n";
 
-proj_params = "+proj=lcc";   //  lcc means lambert conformal conic
+      if(!find_grid_by_name(GridName.c_str(), g)) {
+         mlog << Error << "\n" << program_name << " -> "
+              << "Legacy Grid name \"" << GridName << "\" not found!\n\n";
+         exit(1);
+      }
 
-   //   radius of earth
-
-sprintf(junk, "+R=%.3f", data.r_km);   //  earth radius
-
-proj_params << ' ' << junk;
-
-   //  the two secant latitudes for the lambert projection
-
-sprintf(junk, "+lat_1=%.3f +lat_2=%.3f", data.scale_lat_1, data.scale_lat_2);
-
-proj_params << ' ' << junk;
-
-   //  orientation longitude
-   //    note: PROJ calls it "Longitude of projection center"
+      if(!find_grid_by_name(GridName.c_str(), p))  {
+         mlog << Error << "\n" << program_name << " -> "
+              << "Proj Grid name \"" << GridName << "\" not found!\n\n";
+         exit(1);
+      }
+   }
+   else {
+      mlog << Debug(1) << "Reading default grid.\n";
+      g.set(default_data);
+      p.set(default_data);
+   }
 
-sprintf(junk, "+lon_0=%.3f", -(data.lon_orient));   //  note minus sign
+   if(mlog.verbosity_level() >= 2) {
+      mlog << Debug(2) << "\nLegacy Grid Definition:\n\n";
+      g.dump(cout, 1);
+      mlog << Debug(2) << "\nProj Grid Definition:\n\n";
+      p.dump(cout, 1);
+   }
 
-proj_params << ' ' << junk;
+   //
+   // find max difference over grid
+   //
+   double px, py, dx, dy, lat, lon, dist;
+   double max_dist = 0.0;
 
-   //  false easting and northing
-   //     (on second thought, we'll let the Affine class handle this)
+   for(int gx=0; gx<(g.nx()); ++gx) {
+      for(int gy=0; gy<(g.ny()); ++gy) {
 
-   //  ok
+         g.xy_to_latlon(gx, gy, lat, lon);
+         p.latlon_to_xy(lat, lon, px, py);
 
-cout << "proj_params = \"" << proj_params << "\"\n\n" << flush;
+         dx = (double) gx - px;
+         dy = (double) gy - py;
 
-p.set_proj(proj_params.c_str());
+         dist = sqrt(dx*dx + dy*dy);
 
-   //
-   //  create the affine part of the transformation
-   //
+         if(mlog.verbosity_level() >= 3) {
+            mlog << Debug(3)
+                 << "(grid_x, grid_y) = (" << gx << ", " << gy << "), "
+                 << "(lat, lon) = (" << lat << ", " << lon << "), "
+                 << "(proj_x, proj_y) = (" << px << ", " << py << "), "
+                 << "dist = " << dist << "\n";
+         }
 
-const double s = 1.0/(data.d_km);   //  scale factor
+         if(dist > max_dist) max_dist = dist;
+      }
+   }
 
-aff.set_mb(s, 0.0, 0.0, s, 0.0, 0.0);
+   if ( !is_eq(max_dist, 0.0) ) {
+      mlog << Warning << "\nMaximum distance = " << max_dist << " > 0\n\n";
+   }
+   else {
+      mlog << Debug(1) << "Maximum distance = " << max_dist << "\n";
+   }
 
-p.set_affine(aff);
+   return(0);
+}
 
+////////////////////////////////////////////////////////////////////////
 
-p.latlon_to_xy(data.lat_pin, data.lon_pin, xx, yy);
+const string get_tool_name() {
+   return "grid_proj_comparison";
+}
 
-aff.set_b(data.x_pin - xx, data.y_pin - yy);
+////////////////////////////////////////////////////////////////////////
 
-p.set_affine(aff);   //  again
+void process_command_line(int argc, char **argv) {
+   CommandLine cline;
 
    //
-   //  misc
+   // parse the command line into tokens
    //
+   cline.set(argc, argv);
 
-p.set_size(data.nx, data.ny);
-
-p.set_name(data.name);
+   //
+   // set the usage function
+   //
+   cline.set_usage(usage);
 
    //
-   //  find max error over grid
+   // add the options function calls
    //
+   cline.add(set_grid, "-grid", 1);
 
-int i, j;
-double er, max_er;
-double lat, lon;
-double px, py;
+   //
+   // parse the command line
+   //
+   cline.parse();
 
+   return;
+}
 
-max_er = 0.0;
+////////////////////////////////////////////////////////////////////////
 
-for (i=0; i<(g.nx()); ++i)  {
+void usage() {
 
-   for (j=0; j<(g.ny()); ++j)  {
+   cout << "\n*** Model Evaluation Tools (MET" << met_version
+        << ") ***\n\n"
 
-      g.xy_to_latlon(i, j, lat, lon);
+        << "Usage: " << program_name << "\n"
 
-      p.latlon_to_xy(lat, lon, px, py);
+        << "\t[-grid name]\n"
+        << "\t[-log file]\n"
+        << "\t[-v level]\n\n"
 
-      er = fabs(i - px) + fabs(j - py);
+        << "\twhere\t\"-grid name\" overrides the default name of "
+        << "grid to be analyzed (optional).\n"
 
-      if ( er > max_er )  max_er = er;
+        << "\t\t\"-log file\" outputs log messages to the specified "
+        << "file (optional).\n"
 
-   }
+        << "\t\t\"-v level\" overrides the default level of logging ("
+        << mlog.verbosity_level() << ") (optional).\n\n" << flush;
 
+   exit(1);
 }
 
+////////////////////////////////////////////////////////////////////////
 
-cout << "max error = " << max_er << "\n\n";
-
-
-
-   //
-   //  done
-   //
-
-return ( 0 );
-
+void set_grid(const StringArray &a) {
+   GridName = a[0];
 }
 
-
 ////////////////////////////////////////////////////////////////////////
-
-
-
-
-