Refactor/projection

src/cpp/fs/Grid.cpp

@@ -5,8 +5,8 @@
 #include <tiffio.h>
 #include <xtiffio.h>
 #include "Log.h"
+#include "projection.h"
 #include "tiff.h"
-#include "UTM.h"
 using fs::Idx;
 namespace fs
 {
@@ -391,34 +391,36 @@ 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);
+  // 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<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);
+  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(
       "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,
+      grid_north.x.value,
+      grid_south.x.value,
       deviation,
       MAX_DEVIATION
     );
@@ -434,10 +436,10 @@ std::optional<FullCoordinates> 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<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/UTM.cpp → src/cpp/fs/projection.cpp

@@ -1,5 +1,5 @@
 /* SPDX-License-Identifier: AGPL-3.0-or-later */
-#include "UTM.h"
+#include "projection.h"
 #include <proj.h>
 #include "Log.h"
 #include "Point.h"
@@ -8,6 +8,62 @@
 #include "Util.h"
 namespace fs
 {
+std::optional<FullCoordinates> 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)

src/cpp/fs/UTM.h → 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<FullCoordinates> 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