From b0081ce2ec110a0000aa121ee6e8c1aa24427b14 Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Thu, 11 Jun 2026 14:50:49 +0000 Subject: [PATCH 1/2] combine UTM/project.h/cpp into projection.h/cpp --- src/cpp/fs/project.cpp | 62 -------------------------- src/cpp/fs/project.h | 17 ------- src/cpp/fs/{UTM.cpp => projection.cpp} | 58 +++++++++++++++++++++++- src/cpp/fs/{UTM.h => projection.h} | 13 +++++- 4 files changed, 68 insertions(+), 82 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} (74%) rename src/cpp/fs/{UTM.h => projection.h} (80%) diff --git a/src/cpp/fs/project.cpp b/src/cpp/fs/project.cpp deleted file mode 100644 index 719bc7243ac..00000000000 --- 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 f6699d1f336..00000000000 --- 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 74% rename from src/cpp/fs/UTM.cpp rename to src/cpp/fs/projection.cpp index 7042af42780..f299c84d5c2 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 diff --git a/src/cpp/fs/UTM.h b/src/cpp/fs/projection.h similarity index 80% rename from src/cpp/fs/UTM.h rename to src/cpp/fs/projection.h index 6fb60b4b166..2356366fb47 100644 --- a/src/cpp/fs/UTM.h +++ b/src/cpp/fs/projection.h @@ -1,9 +1,18 @@ /* 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 "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 From d72615c0144cc8852677062507afcd33fedfad0e Mon Sep 17 00:00:00 2001 From: Jordan Evens Date: Wed, 17 Jun 2026 13:20:58 +0000 Subject: [PATCH 2/2] return Point/XYPos from functions instead of passing x/y pointers --- src/cpp/fs/Grid.cpp | 44 ++++++++++++++++++++------------------- src/cpp/fs/projection.cpp | 12 +++++------ src/cpp/fs/projection.h | 5 +++-- 3 files changed, 32 insertions(+), 29 deletions(-) diff --git a/src/cpp/fs/Grid.cpp b/src/cpp/fs/Grid.cpp index ea9aa2cddc5..70176823f2f 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/projection.cpp b/src/cpp/fs/projection.cpp index f299c84d5c2..2a0993e7dd8 100644 --- a/src/cpp/fs/projection.cpp +++ b/src/cpp/fs/projection.cpp @@ -97,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 @@ -109,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( @@ -127,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 @@ -137,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/projection.h b/src/cpp/fs/projection.h index 2356366fb47..ee4416fc143 100644 --- a/src/cpp/fs/projection.h +++ b/src/cpp/fs/projection.h @@ -2,6 +2,7 @@ #ifndef FS_PROJECTION_H #define FS_PROJECTION_H #include "stdafx.h" +#include "Location.h" #include "Point.h" namespace fs { @@ -32,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