RasterProcessTool/RasterMainWidgetGUI/RasterMainWidget/crs.cpp

//#include <QDebug>
//#include <proj_api.h>
//
//#include <crs.h>
//#pragma execution_character_set("utf-8")
//
//namespace LAMPMainWidget {
//
//PointXY
//CRS::forward(const PointXY &point) const {
//  auto pjCtx = pj_ctx_alloc();
//  auto pjLonlat = pj_init_plus_ctx(pjCtx, "+proj=longlat +datum=WGS84 +no_defs");
//  if (!pjLonlat) {
//    qWarning() << "初始化wgs84坐标系失败=>" << pj_strerrno(pj_ctx_get_errno(pjCtx));
//    pj_ctx_free(pjCtx);
//    return PointXY{};
//  }
//  auto pjDest = pj_init_plus_ctx(pjCtx, proj4Def().toStdString().c_str());
//  if (!pjDest) {
//    qWarning() << "初始化目标坐标系失败=>" << pj_strerrno(pj_ctx_get_errno(pjCtx));
//    pj_free(pjLonlat);
//    pj_ctx_free(pjCtx);
//    return PointXY{};
//  }
//
//  double x{point.x() * DEG_TO_RAD}, y{point.y() * DEG_TO_RAD};
//  auto err = pj_transform(pjLonlat, pjDest, 1, 1, &x, &y, nullptr);
//  if (err) {
//    qWarning() << "坐标转换失败=>" << point;
//    pj_free(pjDest);
//    pj_free(pjLonlat);
//    pj_ctx_free(pjCtx);
//    return PointXY{};
//  }
//
//  pj_free(pjDest);
//  pj_free(pjLonlat);
//  pj_ctx_free(pjCtx);
//
//  return PointXY{x, y};
//}
//
//PointXY
//CRS::inverse(const PointXY &point) const {
//  auto pjCtx = pj_ctx_alloc();
//  auto pjLonlat = pj_init_plus_ctx(pjCtx, "+proj=longlat +datum=WGS84 +no_defs");
//  if (!pjLonlat) {
//    qWarning() << "初始化wgs84坐标系失败=>" << pj_strerrno(pj_ctx_get_errno(pjCtx));
//    pj_ctx_free(pjCtx);
//    return PointXY{};
//  }
//  auto pjDest = pj_init_plus_ctx(pjCtx, proj4Def().toStdString().c_str());
//  if (!pjDest) {
//    qWarning() << "初始化目标坐标系失败=>" << pj_strerrno(pj_ctx_get_errno(pjCtx));
//    pj_free(pjLonlat);
//    pj_ctx_free(pjCtx);
//    return PointXY{};
//  }
//
//  double x{point.x()}, y{point.y()};
//  auto err = pj_transform(pjDest, pjLonlat, 1, 1, &x, &y, nullptr);
//  if (err) {
//    qWarning() << "坐标转换失败=>" << point;
//    pj_free(pjDest);
//    pj_free(pjLonlat);
//    pj_ctx_free(pjCtx);
//    return PointXY{};
//  }
//
//  pj_free(pjDest);
//  pj_free(pjLonlat);
//  pj_ctx_free(pjCtx);
//
//  return PointXY{x * RAD_TO_DEG, y * RAD_TO_DEG};
//}
//
//}



#include <QDebug>
#include <proj.h>

#include <crs.h>
#pragma execution_character_set("utf-8")
namespace LAMPMainWidget {

    PointXY
        CRS::forward(const PointXY& point) const {
        const double DEG_TO_RAD = 0.017453292519943295;
        PJ_CONTEXT* pjCtx = proj_context_create();
        if (!pjCtx) {
            qWarning() << u8"Create proj context failed";
            return PointXY{};
        }


       // PJ* target_crs = proj_create(pjCtx, proj4Def().toStdString().c_str());
       // if (!target_crs) {
       //     qWarning() << "Failed to create target CRS from PROJ string";
       //     return PointXY{};
       // }
       // 
       //
       //const char* wkt = proj_as_wkt(pjCtx, target_crs, PJ_WKT2_2019, nullptr);
       // if (wkt) {
       //     printf("CRS (WKT):\n%s\n", wkt);
       //     //return PointXY{};
       // }
        qWarning() << "create target_crs sucessfully";

        PJ* pjLonlat = proj_create_crs_to_crs(pjCtx,
            "EPSG:4326", // WGS84
            proj4Def().toStdString().c_str(),
            nullptr);

        if (!pjLonlat) {
            qWarning() << u8"convert unsucessfully from  EPSG:4326 to target project coodination =>" << proj_errno_string(proj_context_errno(pjCtx));
            //proj_destroy(target_crs);
            proj_context_destroy(pjCtx);
            return PointXY{};
        }

        PJ_COORD coordIn, coordOut;
        coordIn.lpzt.lam = point.x() * DEG_TO_RAD;
        coordIn.lpzt.phi = point.y() * DEG_TO_RAD;

        coordOut = proj_trans(pjLonlat, PJ_FWD, coordIn);
        if (proj_errno(pjLonlat)) {
            qWarning() << u8"convert unsucessfully point from  EPSG:4326 to point in EPSG:4326=>" << point;
            proj_destroy(pjLonlat);
            //proj_destroy(target_crs);
            proj_context_destroy(pjCtx);
            return PointXY{};
        }

        proj_destroy(pjLonlat);
        //proj_destroy(target_crs);
        proj_context_destroy(pjCtx);

        return PointXY{ coordOut.xy.x, coordOut.xy.y };
    }

    PointXY
        CRS::inverse(const PointXY& point) const {
        const double RAD_TO_DEG = 57.29577951308232;
        PJ_CONTEXT* pjCtx = proj_context_create();
        if (!pjCtx) {
            qWarning() << u8"create proj context ";
            return PointXY{};
        }

        //PJ* target_crs = proj_create(pjCtx, proj4Def().toStdString().c_str());
        //if (!target_crs) {
        //    qWarning() << "Failed to create target CRS from PROJ string";
        //    return PointXY{};
        //}

        //qWarning() << "create target_crs sucessfully";

        //const char* wkt = proj_as_wkt(pjCtx, target_crs, PJ_WKT2_2019, nullptr);
        //if (wkt) {
        //    printf("CRS (WKT):\n%s\n", wkt);
        //    return PointXY{};
        //}
        qWarning() << "create target_crs sucessfully";

        PJ* pjLonlat = proj_create_crs_to_crs(pjCtx,
            "EPSG:4326", // WGS84
            proj4Def().toStdString().c_str(),
            nullptr);

        if (!pjLonlat) {
            qWarning() << u8"=>" << proj_errno_string(proj_context_errno(pjCtx));
            //proj_destroy(target_crs);
            proj_context_destroy(pjCtx);
            return PointXY{};
        }

        PJ_COORD coordIn, coordOut;
        coordIn.xy.x = point.x();
        coordIn.xy.y = point.y();

        coordOut = proj_trans(pjLonlat, PJ_INV, coordIn);
        if (proj_errno(pjLonlat)) {
            qWarning() << u8"convert unsucessfully from  point in  target project coodination to point in EPSG:4326=>" << point;
            proj_destroy(pjLonlat);
            //proj_destroy(target_crs);
            proj_context_destroy(pjCtx);
            return PointXY{};
        }
        //proj_destroy(target_crs);
        proj_destroy(pjLonlat);
        proj_context_destroy(pjCtx);

        return PointXY{ coordOut.lpzt.lam * RAD_TO_DEG, coordOut.lpzt.phi * RAD_TO_DEG };
    }

}