coordinate_system.h

/* The MIT License (MIT)
 *
 * (c) Jürgen Simon 2014 (juergen.simon@uni-bonn.de)
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef RADOLAN_COORDINATE_SYSTEM_H
#define RADOLAN_COORDINATE_SYSTEM_H

//C-headers
//C++ headers
//3rd-party headers
#include <radolan/types.h>

// TODO: refactor using RDGridPoint and RDGeographicalPointRad etc

namespace Radolan {

    typedef struct
    {
        double latitude;
        double longitude;
    } RDGeographicalPointRad;

    inline RDGeographicalPointRad rdGeographicalPointRad(double longitude, double latitude) {
        RDGeographicalPointRad p;
        p.longitude = longitude;
        p.latitude = latitude;
        return p;
    }

    inline bool RDEquals(RDGeographicalPointRad a, RDGeographicalPointRad b) {
        return a.latitude == b.latitude && a.longitude == b.longitude;
    }

#define RDGeographicalPointRadUndefined rdGeographicalPointRad(NAN,NAN)

    typedef struct
    {
        double latitude;
        double longitude;
    } RDGeographicalPoint;

    inline RDGeographicalPoint rdGeographicalPoint(double longitude, double latitude) {
        RDGeographicalPoint p;
        p.longitude = longitude;
        p.latitude = latitude;
        return p;
    }

    inline bool RDEquals(RDGeographicalPoint a, RDGeographicalPoint b) {
        return a.latitude == b.latitude && a.longitude == b.longitude;
    }

#define RDGeographicalPointUndefined rdGeographicalPoint(NAN,NAN)

    typedef struct
    {
        double x;
        double y;
    } RDCartesianPoint;

    inline RDCartesianPoint rdCartesianPoint(double x, double y) {
        RDCartesianPoint p;
        p.x = x;
        p.y = y;
        return p;
    }

    inline bool RDEquals(RDCartesianPoint a, RDCartesianPoint b) {
        return a.x == b.x && a.y == b.y;
    }

#define RDCartesianPointUndefined rdCartesianPoint(NAN,NAN)

    typedef struct
    {
        int ix;
        int iy;
    } RDGridPoint;

    inline RDGridPoint rdGridPoint(int ix, int iy) {
        RDGridPoint p;
        p.ix = ix;
        p.iy = iy;
        return p;
    }

    inline bool RDEquals(RDGridPoint a, RDGridPoint b) {
        return a.ix == b.ix && a.iy == b.iy;
    }

#define RDGridPointUndefined rdGridPoint(-1, -1)

    inline RDDataType RDValueAt(RDScan *scan, RDGridPoint gridPoint) {
        return scan->data[gridPoint.iy * scan->dimLon + gridPoint.ix];
    }

    //  1 | 2
    //  - + -
    //  0 | 3

    typedef enum
    {
        RDLowerLeft = 0,
        RDUpperLeft = 1,
        RDUpperRight = 2,
        RDLowerRight = 3
    } RDGridQuadrant;

    class RDCoordinateSystem
    {
    public:

        RDCoordinateSystem(RDScanType type);

        virtual ~RDCoordinateSystem();

        void setScanType(RDScanType type);

        RDCartesianPoint cartesianCoordinate(RDGeographicalPointRad coord);

        RDCartesianPoint cartesianCoordinate(RDGeographicalPoint coord);

        RDCartesianPoint cartesianCoordinate(RDGridPoint gridpoint);

        RDGeographicalPoint geographicalCoordinate(RDCartesianPoint p);

        RDGeographicalPoint geographicalCoordinate(RDGridPoint p);

        RDGeographicalPointRad geographicalCoordinateRad(RDCartesianPoint p);

        RDGeographicalPointRad geographicalCoordinateRad(RDGridPoint p);

        RDGeographicalPointRad toRad(RDGeographicalPoint p);

        RDGeographicalPoint toDeg(RDGeographicalPointRad p);

        RDGeographicalPoint *geographicalPolygonForGridpoint(RDGridPoint p, int &num);

        RDGridPoint gridPoint(RDCartesianPoint p, bool &isInside);

        RDGridPoint gridPoint(RDGeographicalPoint p, bool &isInside);

        RDGridPoint gridPoint(RDGeographicalPointRad p, bool &isInside);

        RDGridQuadrant RDQuadrant(RDGridPoint p);

        double polarStereographicScalingFactor(double phi0, double phi);

    private:

        // these parameters depend on the radolan product
        int m_radolanGridCountHorizontal;
        int m_radolanGridCountVertical;

        // Origin of the coordinate system in geographical coordinates
        RDGeographicalPoint m_originGeographical;

        // Origin in cartesian coordinates
        RDCartesianPoint m_originCartesian;

        // Offset of the origin from lower left corner (example: EX: -600, -800)
        RDCartesianPoint m_offset;

        RDScanType m_scanType;

        // These are parameters of the projection itself and don't change
        static const double LAMBDA_0;
        static const double PHI_0;
        static const double R_EARTH;
        static const double MESH_WIDTH;
        static const double RADOLAN_GRID_COUNT;
        static const double LAMBDA_MIN_BOTTOM;
        static const double PHI_MIN_BOTTOM;
        static const double LAMBDA_MAX_TOP;
        static const double PHI_MAX_TOP;

        double deg(double rad);

        double rad(double deg);

        // Called after changing scan type to update internals
        void updateGridInfo();

    }; //class RDCoordinateSystem
}

#endif /* Header Guard */