Fix north alignment check

src/cpp/fs/Environment.cpp

@@ -59,7 +59,7 @@ Environment Environment::loadEnvironment(
   const YearSize year
 )
 {
-  logging::note("Using ignition point ({:f}, {:f})", point.latitude(), point.longitude());
+  logging::note("Using ignition point {}", point);
   logging::info("Running using inputs directory '{:s}'", string(path));
   auto rasters = find_rasters(path, year);
   auto best_score = numeric_limits<MathSize>::min();
@@ -142,18 +142,8 @@ Environment Environment::loadEnvironment(
     logging::note("Loading info for fuel {:s}", best_fuel);
     env_info = EnvironmentInfo::loadInfo(best_fuel, best_elevation);
   }
-  logging::check_fatal(
-    nullptr == env_info,
-    "Could not find an environment to use for ({:f}, {:f})",
-    point.latitude(),
-    point.longitude()
-  );
-  logging::debug(
-    "Best match for ({:f}, {:f}) has projection '{:s}'",
-    point.latitude(),
-    point.longitude(),
-    env_info->proj4()
-  );
+  logging::check_fatal(nullptr == env_info, "Could not find an environment to use for {}", point);
+  logging::debug("Best match for {} has projection '{:s}'", point, env_info->proj4());
   logging::note("Projection is {:s}", env_info->proj4());
   // envInfo should get deleted automatically because it uses unique_ptr
   return env_info->load(point);

src/cpp/fs/Grid.cpp

@@ -5,8 +5,9 @@
 #include <tiffio.h>
 #include <xtiffio.h>
 #include "Log.h"
+#include "projection.h"
 #include "tiff.h"
-#include "UTM.h"
+#include "unstable.h"
 using fs::Idx;
 namespace fs
 {
@@ -28,8 +29,10 @@ string saveToTiffFile(
   const int nodata_as_int
 )
 {
-  uint32_t tileWidth = min(static_cast<int>(width), 256);
-  uint32_t tileHeight = min(static_cast<int>(height), 256);
+  // https://gdal.org/en/stable/drivers/raster/gtiff.html
+  // tile width/height must be divisible by 16
+  uint32_t tileWidth = min(static_cast<int>(ceil(width / 16.0) * 16), 256);
+  uint32_t tileHeight = min(static_cast<int>(ceil(height / 16.0) * 16), 256);
   auto [min_x, min_y, max_x, max_y] = bounds;
   // auto min_x = std::get<0>(bounds);
   // auto min_y = std::get<1>(bounds);
@@ -84,10 +87,16 @@ string saveToTiffFile(
   const auto width_calc = static_cast<size_t>(max_x - min_x);
   // ensure this is always divisible by tile size
   logging::check_fatal(
-    0 != (height_calc % tileWidth), "Height {:d} not divisible by tiles", height_calc
+    0 != (height_calc % tileHeight),
+    "Height {:d} not divisible by tile height {:d}",
+    height_calc,
+    tileHeight
   );
   logging::check_fatal(
-    0 != (width_calc % tileHeight), "Width {:d} not divisible by tiles", width_calc
+    0 != (width_calc % tileWidth),
+    "Width {:d} not divisible by tile width {:d}",
+    width_calc,
+    tileWidth
   );
   const auto filename = create_file_name(dir, base_name, "tif");
   GeoTiff geotiff{filename, "w"};
@@ -383,53 +392,33 @@ std::optional<Coordinates> GridBase::findCoordinates(const Point& point, const b
 std::optional<FullCoordinates> GridBase::findFullCoordinates(const Point& point, const bool flipped)
   const
 {
-  MathSize x;
-  MathSize y;
-  from_lat_long(this->proj4_, point, &x, &y);
-  logging::debug(
-    "Coordinates ({:f}, {:f}) converted to ({:f}, {:f})", point.latitude(), point.longitude(), x, y
-  );
-  // check that north is the top of the raster at least along center
-  const auto x_mid = (xllcorner_ + xurcorner_) / 2.0;
-  Point south = to_lat_long(proj4_, x_mid, yllcorner_);
-  Point north = to_lat_long(proj4_, x_mid, yurcorner_);
-  auto x_s = static_cast<MathSize>(0.0);
-  auto y_s = static_cast<MathSize>(0.0);
-  from_lat_long(this->proj4_, south, &x_s, &y_s);
-  auto x_n = static_cast<MathSize>(0.0);
-  auto y_n = static_cast<MathSize>(0.0);
-  from_lat_long(this->proj4_, north, &x_n, &y_n);
+  const auto xy = from_lat_long(this->proj4_, point);
+  logging::debug("Coordinates {} converted to ({:f}, {:f})", point, xy.x.value, xy.y.value);
   // FIX: how different is too much?
-  constexpr MathSize MAX_DEVIATION = 0.001;
-  const auto deviation = x_n - x_s;
-  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 ({:f}, {:f}) with proj4 '{:s}' - {:f} vs {:f} gives deviation of {:f} degrees which exceeds maximum of {:f} degrees",
-      point.latitude(),
-      point.longitude(),
-      this->proj4_,
-      x_n,
-      x_s,
-      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
-  const auto actual_x = (x - this->xllcorner_) / this->cell_size_;
+  const auto actual_x = (xy.x.value - this->xllcorner_) / this->cell_size_;
   // these are already flipped across the y-axis on reading, so it's the same as for x now
-  auto actual_y =
-    (!flipped) ? (y - this->yllcorner_) / this->cell_size_ : (yurcorner_ - y) / cell_size_;
+  auto actual_y = (!flipped) ? (xy.y.value - this->yllcorner_) / this->cell_size_
+                             : (yurcorner_ - xy.y.value) / cell_size_;
   const auto x1 = static_cast<FullIdx>(actual_x);
   const auto y1 = static_cast<FullIdx>(round(actual_y - 0.5));
   if (0 > x1 || x1 >= calculateWidth() || 0 > y1 || y1 >= calculateHeight())

src/cpp/fs/Log.cpp

@@ -1,5 +1,6 @@
 /* SPDX-License-Identifier: AGPL-3.0-or-later */
 #include "Log.h"
+#include <tiffio.h>
 #ifdef NDEBUG
 #ifdef _WIN32
 #include "TimeUtil.h"
@@ -54,6 +55,21 @@ ofstream* out_file() noexcept
 }
 mutex mutex_{};
 ostream* output_{pre_file_log()};
+// HACK: rely on error handler being set when log is opened
+void handle_tiff_error(const char* module, const char* fmt, va_list args)
+{
+  const string tmp;
+  stringstream iss(tmp);
+  // HACK: just use a static size instead of trying to determine
+  char buffer[1024]{0};
+  vsnprintf(buffer, std::size(buffer), fmt, args);
+  if (module && 0 < strlen(module))
+  {
+    iss << module << ": ";
+  }
+  iss << buffer;
+  std::ignore = output(logging::level::error, "{}", iss.str());
+}
 int open_log_file(const char* filename) noexcept
 {
   lock_guard<mutex> lock(mutex_);
@@ -73,6 +89,8 @@ int open_log_file(const char* filename) noexcept
     prefile->clear();
   }
   output_ = outfile;
+  // HACK: redirect TIFF errors here since log is always opened
+  TIFFSetErrorHandler(&handle_tiff_error);
   return true;
 }
 int close_log_file() noexcept

src/cpp/fs/Model.cpp

@@ -37,7 +37,7 @@ Model::Model(
     active_simulations_still_required_(settings.minimum_active_simulation_count),
     latitude_(start_point.latitude()), longitude_(start_point.longitude())
 {
-  logging::debug("Calculating for ({:f}, {:f})", start_point.latitude(), start_point.longitude());
+  logging::debug("Calculating for {}", start_point);
   const auto nd_for_point = calculate_nd_ref_for_point(env->elevation(), start_point);
   for (auto day = 0; day < MAX_DAYS; ++day)
   {
@@ -853,13 +853,7 @@ map<DurationSize, shared_ptr<ProbabilityMap>> Model::runIterations(
   static_assert(std::numeric_limits<size_t>::digits10 >= precision);
   const auto lat = static_cast<size_t>(abs(start_point.latitude()) * pow(10, precision));
   const auto lon = static_cast<size_t>(abs(start_point.longitude()) * pow(10, precision));
-  logging::debug(
-    "lat/long ({:f}, {:f}) converted to ({:d}, {:d})",
-    start_point.latitude(),
-    start_point.longitude(),
-    lat,
-    lon
-  );
+  logging::debug("lat/long {} converted to ({:d}, {:d})", start_point, lat, lon);
   const size_t base_salt = settings.salt;
   auto make_seed = [&](const char* name, const size_t salt) {
     const auto d = static_cast<size_t>(start_day);

src/cpp/fs/Point.h

@@ -45,4 +45,25 @@ class Point
   MathSize longitude_;
 };
 }
+template <typename T>
+concept LatLong = requires(const T& p) {
+  { p.latitude() } noexcept -> std::same_as<fs::MathSize>;
+  { p.longitude() } noexcept -> std::same_as<fs::MathSize>;
+};
+template <class T>
+  requires LatLong<T>
+struct std::formatter<T> : std::formatter<string_view>
+{
+  constexpr auto parse(std::format_parse_context& ctx) { return ctx.begin(); }
+  // apple clang didn't work with std::format_context
+  template <class FormatContext>
+  auto format(const T& p, FormatContext& ctx) const
+  {
+    std::string tmp{};
+    std::format_to(
+      std::back_inserter(tmp), "({:f}\u00b0, {:f}\u00b0)", p.latitude(), p.longitude()
+    );
+    return std::formatter<string_view>::format(tmp, ctx);
+  }
+};
 #endif

src/cpp/fs/SimpleFBP.cpp

@@ -532,10 +532,9 @@ vector<int> find_nd_values()
         const auto nd_ref = calculate_nd_ref_for_point(elevation, pt);
         nd_ref_values.emplace(nd_ref);
         logging::verbose(
-          "now have {:d} values for nd: {:g}, {:g}, {:g} gives nd_ref {:d}",
+          "now have {:d} values for nd: {}, {:g} gives nd_ref {:d}",
           nd_ref_values.size(),
-          latitude,
-          longitude,
+          pt,
           elevation,
           nd_ref
         );

src/cpp/fs/Statistics.h

@@ -105,7 +105,8 @@ class Statistics
    * \brief Calculates statistics on a vector of values
    * \param values Values to use for calculation
    */
-  explicit Statistics(vector<MathSize> values)
+  template <class T>
+  explicit Statistics(vector<T> values)
     // values should already be sorted
     : n_(values.size()), min_(values.empty() ? std::numeric_limits<MathSize>::max() : values[0]),
       max_(values.empty() ? std::numeric_limits<MathSize>::min() : values[n_ - 1]), mean_([&]() {