From 12cf96a658c9b56aa1d9e3c297db4f3e5bf2a536 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 13:33:36 +0000 Subject: [PATCH 1/6] 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 8a719866bc10b8415f625d25fcab821308d6a02b Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 14:07:10 +0000 Subject: [PATCH 2/6] 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 c287c14e7c6216b7d11ba3a792a110d7e51f04b9 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 16:44:35 +0000 Subject: [PATCH 3/6] 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 0a0dcdb5447417f11356f2bb313d48f4d1ed628e Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 18:39:03 +0000 Subject: [PATCH 4/6] 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 0fe8130f19b038232c6fb673a55f8bff80b93b4b Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 18 Jun 2026 18:43:19 +0000 Subject: [PATCH 5/6] 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 7320671466baf7d9e0a309f109fceae2df3d3329 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Mon, 22 Jun 2026 12:10:02 +0000 Subject: [PATCH 6/6] 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 {}; }