feat(core): CrsTransform(PROJ 封装, UTM/WGS84/WebMercator 互转)

- RAII 管理 PJ_CONTEXT/PJ;normalize_for_visualization 统一轴序为 (x=经度/东, y=纬度/北)
- vcpkg 加 proj 依赖;core 链接 PROJ::proj(保持 core 纯净,无 Qt/VTK)
- 测试经 CMake gtest_discover_tests 注入 PROJ_DATA,ctest 开箱即用
- 修正用例期望值:UTM49N(516868) 实际经度约 111.16°E(中央经线 111°E);
  WebMercator tfw 原点 114.16°E 在 49N 的 easting 约 825km,均与 PROJ 数据库一致
This commit is contained in:
gaozheng 2026-06-07 20:06:37 +08:00
parent e5a48c5af7
commit 4fdc6f7b86
6 changed files with 102 additions and 3 deletions

View File

@ -1,13 +1,15 @@
find_package(Eigen3 CONFIG REQUIRED)
find_package(PROJ CONFIG REQUIRED)
add_library(geopro_core STATIC
geo/LocalFrame.cpp
geo/CrsTransform.cpp
model/ColorScale.cpp
algo/IdwInterpolator.cpp
)
target_include_directories(geopro_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
target_link_libraries(geopro_core PUBLIC Eigen3::Eigen)
target_link_libraries(geopro_core PUBLIC Eigen3::Eigen PROJ::proj)
target_compile_features(geopro_core PUBLIC cxx_std_17)
# core Qt / VTK

View File

@ -0,0 +1,37 @@
#include "geo/CrsTransform.hpp"
#include <proj.h>
#include <stdexcept>
namespace geopro::core {
struct CrsTransform::Impl {
PJ_CONTEXT* ctx = nullptr;
PJ* pj = nullptr;
};
CrsTransform::CrsTransform(const std::string& src, const std::string& dst)
: impl_(std::make_unique<Impl>()) {
impl_->ctx = proj_context_create();
impl_->pj = proj_create_crs_to_crs(impl_->ctx, src.c_str(), dst.c_str(), nullptr);
if (!impl_->pj) throw std::runtime_error("CrsTransform: failed to create " + src + "->" + dst);
PJ* norm = proj_normalize_for_visualization(impl_->ctx, impl_->pj);
if (norm) { proj_destroy(impl_->pj); impl_->pj = norm; }
}
CrsTransform::~CrsTransform() {
if (impl_->pj) proj_destroy(impl_->pj);
if (impl_->ctx) proj_context_destroy(impl_->ctx);
}
Xy CrsTransform::forward(double x, double y) const {
PJ_COORD c = proj_coord(x, y, 0, 0);
PJ_COORD r = proj_trans(impl_->pj, PJ_FWD, c);
return Xy{r.xy.x, r.xy.y};
}
Xy CrsTransform::inverse(double x, double y) const {
PJ_COORD c = proj_coord(x, y, 0, 0);
PJ_COORD r = proj_trans(impl_->pj, PJ_INV, c);
return Xy{r.xy.x, r.xy.y};
}
} // namespace geopro::core

View File

@ -0,0 +1,20 @@
#pragma once
#include <string>
#include <memory>
namespace geopro::core {
struct Xy { double x, y; };
// PROJ 封装:源 CRS <-> 目标 CRS。用于世界系<->GIS<->经纬度、影像异源 CRS 重投影(设计 §5)。
class CrsTransform {
public:
CrsTransform(const std::string& srcCrs, const std::string& dstCrs);
~CrsTransform();
Xy forward(double x, double y) const; // src -> dst
Xy inverse(double x, double y) const; // dst -> src
private:
struct Impl;
std::unique_ptr<Impl> impl_;
};
} // namespace geopro::core

View File

@ -7,12 +7,25 @@ add_executable(geopro_tests smoke_test.cpp)
target_link_libraries(geopro_tests PRIVATE GTest::gtest GTest::gtest_main)
include(GoogleTest)
gtest_discover_tests(geopro_tests)
# PROJ proj.db vcpkg share/proj
# PROJ_DATA 使 `ctest -R CrsTransform` set
file(GLOB _proj_data_dirs
"${CMAKE_BINARY_DIR}/vcpkg_installed/*/share/proj"
)
if(_proj_data_dirs)
list(GET _proj_data_dirs 0 GEOPRO_PROJ_DATA)
gtest_discover_tests(geopro_tests
PROPERTIES ENVIRONMENT "PROJ_DATA=${GEOPRO_PROJ_DATA}")
else()
gtest_discover_tests(geopro_tests)
endif()
target_sources(geopro_tests PRIVATE core/test_local_frame.cpp)
target_sources(geopro_tests PRIVATE core/test_model.cpp)
target_sources(geopro_tests PRIVATE core/test_color_scale.cpp)
target_sources(geopro_tests PRIVATE core/test_idw.cpp)
target_sources(geopro_tests PRIVATE core/test_crs_transform.cpp)
target_link_libraries(geopro_tests PRIVATE geopro_core)
add_subdirectory(spike) # spike S3: banded contour

View File

@ -0,0 +1,26 @@
#include <gtest/gtest.h>
#include "geo/CrsTransform.hpp"
using namespace geopro::core;
// UTM 49N (CM 111°E, false easting 500000) <-> WGS84 经纬度往返。
// 输入 easting≈516868 仅比中央经线偏东 ~16.9km,故经度落在 ~111.16°E非 114°
// normalize_for_visualization 后forward 输出 (lon, lat)inverse 输入 (lon, lat)。
TEST(CrsTransform, Utm49nToWgs84RoundTrip) {
CrsTransform t("EPSG:32649", "EPSG:4326");
auto ll = t.forward(516868.0, 2494259.0); // (east,north) -> (lon,lat)
EXPECT_NEAR(ll.x, 111.16, 0.05); // lon
EXPECT_NEAR(ll.y, 22.55, 0.05); // lat
auto en = t.inverse(ll.x, ll.y); // (lon,lat) -> (east,north)
EXPECT_NEAR(en.x, 516868.0, 1.0);
EXPECT_NEAR(en.y, 2494259.0, 1.0);
}
// WebMercator -> UTM 49N 异源重投影(设计 §5 影像重投影)。
// tfw 原点经度真实约 114.16°E距 49N 中央经线 111°E 偏东 ~3.16°,
// 故其 49N easting 远大于 false easting落在 ~825km 处(合理且 < UTM 带宽上限)。
TEST(CrsTransform, WebMercatorToUtm) {
CrsTransform t("EPSG:3857", "EPSG:32649");
auto p = t.forward(12708343.88, 2577685.90); // tfw 原点
EXPECT_GT(p.x, 800000.0); EXPECT_LT(p.x, 850000.0); // UTM 49N 合理 Easting
EXPECT_GT(p.y, 2400000.0); EXPECT_LT(p.y, 2600000.0); // 合理 Northing
}

View File

@ -5,6 +5,7 @@
"dependencies": [
"eigen3",
"gtest",
"nlohmann-json"
"nlohmann-json",
"proj"
]
}