Projections.h

//
// Source Author: Piyush Agram
// Co-Author: Joshua Cohen
// Copyright 2017
//

#pragma once

#include <iostream>
#include <memory>

#include "Constants.h"
#include "Ellipsoid.h"

namespace isce3 { namespace core {

    class ProjectionBase {
        Ellipsoid _ellipse;

        int _epsgcode;

    public:

        ProjectionBase(int code) : _ellipse(6378137.,.0066943799901), _epsgcode(code) {}

        inline int code() const { return _epsgcode; }

        inline const Ellipsoid& ellipsoid() const { return _ellipse; }

        virtual void print() const = 0;

        virtual int forward(const cartesian_t& llh, cartesian_t& xyz) const = 0 ;

        virtual int inverse(const cartesian_t& xyz, cartesian_t& llh) const = 0 ;
        inline Vec3 inverse(const Vec3& native) const {
            Vec3 llh;
            inverse(native, llh);
            return llh;
        }

        virtual ~ProjectionBase() {}
    };

    class LonLat : public ProjectionBase {
    public:
        // Value constructor
        LonLat() : ProjectionBase(4326) {}

        inline void print() const;
        // This will be a pass through for Lat/Lon
        inline int forward(const cartesian_t&, cartesian_t&) const;
        // This will also be a pass through for Lat/Lon
        inline int inverse(const cartesian_t&, cartesian_t&) const;
    };

    inline void LonLat::print() const {
        std::cout << "Projection: LatLon" << std::endl << "EPSG: " << code() << std::endl;
    }

    inline int LonLat::forward(const cartesian_t &in, cartesian_t &out) const {
        out[0] = in[0] * 180.0/M_PI;
        out[1] = in[1] * 180.0/M_PI;
        out[2] = in[2];
        return 0;
    }

    inline int LonLat::inverse(const cartesian_t &in, cartesian_t &out) const {
        out[0] = in[0] * M_PI/180.0;
        out[1] = in[1] * M_PI/180.0;
        out[2] = in[2];
        return 0;
    }

    class Geocent : public ProjectionBase {
    public:
        // Value constructor
        Geocent() : ProjectionBase(4978) {}

        inline void print() const;
        int forward(const cartesian_t& llh,cartesian_t& xyz) const;

        int inverse(const cartesian_t& xyz,cartesian_t& llh) const;
    };

    inline void Geocent::print() const {
        std::cout << "Projection: Geocent" << std::endl << "EPSG: " << code() << std::endl;
    }

    class UTM : public ProjectionBase {
        // Constants related to the projection system
        double lon0;
        int zone;
        bool isnorth;
        // Parameters from Proj.4
        double cgb[6], cbg[6], utg[6], gtu[6];
        double Qn, Zb;

    public:
        // Value constructor
        UTM(int);

        inline void print() const;
        int forward(const cartesian_t& llh, cartesian_t& xyz) const;

        int inverse(const cartesian_t& xyz, cartesian_t& llh) const;
    };

    inline void UTM::print() const {
        std::cout << "Projection: UTM" << std::endl
            << "Zone: " << zone << (isnorth ? "N" : "S") << std::endl
            << "EPSG: " << code() << std::endl;
    }

    class PolarStereo : public ProjectionBase {
        // Constants related to projection system
        double lat0, lon0, lat_ts, akm1, e;
        bool isnorth;

    public:
        // Value constructor
        PolarStereo(int);

        inline void print() const;
        int forward(const cartesian_t&,cartesian_t&) const;

        int inverse(const cartesian_t&,cartesian_t&) const;
    };

    inline void PolarStereo::print() const {
        std::cout << "Projection: " << (isnorth ? "North" : "South")
            << " Polar Stereographic" << std::endl
            << "EPSG: " << code() << std::endl;
    }

    class CEA: public ProjectionBase {
        // Constants related to projection system
        double apa[3];
        double lat_ts, k0, e, one_es, qp;

    public:
        // Value constructor
        CEA();

        inline void print() const;

        int forward(const cartesian_t& llh,cartesian_t& xyz) const;

        int inverse(const cartesian_t& xyz,cartesian_t& llh) const;
    };

    inline void CEA::print() const {
        std::cout << "Projection: Cylindrical Equal Area" << std::endl
            << "EPSG: " << code() << std::endl;
    }

    //This is to create a projection system from the EPSG code
    ProjectionBase* createProj(int epsg);

    inline
    std::unique_ptr<ProjectionBase> makeProjection(int epsg)
    {
        return std::unique_ptr<ProjectionBase>(createProj(epsg));
    }

    // This is to transform a point from one coordinate system to another
    int projTransform(ProjectionBase* in, ProjectionBase *out, const Vec3& inpts,
                      Vec3& outpts);
}}