diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index c5885e0c28..295407765b 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,41 +392,26 @@ 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; - 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 diff --git a/src/cpp/fs/projection.cpp b/src/cpp/fs/projection.cpp index 7192f7a134..a1677fa7bd 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,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; +}; } diff --git a/src/cpp/fs/projection.h b/src/cpp/fs/projection.h index ee4416fc14..b6382806ae 100644 --- a/src/cpp/fs/projection.h +++ b/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