From f8a0f4773141cc978e09d75f1ee01ac322083522 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Wed, 17 Jun 2026 10:16:18 -0400 Subject: [PATCH 01/11] FIX: tile width/height not divisible by 16 when input raster is small (#29) Fix #30 * fix check for being divisible by tile width/height * FIX: tile width/height not divisible by 16 when input raster is small --- src/cpp/fs/Grid.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index a44f3d0ebd..ea9aa2cddc 100644 --- a/src/cpp/fs/Grid.cpp +++ b/src/cpp/fs/Grid.cpp @@ -28,8 +28,10 @@ string saveToTiffFile( const int nodata_as_int ) { - uint32_t tileWidth = min(static_cast(width), 256); - uint32_t tileHeight = min(static_cast(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(ceil(width / 16.0) * 16), 256); + uint32_t tileHeight = min(static_cast(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 +86,16 @@ string saveToTiffFile( const auto width_calc = static_cast(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"}; From 1283e44e62f1c310762ca826039d50fd9131c8a8 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Wed, 17 Jun 2026 10:58:35 -0400 Subject: [PATCH 02/11] FIX: Some errors are not redirected to the log file #31 --- src/cpp/fs/Log.cpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/cpp/fs/Log.cpp b/src/cpp/fs/Log.cpp index 077c033d10..47056a4a8c 100644 --- a/src/cpp/fs/Log.cpp +++ b/src/cpp/fs/Log.cpp @@ -1,5 +1,6 @@ /* SPDX-License-Identifier: AGPL-3.0-or-later */ #include "Log.h" +#include #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 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 From e0806a6aa7e9f1a792773fb9d89714cd6d36141b Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Wed, 10 Jun 2026 19:22:15 +0000 Subject: [PATCH 03/11] add template to Statistics constructor --- src/cpp/fs/Statistics.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cpp/fs/Statistics.h b/src/cpp/fs/Statistics.h index 7713c53588..66b2dcbe99 100644 --- a/src/cpp/fs/Statistics.h +++ b/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 values) + template + explicit Statistics(vector values) // values should already be sorted : n_(values.size()), min_(values.empty() ? std::numeric_limits::max() : values[0]), max_(values.empty() ? std::numeric_limits::min() : values[n_ - 1]), mean_([&]() { From 98dfe745de98df0af86df01aa1fa189b40f69c69 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 09:50:04 -0400 Subject: [PATCH 04/11] Refactor/projection (#33) * combine UTM/project.h/cpp into projection.h/cpp * return Point/XYPos from functions instead of passing x/y pointers --- src/cpp/fs/Grid.cpp | 44 ++++++++-------- src/cpp/fs/project.cpp | 62 ----------------------- src/cpp/fs/project.h | 17 ------- src/cpp/fs/{UTM.cpp => projection.cpp} | 70 +++++++++++++++++++++++--- src/cpp/fs/{UTM.h => projection.h} | 18 +++++-- 5 files changed, 100 insertions(+), 111 deletions(-) delete mode 100644 src/cpp/fs/project.cpp delete mode 100644 src/cpp/fs/project.h rename src/cpp/fs/{UTM.cpp => projection.cpp} (70%) rename src/cpp/fs/{UTM.h => projection.h} (66%) diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index ea9aa2cddc..70176823f2 100644 --- a/src/cpp/fs/Grid.cpp +++ b/src/cpp/fs/Grid.cpp @@ -5,8 +5,8 @@ #include #include #include "Log.h" +#include "projection.h" #include "tiff.h" -#include "UTM.h" using fs::Idx; namespace fs { @@ -391,25 +391,27 @@ std::optional GridBase::findCoordinates(const Point& point, const b std::optional GridBase::findFullCoordinates(const Point& point, const bool flipped) const { - MathSize x; - MathSize y; - from_lat_long(this->proj4_, point, &x, &y); + // 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 ({:f}, {:f}) converted to ({:f}, {:f})", point.latitude(), point.longitude(), x, y + "Coordinates ({:f}, {:f}) converted to ({:f}, {:f})", + point.latitude(), + point.longitude(), + xy.x.value, + xy.y.value ); - // 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(0.0); - auto y_s = static_cast(0.0); - from_lat_long(this->proj4_, south, &x_s, &y_s); - auto x_n = static_cast(0.0); - auto y_n = static_cast(0.0); - from_lat_long(this->proj4_, north, &x_n, &y_n); + 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 = x_n - x_s; + const auto deviation = (grid_north.x - grid_south.x).value; if (abs(deviation) > MAX_DEVIATION) { logging::note( @@ -417,8 +419,8 @@ std::optional GridBase::findFullCoordinates(const Point& point, point.latitude(), point.longitude(), this->proj4_, - x_n, - x_s, + grid_north.x.value, + grid_south.x.value, deviation, MAX_DEVIATION ); @@ -434,10 +436,10 @@ std::optional GridBase::findFullCoordinates(const Point& point, } 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(actual_x); const auto y1 = static_cast(round(actual_y - 0.5)); if (0 > x1 || x1 >= calculateWidth() || 0 > y1 || y1 >= calculateHeight()) diff --git a/src/cpp/fs/project.cpp b/src/cpp/fs/project.cpp deleted file mode 100644 index 719bc7243a..0000000000 --- a/src/cpp/fs/project.cpp +++ /dev/null @@ -1,62 +0,0 @@ -/* SPDX-License-Identifier: AGPL-3.0-or-later */ -#include "project.h" -#include -namespace fs -{ -std::optional to_proj4( - const string& proj4, - const fs::Point& point, - MathSize* x, - MathSize* y -) -{ - PJ_CONTEXT* C; - PJ* P; - PJ* norm; - PJ_COORD a, b, c; - /* or you may set C=PJ_DEFAULT_CTX if you are sure you will */ - /* use PJ objects from only one thread */ - C = proj_context_create(); - P = proj_create_crs_to_crs(C, "EPSG:4326", proj4.c_str(), nullptr); - if (nullptr == P) - { - std::cerr << "Failed to create transformation object.\n"; - return {}; - } - /* This will ensure that the order of coordinates for the input CRS */ - /* will be longitude, latitude, whereas EPSG:4326 mandates latitude, */ - /* longitude */ - norm = proj_normalize_for_visualization(C, P); - if (nullptr == norm) - { - std::cerr << "Failed to normalize transformation object.\n"; - return {}; - } - proj_destroy(P); - P = norm; - /* Given that we have used proj_normalize_for_visualization(), the order */ - /* of coordinates is longitude, latitude, and values are expressed in */ - /* degrees. */ - a = proj_coord(point.longitude(), point.latitude(), 0, 0); - /* transform to UTM, then back to geographical */ - b = proj_trans(P, PJ_FWD, a); - *x = b.enu.e; - *y = b.enu.n; - c = proj_trans(P, PJ_INV, b); - cout << std::format( - "longitude: {:f}, latitude: {:f} => easting: {:.3f}, northing: {:.3f} => x: {:.3f}, y: {:.3f} => longitude: {:g}, latitude: {:g}\n", - point.longitude(), - point.latitude(), - b.enu.e, - b.enu.n, - b.xy.x, - b.xy.y, - c.lp.lam, - c.lp.phi - ); - /* Clean up */ - proj_destroy(P); - proj_context_destroy(C); /* may be omitted in the single threaded case */ - return {}; -} -} diff --git a/src/cpp/fs/project.h b/src/cpp/fs/project.h deleted file mode 100644 index f6699d1f33..0000000000 --- a/src/cpp/fs/project.h +++ /dev/null @@ -1,17 +0,0 @@ -/* SPDX-License-Identifier: AGPL-3.0-or-later */ -#ifndef FS_PROJECT_H -#define FS_PROJECT_H -#include "stdafx.h" -#include "Point.h" -namespace fs -{ -using fs::FullCoordinates; -using fs::MathSize; -std::optional to_proj4( - const string& proj4, - const fs::Point& point, - MathSize* x, - MathSize* y -); -} -#endif diff --git a/src/cpp/fs/UTM.cpp b/src/cpp/fs/projection.cpp similarity index 70% rename from src/cpp/fs/UTM.cpp rename to src/cpp/fs/projection.cpp index 7042af4278..2a0993e7dd 100644 --- a/src/cpp/fs/UTM.cpp +++ b/src/cpp/fs/projection.cpp @@ -1,5 +1,5 @@ /* SPDX-License-Identifier: AGPL-3.0-or-later */ -#include "UTM.h" +#include "projection.h" #include #include "Log.h" #include "Point.h" @@ -8,6 +8,62 @@ #include "Util.h" namespace fs { +std::optional to_proj4( + const string& proj4, + const fs::Point& point, + MathSize* x, + MathSize* y +) +{ + PJ_CONTEXT* C; + PJ* P; + PJ* norm; + PJ_COORD a, b, c; + /* or you may set C=PJ_DEFAULT_CTX if you are sure you will */ + /* use PJ objects from only one thread */ + C = proj_context_create(); + P = proj_create_crs_to_crs(C, "EPSG:4326", proj4.c_str(), nullptr); + if (nullptr == P) + { + std::cerr << "Failed to create transformation object.\n"; + return {}; + } + /* This will ensure that the order of coordinates for the input CRS */ + /* will be longitude, latitude, whereas EPSG:4326 mandates latitude, */ + /* longitude */ + norm = proj_normalize_for_visualization(C, P); + if (nullptr == norm) + { + std::cerr << "Failed to normalize transformation object.\n"; + return {}; + } + proj_destroy(P); + P = norm; + /* Given that we have used proj_normalize_for_visualization(), the order */ + /* of coordinates is longitude, latitude, and values are expressed in */ + /* degrees. */ + a = proj_coord(point.longitude(), point.latitude(), 0, 0); + /* transform to UTM, then back to geographical */ + b = proj_trans(P, PJ_FWD, a); + *x = b.enu.e; + *y = b.enu.n; + c = proj_trans(P, PJ_INV, b); + cout << std::format( + "longitude: {:f}, latitude: {:f} => easting: {:.3f}, northing: {:.3f} => x: {:.3f}, y: {:.3f} => longitude: {:g}, latitude: {:g}\n", + point.longitude(), + point.latitude(), + b.enu.e, + b.enu.n, + b.xy.x, + b.xy.y, + c.lp.lam, + c.lp.phi + ); + /* Clean up */ + proj_destroy(P); + proj_context_destroy(C); /* may be omitted in the single threaded case */ + return {}; +} PJ_CONTEXT* get_context() { // HACK: resolve once and fail if not set already @@ -41,7 +97,7 @@ PJ* normalized_context(PJ_CONTEXT* C, const string_view proj4_from, const string proj_destroy(P); return P_norm; } -void from_lat_long(const string_view proj4, const fs::Point& point, MathSize* x, MathSize* y) +fs::XYPos from_lat_long(const string_view proj4, const fs::Point& point) { // see https://proj.org/en/stable/development/quickstart.html // do this in a function so we can hide and clean up intial context @@ -53,8 +109,7 @@ void from_lat_long(const string_view proj4, const fs::Point& point, MathSize* x, const PJ_COORD a = proj_coord(point.longitude(), point.latitude(), 0, 0); // transform to UTM, then back to geographical const PJ_COORD b = proj_trans(P, PJ_FWD, a); - *x = b.enu.e; - *y = b.enu.n; + XYPos result{XPos{b.enu.e}, YPos{b.enu.n}}; // #ifdef DEBUG_PROJ PJ_COORD c = proj_trans(P, PJ_INV, b); fs::logging::verbose( @@ -71,8 +126,9 @@ void from_lat_long(const string_view proj4, const fs::Point& point, MathSize* x, // #endif proj_destroy(P); proj_context_destroy(C); + return result; } -fs::Point to_lat_long(const string_view proj4, const MathSize x, const MathSize y) +fs::Point to_lat_long(const string_view proj4, const fs::XYPos xy) { // see https://proj.org/en/stable/development/quickstart.html // do this in a function so we can hide and clean up intial context @@ -81,8 +137,8 @@ fs::Point to_lat_long(const string_view proj4, const MathSize x, const MathSize // Given that we have used proj_normalize_for_visualization(), the order // of coordinates is longitude, latitude, and values are expressed in // degrees. - logging::verbose("proj_coord({:f}, {:f}, 0, 0)", x, y); - const PJ_COORD a = proj_coord(x, y, 0, 0); + logging::verbose("proj_coord({:f}, {:f}, 0, 0)", xy.x.value, xy.y.value); + const PJ_COORD a = proj_coord(xy.x.value, xy.y.value, 0, 0); // transform to geographical const PJ_COORD b = proj_trans(P, PJ_FWD, a); // Point is (lat, lon) diff --git a/src/cpp/fs/UTM.h b/src/cpp/fs/projection.h similarity index 66% rename from src/cpp/fs/UTM.h rename to src/cpp/fs/projection.h index 6fb60b4b16..ee4416fc14 100644 --- a/src/cpp/fs/UTM.h +++ b/src/cpp/fs/projection.h @@ -1,9 +1,19 @@ /* SPDX-License-Identifier: AGPL-3.0-or-later */ -#ifndef FS_UTM_H -#define FS_UTM_H +#ifndef FS_PROJECTION_H +#define FS_PROJECTION_H #include "stdafx.h" +#include "Location.h" +#include "Point.h" namespace fs { +using fs::FullCoordinates; +using fs::MathSize; +std::optional to_proj4( + const string& proj4, + const fs::Point& point, + MathSize* x, + MathSize* y +); class Point; /** * \brief Calculate the UTM zone for the given meridian @@ -23,8 +33,8 @@ class Point; { return -183.0 + zone * 6.0; } -void from_lat_long(const string_view proj4, const fs::Point& point, MathSize* x, MathSize* y); -fs::Point to_lat_long(const string_view proj4, const MathSize x, const MathSize y); +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); } #endif From 4a2b58d62e226d67d669822fb043195e7bcdc011 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 13:38:54 -0400 Subject: [PATCH 05/11] add std::formatter for lat/long with degrees symbol (#34) * add std::formatter for Point/StartPoint that uses pair with degrees symbol * use concept for formatter * template to make clang work --- src/cpp/fs/Environment.cpp | 16 +++------------- src/cpp/fs/Grid.cpp | 13 +++---------- src/cpp/fs/Model.cpp | 10 ++-------- src/cpp/fs/Point.h | 21 +++++++++++++++++++++ src/cpp/fs/SimpleFBP.cpp | 5 ++--- src/cpp/fs/projection.cpp | 10 ++++------ 6 files changed, 35 insertions(+), 40 deletions(-) diff --git a/src/cpp/fs/Environment.cpp b/src/cpp/fs/Environment.cpp index 159e787417..e1fc1b0ef3 100644 --- a/src/cpp/fs/Environment.cpp +++ b/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::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); diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index 70176823f2..c5885e0c28 100644 --- a/src/cpp/fs/Grid.cpp +++ b/src/cpp/fs/Grid.cpp @@ -398,13 +398,7 @@ std::optional GridBase::findFullCoordinates(const Point& point, 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 ({:f}, {:f}) converted to ({:f}, {:f})", - point.latitude(), - point.longitude(), - xy.x.value, - xy.y.value - ); + 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); @@ -415,9 +409,8 @@ std::optional GridBase::findFullCoordinates(const Point& point, if (abs(deviation) > 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(), + "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, diff --git a/src/cpp/fs/Model.cpp b/src/cpp/fs/Model.cpp index 88323c605e..b9f4992373 100644 --- a/src/cpp/fs/Model.cpp +++ b/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> Model::runIterations( static_assert(std::numeric_limits::digits10 >= precision); const auto lat = static_cast(abs(start_point.latitude()) * pow(10, precision)); const auto lon = static_cast(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(start_day); diff --git a/src/cpp/fs/Point.h b/src/cpp/fs/Point.h index 88fc2195fc..25e1ff807f 100644 --- a/src/cpp/fs/Point.h +++ b/src/cpp/fs/Point.h @@ -45,4 +45,25 @@ class Point MathSize longitude_; }; } +template +concept LatLong = requires(const T& p) { + { p.latitude() } noexcept -> std::same_as; + { p.longitude() } noexcept -> std::same_as; +}; +template + requires LatLong +struct std::formatter : std::formatter +{ + constexpr auto parse(std::format_parse_context& ctx) { return ctx.begin(); } + // apple clang didn't work with std::format_context + template + 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::format(tmp, ctx); + } +}; #endif diff --git a/src/cpp/fs/SimpleFBP.cpp b/src/cpp/fs/SimpleFBP.cpp index 78f31bbc3d..b2f3b4bb40 100644 --- a/src/cpp/fs/SimpleFBP.cpp +++ b/src/cpp/fs/SimpleFBP.cpp @@ -532,10 +532,9 @@ vector 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 ); diff --git a/src/cpp/fs/projection.cpp b/src/cpp/fs/projection.cpp index 2a0993e7dd..7192f7a134 100644 --- a/src/cpp/fs/projection.cpp +++ b/src/cpp/fs/projection.cpp @@ -49,9 +49,8 @@ std::optional to_proj4( *y = b.enu.n; c = proj_trans(P, PJ_INV, b); cout << std::format( - "longitude: {:f}, latitude: {:f} => easting: {:.3f}, northing: {:.3f} => x: {:.3f}, y: {:.3f} => longitude: {:g}, latitude: {:g}\n", - point.longitude(), - point.latitude(), + "{} => easting: {:.3f}, northing: {:.3f} => x: {:.3f}, y: {:.3f} => longitude: {:g}, latitude: {:g}\n", + point, b.enu.e, b.enu.n, b.xy.x, @@ -113,9 +112,8 @@ fs::XYPos from_lat_long(const string_view proj4, const fs::Point& point) // #ifdef DEBUG_PROJ PJ_COORD c = proj_trans(P, PJ_INV, b); fs::logging::verbose( - "longitude: {:f}, latitude: {:f} => easting: {:.3f}, northing: {:.3f} => x: {:.3f}, y: {:.3f} => longitude: {:g}, latitude: {:g}", - point.longitude(), - point.latitude(), + "{} => easting: {:.3f}, northing: {:.3f} => x: {:.3f}, y: {:.3f} => longitude: {:g}, latitude: {:g}", + point, b.enu.e, b.enu.n, b.xy.x, From 43e13b61052f72e7f69803e98f513b5a42915380 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 13:33:36 +0000 Subject: [PATCH 06/11] use atan2 for direction check --- src/cpp/fs/Grid.cpp | 17 +++-------------- src/cpp/fs/projection.cpp | 17 +++++++++++++++++ src/cpp/fs/projection.h | 1 + 3 files changed, 21 insertions(+), 14 deletions(-) diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index c5885e0c28..a27e1ffaa8 100644 --- a/src/cpp/fs/Grid.cpp +++ b/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,29 +392,17 @@ std::optional GridBase::findCoordinates(const Point& point, const b std::optional 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; + auto deviation = find_north_south_deviation(this->proj4_, point).value; if (abs(deviation) > 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", + "Due north is not the top of the raster for {} with proj4 '{:s}' 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 ); diff --git a/src/cpp/fs/projection.cpp b/src/cpp/fs/projection.cpp index 7192f7a134..f6c79f7709 100644 --- a/src/cpp/fs/projection.cpp +++ b/src/cpp/fs/projection.cpp @@ -3,6 +3,7 @@ #include #include "Log.h" #include "Point.h" +#include "Radians.h" #include "Settings.h" #include "unstable.h" #include "Util.h" @@ -191,4 +192,20 @@ 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{10}; + 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); + // angle is going to be how far off North we are + auto deviation = + Radians{atan2(grid1.y.value - grid0.y.value, grid1.x.value - grid0.x.value)}.asDegrees(); + logging::note( + "Deviation for {} to {} with proj4 '{:s}' - is {:f} degrees", p0, p1, proj4, deviation.value + ); + return deviation; +}; } diff --git a/src/cpp/fs/projection.h b/src/cpp/fs/projection.h index ee4416fc14..53d77b3943 100644 --- a/src/cpp/fs/projection.h +++ b/src/cpp/fs/projection.h @@ -36,5 +36,6 @@ 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); } #endif From 9db9ecf75bed743cb11121eb770247be2b46fff5 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 14:07:10 +0000 Subject: [PATCH 07/11] show points used for deviation calculation --- src/cpp/fs/Grid.cpp | 2 +- src/cpp/fs/projection.cpp | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index a27e1ffaa8..3a658e53a3 100644 --- a/src/cpp/fs/Grid.cpp +++ b/src/cpp/fs/Grid.cpp @@ -395,7 +395,7 @@ std::optional GridBase::findFullCoordinates(const Point& point, 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; + constexpr MathSize MAX_DEVIATION = 0.1; auto deviation = find_north_south_deviation(this->proj4_, point).value; if (abs(deviation) > MAX_DEVIATION) { diff --git a/src/cpp/fs/projection.cpp b/src/cpp/fs/projection.cpp index f6c79f7709..a6629db345 100644 --- a/src/cpp/fs/projection.cpp +++ b/src/cpp/fs/projection.cpp @@ -195,17 +195,19 @@ string try_fix_meridian(const string_view proj4) Degrees find_north_south_deviation(const string_view proj4, const Point& p0) { constexpr auto MAX_LATITUDE{89.9}; - constexpr auto LATITUDE_DISTANCE{10}; + 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::note("{} => ({:f}, {:f})", p0, grid0.x.value, grid0.y.value); + logging::note("{} => ({:f}, {:f})", p1, grid1.x.value, grid1.y.value); // angle is going to be how far off North we are auto deviation = Radians{atan2(grid1.y.value - grid0.y.value, grid1.x.value - grid0.x.value)}.asDegrees(); logging::note( "Deviation for {} to {} with proj4 '{:s}' - is {:f} degrees", p0, p1, proj4, deviation.value ); - return deviation; + return deviation - Radians::D_090().asDegrees(); }; } From 4da99a37b870d44631f409e0ce8cd914877fe915 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 16:44:35 +0000 Subject: [PATCH 08/11] correct deviation from math direction --- src/cpp/fs/projection.cpp | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/cpp/fs/projection.cpp b/src/cpp/fs/projection.cpp index a6629db345..532b052971 100644 --- a/src/cpp/fs/projection.cpp +++ b/src/cpp/fs/projection.cpp @@ -200,14 +200,23 @@ Degrees find_north_south_deviation(const string_view proj4, const Point& p0) Point p1{lat1, p0.longitude()}; auto grid0 = from_lat_long(proj4, p0); auto grid1 = from_lat_long(proj4, p1); - logging::note("{} => ({:f}, {:f})", p0, grid0.x.value, grid0.y.value); - logging::note("{} => ({:f}, {:f})", p1, grid1.x.value, grid1.y.value); + 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{atan2(grid1.y.value - grid0.y.value, grid1.x.value - grid0.x.value)}.asDegrees(); - logging::note( + Radians{atan2(grid1.y.value - grid0.y.value, grid1.x.value - grid0.x.value)}.asDegrees() + - Radians::D_090().asDegrees(); + logging::debug( "Deviation for {} to {} with proj4 '{:s}' - is {:f} degrees", p0, p1, proj4, deviation.value ); - return deviation - Radians::D_090().asDegrees(); + return deviation; }; } From 2d0d4a8fa2b15f4433ccf2a79fa3e150c9180f59 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 18:39:03 +0000 Subject: [PATCH 09/11] improve due north detection --- src/cpp/fs/Grid.cpp | 28 ++++++++++-------------- src/cpp/fs/projection.cpp | 46 +++++++++++++++++++++++++++++++++++++-- src/cpp/fs/projection.h | 6 +++++ 3 files changed, 62 insertions(+), 18 deletions(-) diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index 3a658e53a3..276c60c4a6 100644 --- a/src/cpp/fs/Grid.cpp +++ b/src/cpp/fs/Grid.cpp @@ -395,26 +395,22 @@ std::optional GridBase::findFullCoordinates(const Point& point, 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.1; - auto deviation = find_north_south_deviation(this->proj4_, point).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}' gives deviation of {:f} degrees which exceeds maximum of {:f} degrees", - point, - this->proj4_, - 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 {}; + } + const auto upper_right = to_lat_long(this->proj4_, XYPos{XPos{xurcorner_}, YPos{yurcorner_}}); + if (!check_deviation("upper right", proj4_, upper_right, MAX_DEVIATION_EDGES)) + { + return {}; } logging::verbose("Lower left is ({:f}, {:f})", this->xllcorner_, this->yllcorner_); // convert coordinates into cell position diff --git a/src/cpp/fs/projection.cpp b/src/cpp/fs/projection.cpp index 532b052971..5cfd6eb358 100644 --- a/src/cpp/fs/projection.cpp +++ b/src/cpp/fs/projection.cpp @@ -214,9 +214,51 @@ Degrees find_north_south_deviation(const string_view proj4, const Point& p0) auto deviation = Radians{atan2(grid1.y.value - grid0.y.value, grid1.x.value - grid0.x.value)}.asDegrees() - Radians::D_090().asDegrees(); - logging::debug( - "Deviation for {} to {} with proj4 '{:s}' - is {:f} degrees", p0, p1, proj4, deviation.value + 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( + "Due North of {:s} {} deviates from grid 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( + "Due North of {:s} {} deviates from grid North by {:g} degrees which is near maximum of {:f} degrees", + what, + p, + deviation, + max_deviation + ); + } + else + { + logging::info( + "Due North of {:s} {} deviates from grid North by {:g} degrees", + what, + p, + deviation, + max_deviation + ); + } + return true; +}; } diff --git a/src/cpp/fs/projection.h b/src/cpp/fs/projection.h index 53d77b3943..b6382806ae 100644 --- a/src/cpp/fs/projection.h +++ b/src/cpp/fs/projection.h @@ -37,5 +37,11 @@ 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 From d908b19d50f3a2eb58054bcff4ca074c5a2f975b Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 18:43:19 +0000 Subject: [PATCH 10/11] reword error message --- src/cpp/fs/projection.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cpp/fs/projection.cpp b/src/cpp/fs/projection.cpp index 5cfd6eb358..a1677fa7bd 100644 --- a/src/cpp/fs/projection.cpp +++ b/src/cpp/fs/projection.cpp @@ -212,8 +212,8 @@ Degrees find_north_south_deviation(const string_view proj4, const Point& p0) // angle is going to be how far off North we are // -90 to convert from math direction auto deviation = - Radians{atan2(grid1.y.value - grid0.y.value, grid1.x.value - grid0.x.value)}.asDegrees() - - Radians::D_090().asDegrees(); + 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 ); @@ -230,7 +230,7 @@ bool check_deviation( if (abs(deviation) > max_deviation) { logging::error( - "Due North of {:s} {} deviates from grid North by {:g} degrees which exceeds maximum of {:g} degrees", + "Grid North of {:s} {} deviates from true North by {:g} degrees which exceeds maximum of {:g} degrees", what, p, deviation, @@ -242,7 +242,7 @@ bool check_deviation( { // if close to max deviation then warn about it logging::warning( - "Due North of {:s} {} deviates from grid North by {:g} degrees which is near maximum of {:f} degrees", + "Grid North of {:s} {} deviates from true North by {:g} degrees which is near maximum of {:f} degrees", what, p, deviation, @@ -252,7 +252,7 @@ bool check_deviation( else { logging::info( - "Due North of {:s} {} deviates from grid North by {:g} degrees", + "Grid North of {:s} {} deviates from true North by {:g} degrees", what, p, deviation, From ce833e141e94b428e50b6cb51b8603c57b2b7463 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Mon, 22 Jun 2026 12:10:02 +0000 Subject: [PATCH 11/11] check deviation of lower right instead of upper right --- src/cpp/fs/Grid.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index 276c60c4a6..295407765b 100644 --- a/src/cpp/fs/Grid.cpp +++ b/src/cpp/fs/Grid.cpp @@ -407,8 +407,9 @@ std::optional GridBase::findFullCoordinates(const Point& point, { return {}; } - const auto upper_right = to_lat_long(this->proj4_, XYPos{XPos{xurcorner_}, YPos{yurcorner_}}); - if (!check_deviation("upper right", proj4_, upper_right, MAX_DEVIATION_EDGES)) + // 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 {}; }