Projections.h

#pragma once
 
#include <iostream>
#include <memory>
 
#include "Constants.h"
#include "Ellipsoid.h"
#include "Vector.h"
 
#include <isce3/except/Error.h>
 
namespace isce3 { namespace core {

class ProjectionBase {
    Ellipsoid _ellipse;

    int _epsgcode;
 
public:
    ProjectionBase(int code)
        : _ellipse(), _epsgcode(code)
    {}

    int code() const { return _epsgcode; }

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

    virtual void print() const = 0;

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

    Vec3 forward(const Vec3& llh) const
    {
        Vec3 xyz;
        const int ret = forward(llh, xyz);
        if (ret != 0) {
            throw isce3::except::RuntimeError(ISCE_SRCINFO(),
                    "Forward projection transformation failed");
        }
        return xyz;
    }

    virtual int inverse(const Vec3& xyz, Vec3& llh) const = 0;

    Vec3 inverse(const Vec3& xyz) const
    {
        Vec3 llh;
        const int ret = inverse(xyz, llh);
        if (ret != 0) {
            throw isce3::except::RuntimeError(ISCE_SRCINFO(),
                    "Inverse projection transformation failed");
        }
        return llh;
    }
 
    virtual ~ProjectionBase() = default;
};

class LonLat : public ProjectionBase {
public:
    LonLat() : ProjectionBase(4326) {}

    void print() const override;
 
    int forward(const Vec3&, Vec3&) const override;
 
    int inverse(const Vec3&, Vec3&) const override;
};
 
inline void LonLat::print() const
{
    std::cout << "Projection: LatLon" << std::endl
              << "EPSG: " << code() << std::endl;
}
 
inline int LonLat::forward(const Vec3& in, Vec3& 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 Vec3& in, Vec3& 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:
    Geocent() : ProjectionBase(4978) {}

    void print() const;

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

    int inverse(const Vec3& xyz, Vec3& 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:
    UTM(int);

    void print() const override;

    int forward(const Vec3& llh, Vec3& xyz) const override;

    int inverse(const Vec3& xyz, Vec3& llh) const override;
};
 
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:
    PolarStereo(int);

    void print() const override;

    int forward(const Vec3&, Vec3&) const override;

    int inverse(const Vec3&, Vec3&) const override;
};
 
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:
    CEA();
 
    void print() const override;

    int forward(const Vec3& llh, Vec3& xyz) const override;

    int inverse(const Vec3& xyz, Vec3& llh) const override;
};
 
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);
 
}} // namespace isce3::core