proj6及后续版本投影转换

#ifndef COORDINATE_PROJ_HPP
#define COORDINATE_PROJ_HPP

#include "proj.h"
#include <stdlib.h>
#include <stdio.h>
#include <exception>

class proj_conv
{
public:
    proj_conv() {
        C = nullptr;

        Pmerc = nullptr;
        P_for_merc = nullptr;

        Pecef = nullptr;
        P_for_ecef = nullptr;
    }

    inline static int get_zone(double lon) { return int((lon + 1.5)/3.0); }
    inline static int get_lon_0(double lon) { return get_zone(lon) * 3; }

    bool init(double lon) {
        try {

            release();

            //merc
            C = proj_context_create();
            char psztmp[32] = "";
            sprintf(psztmp, "+proj=tmerc +lon_0=%d +x_0=500000", get_lon_0(lon));
            Pmerc = proj_create_crs_to_crs(C, "EPSG:4326", "+proj=tmerc +lon_0=117 +x_0=500000", nullptr);
            if (nullptr == Pmerc) {
                fprintf(stderr, "Oops\n");
                return false;
            }
            /* This will ensure that the order of coordinates for the input CRS */
            /* will be longitude, latitude, whereas EPSG:4326 mandates latitude, */
            /* longitude */
            if( nullptr == (P_for_merc = proj_normalize_for_visualization(C, Pmerc))) {
                fprintf(stderr, "Oops\n");
                return false;
            }
            proj_destroy(Pmerc);
            Pmerc = P_for_merc;

            //ecef
            Pecef = proj_create_crs_to_crs(C, "EPSG:4326", "+proj=geocent +units=m", nullptr);
            if (nullptr == Pecef) {
                fprintf(stderr, "Oops\n");
                return false;
            }
            if( nullptr == (P_for_ecef = proj_normalize_for_visualization(C, Pecef)))  {
                fprintf(stderr, "Oops\n");
                return false;
            }
            proj_destroy(Pecef);
            Pecef = P_for_ecef;

            return true;
        } catch (std::exception& exc) {

        }
        return false;
    }

    bool blhtoenu(double b, double l, double h, double& e, double& n, double& u){
        try {

            /* coordinates is longitude, latitude, and values are expressed in degrees. */
            /* 40.098393475 116.147123333 */
            incoor = proj_coord (l, b, h, 0);

            /* transform to UTM zone 32, then back to geographical */
            outcoor = proj_trans (Pmerc, PJ_FWD, incoor);

            e = outcoor.enu.e;
            n = outcoor.enu.n;
            u = outcoor.enu.u;
            return true;
        } catch (std::exception& exc) {

        }
        return false;
    }
    bool blhtoecef(double b, double l, double h, double& x, double& y, double& z){
        try {
            /* coordinates is longitude, latitude, and values are expressed in degrees. */
            /* 40.098393475 116.147123333 */
            incoor = proj_coord (l, b, h, 0);

            /* transform to UTM zone 32, then back to geographical */
            outcoor = proj_trans (Pecef, PJ_FWD, incoor);

            x = outcoor.xyz.x;
            y = outcoor.xyz.y;
            z = outcoor.xyz.z;
            return true;
        } catch (std::exception& exc) {

        }
        return false;
    }
    bool eceftoblh(double x, double y, double z, double& b, double& l, double& h){
        try {
            /* coordinates is longitude, latitude, and values are expressed in degrees. */
            /* 40.098393475 116.147123333 */
            incoor = proj_coord (x, y, z, 0);

            /* transform to UTM zone 32, then back to geographical */
            outcoor = proj_trans (Pecef, PJ_INV, incoor);

            l = outcoor.lpz.lam;
            b = outcoor.lpz.phi;
            h = outcoor.lpz.z;
            return true;
        } catch (std::exception& exc) {

        }
        return false;
    }

    void release() {
        /* Clean up */
        if (Pmerc)
            proj_destroy (Pmerc);

        if (Pecef)
            proj_destroy (Pecef);

        if (C)
            proj_context_destroy (C); /* may be omitted in the single threaded case */
    }

    void run(){

        /* coordinates is longitude, latitude, and values are expressed in degrees. */
        /* 40.098393475 116.147123333 */
        double b = 40.098393475, l = 116.147123333, h = 10;

        init(l);

        double e,n,u,x,y,z;
        blhtoenu(b,l,h,e,n,u);
        blhtoecef(b,l,h,x,y,z);
        eceftoblh(x,y,z,b,l,h);

        printf("b:%.9f\r\nl:%.9f\r\nh:%.3f\r\ne:%.3f\r\nn:%.3f\r\nu:%.3f\r\nx:%.3f\r\ny:%.3f\r\nz:%.3f\r\n",
               b,l,h,e,n,u,x,y,z);

        release();
    }

public:
    int zone, lon_0;
    PJ_CONTEXT *C;

    PJ *Pmerc;
    PJ* P_for_merc;

    PJ *Pecef;
    PJ* P_for_ecef;

    PJ_COORD incoor, outcoor;
};

#endif // COORDINATE_PROJ_HPP