Fix grid North alignment check

src/cpp/fs/Grid.cpp

@@ -7,6 +7,7 @@
 #include "Log.h"
 #include "projection.h"
 #include "tiff.h"
+#include "unstable.h"
 using fs::Idx;
 namespace fs
 {
@@ -391,41 +392,26 @@ std::optional<Coordinates> GridBase::findCoordinates(const Point& point, const b
 std::optional<FullCoordinates> GridBase::findFullCoordinates(const Point& point, const bool flipped)
   const
 {
-  // check that north is the top of the raster at least along center
-  const auto x_mid = XPos{(xllcorner_ + xurcorner_) / 2.0};
-  const auto y_bottom = YPos{yllcorner_};
-  const auto y_top = YPos{yurcorner_};
-  XYPos mid_bottom{x_mid, y_bottom};
-  XYPos mid_top{x_mid, y_top};
   const auto xy = from_lat_long(this->proj4_, point);
   logging::debug("Coordinates {} converted to ({:f}, {:f})", point, xy.x.value, xy.y.value);
-  auto south = to_lat_long(proj4_, mid_bottom);
-  auto north = to_lat_long(proj4_, mid_top);
-  auto grid_south = from_lat_long(this->proj4_, south);
-  auto grid_north = from_lat_long(this->proj4_, north);
   // FIX: how different is too much?
-  constexpr MathSize MAX_DEVIATION = 0.001;
-  const auto deviation = (grid_north.x - grid_south.x).value;
-  if (abs(deviation) > MAX_DEVIATION)
+  constexpr MathSize MAX_DEVIATION{1};
+  if (!check_deviation("origin", proj4_, point, MAX_DEVIATION))
   {
-    logging::note(
-      "Due north is not the top of the raster for {} with proj4 '{:s}' - {:f} vs {:f} gives deviation of {:f} degrees which exceeds maximum of {:f} degrees",
-      point,
-      this->proj4_,
-      grid_north.x.value,
-      grid_south.x.value,
-      deviation,
-      MAX_DEVIATION
-    );
     return {};
   }
-  else if (abs(deviation * 10) > MAX_DEVIATION)
+  const auto lower_left = to_lat_long(this->proj4_, XYPos{XPos{xllcorner_}, YPos{yllcorner_}});
+  // how much it can be off by at edges
+  constexpr auto MAX_DEVIATION_EDGES{10 * MAX_DEVIATION};
+  if (!check_deviation("lower left", proj4_, lower_left, MAX_DEVIATION_EDGES))
   {
-    // if we're within an order of magnitude of an unacceptable deviation then warn about it
-    logging::warning(
-      "Due north deviates by {:f} degrees from South to North along the middle of the raster",
-      deviation
-    );
+    return {};
+  }
+  // goint north of upper corners goes way off the grid, so do lower_right
+  const auto lower_right = to_lat_long(this->proj4_, XYPos{XPos{xurcorner_}, YPos{yllcorner_}});
+  if (!check_deviation("lower right", proj4_, lower_right, MAX_DEVIATION_EDGES))
+  {
+    return {};
   }
   logging::verbose("Lower left is ({:f}, {:f})", this->xllcorner_, this->yllcorner_);
   // convert coordinates into cell position

src/cpp/fs/projection.cpp

@@ -3,6 +3,7 @@
 #include <proj.h>
 #include "Log.h"
 #include "Point.h"
+#include "Radians.h"
 #include "Settings.h"
 #include "unstable.h"
 #include "Util.h"
@@ -191,4 +192,73 @@ string try_fix_meridian(const string_view proj4)
     units
   );
 }
+Degrees find_north_south_deviation(const string_view proj4, const Point& p0)
+{
+  constexpr auto MAX_LATITUDE{89.9};
+  constexpr auto LATITUDE_DISTANCE{1};
+  const auto lat1 = min(MAX_LATITUDE, p0.latitude() + LATITUDE_DISTANCE);
+  Point p1{lat1, p0.longitude()};
+  auto grid0 = from_lat_long(proj4, p0);
+  auto grid1 = from_lat_long(proj4, p1);
+  logging::verbose(
+    "Finding deviation between [{} => ({:f}, {:f})] and [{} => ({:f}, {:f})]",
+    p0,
+    grid0.x.value,
+    grid0.y.value,
+    p1,
+    grid1.x.value,
+    grid1.y.value
+  );
+  // angle is going to be how far off North we are
+  // -90 to convert from math direction
+  auto deviation =
+    Radians::D_090().asDegrees()
+    - Radians{atan2(grid1.y.value - grid0.y.value, grid1.x.value - grid0.x.value)}.asDegrees();
+  logging::verbose(
+    "Deviation for {} to {} with proj4 '{:s}' is {:f} degrees", p0, p1, proj4, deviation.value
+  );
+  return deviation;
+};
+bool check_deviation(
+  const string_view what,
+  const string_view proj4,
+  const Point& p,
+  const MathSize max_deviation
+)
+{
+  auto deviation = find_north_south_deviation(proj4, p).value;
+  if (abs(deviation) > max_deviation)
+  {
+    logging::error(
+      "Grid North of {:s} {} deviates from true North by {:g} degrees which exceeds maximum of {:g} degrees",
+      what,
+      p,
+      deviation,
+      max_deviation
+    );
+    return false;
+  }
+  else if (abs(deviation * 2) > max_deviation)
+  {
+    // if close to max deviation then warn about it
+    logging::warning(
+      "Grid North of {:s} {} deviates from true North by {:g} degrees which is near maximum of {:f} degrees",
+      what,
+      p,
+      deviation,
+      max_deviation
+    );
+  }
+  else
+  {
+    logging::info(
+      "Grid North of {:s} {} deviates from true North by {:g} degrees",
+      what,
+      p,
+      deviation,
+      max_deviation
+    );
+  }
+  return true;
+};
 }

src/cpp/fs/projection.h

@@ -36,5 +36,12 @@ class Point;
 fs::XYPos from_lat_long(const string_view proj4, const fs::Point& point);
 fs::Point to_lat_long(const string_view proj4, const fs::XYPos xy);
 string try_fix_meridian(const string_view proj4);
+Degrees find_north_south_deviation(const string_view proj4, const Point& p0);
+bool check_deviation(
+  const string_view what,
+  const string_view proj4,
+  const Point& p,
+  const MathSize max_deviation
+);
 }
 #endif