gridpp.h

Go to the documentation of this file.

#ifndef GRIDPP_API_H
#define GRIDPP_API_H
#include <vector>
#include <string>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/index/rtree.hpp>
#ifdef _OPENMP
    #include <omp.h>
#endif

#define GRIDPP_VERSION "0.4.2"
#define __version__ GRIDPP_VERSION

// Preferred vector types
typedef std::vector<float> vec;
typedef std::vector<vec> vec2;
typedef std::vector<vec2> vec3;
typedef std::vector<int> ivec;
typedef std::vector<ivec> ivec2;

// Only use when double is required
typedef std::vector<double> dvec;
typedef std::vector<dvec> dvec2;

namespace gridpp {
    // Constants
    static const float MV = NAN;
    static const float MV_CML = -999;
    static const float pi = 3.14159265;
    static double radius_earth = 6.37e6;

    class KDTree;
    class Points;
    class Grid;
    class Nearest;
    class StructureFunction;
    class Transform;

    enum Extrapolation {
            OneToOne = 0,      
            MeanSlope = 10,    
            NearestSlope = 20, 
            Zero = 30,         
        };

    enum Statistic {
        Mean      = 0,  
        Min       = 10, 
        Median    = 20, 
        Max       = 30, 
        Quantile  = 40, 
        Std       = 50, 
        Variance  = 60, 
        Sum       = 70, 
        Unknown   = -1  
    };

    Statistic get_statistic(std::string name);

    std::string version();

    /****************************************
     * Data assimilation methods            *
     ****************************************/

    vec2 optimal_interpolation(const gridpp::Grid& bgrid,
            const vec2& background,
            const gridpp::Points& points,
            const vec& pobs,
            const vec& pratios,
            const vec& pbackground,
            const gridpp::StructureFunction& structure,
            int max_points);

    vec2 optimal_interpolation_transform(const gridpp::Grid& bgrid,
            const vec2& background,
            float bsigma,
            const gridpp::Points& points,
            const vec& pobs,
            const vec& psigma,
            const vec& pbackground,
            const gridpp::StructureFunction& structure,
            int max_points,
            const gridpp::Transform& transform);

    vec3 optimal_interpolation_ensi(const gridpp::Grid& bgrid,
            const vec3& input,
            const gridpp::Points& points,
            const vec& pobs,
            const vec& psigmas,
            const vec2& pbackground,
            const gridpp::StructureFunction& structure,
            int max_points);

    /****************************************
     * Neighbourhood methods                *
     ****************************************/

    vec2 neighbourhood(const vec2& input, int radius, Statistic statistic);

    vec2 neighbourhood_ens(const vec3& input, int radius, Statistic statistic);

    vec2 neighbourhood_quantile(const vec2& input, float quantile, int radius);

    vec2 neighbourhood_quantile_ens(const vec3& input, float quantile, int radius);

    vec2 neighbourhood_quantile_fast(const vec2& input, float quantile, int radius, const vec& thresholds);

    vec2 neighbourhood_quantile_ens_fast(const vec3& input, float quantile, int radius, const vec& thresholds);

    vec2 neighbourhood_brute_force(const vec2& input, int radius, Statistic statistic);

    vec get_neighbourhood_thresholds(const vec2& input, int num_thresholds);

    vec get_neighbourhood_thresholds(const vec3& input, int num_thresholds);

    /****************************************
     * Calibration methods                  *
     ****************************************/
    vec2 quantile_mapping(const vec2& input, const vec& x, const vec& y, gridpp::Extrapolation policy);

    vec quantile_mapping(const vec& input, const vec& x, const vec& y, gridpp::Extrapolation policy);

    float quantile_mapping(float input, const vec& ref, const vec& fcst, gridpp::Extrapolation policy);

    vec2 fill(const Grid& igrid, const vec2& input, const Points& points, const vec& radii, float value, bool outside);

    /****************************************
     * Downscaling methods                  *
     ****************************************/
    vec2 nearest(const Grid& igrid, const Grid& ogrid, const vec2 ivalues);
    vec nearest(const Grid& igrid, const Points& opoints, const vec2 ivalues);
    vec2 nearest(const Grid& igrid, const Points& opoints, const vec3 ivalues);

    vec2 bilinear(const Grid& igrid, const Grid& ogrid, const vec2 ivalues);

    vec bilinear(const Grid& igrid, const Points& opoints, const vec2 ivalues);

    vec2 simple_gradient(const Grid& igrid, const Grid& ogrid, const vec2 ivalues, float elev_gradient);
    vec simple_gradient(const Grid& igrid, const Points& opoints, const vec2 ivalues, float elev_gradient);

    /****************************************
     * Diagnosing methods                   *
     ****************************************/

    float wind_speed(float xwind, float ywind);
    vec wind_speed(const vec& xwind, const vec& ywind);
    vec2 wind_speed(const vec2& xwind, const vec2& ywind);

    float dewpoint(float temperature, float relative_humidity);

    vec dewpoint(const vec& temperature, const vec& relative_humidity);

    float relative_humidity(float temperature, float dewpoint);

    float wetbulb(float temperature, float pressure, float relative_humidity);

    float pressure(float ielev, float oelev, float ipressure, float itemperature=288.15);

    float qnh(float pressure, float altitude);

    vec qnh(const vec& pressure, const vec& altitude);

    void set_omp_threads(int num);

    void initialize_omp();

    namespace util {
        // vec2 calc_gradient(const vec2& values, const vec2& aux, int radius);
        // ivec regression(const vec& x, const vec& y);
        double clock();
        void debug(std::string string);
        void warning(std::string string);
        void error(std::string string);
        bool is_valid(float value);
        float calc_statistic(const vec& array, Statistic statistic);
        float calc_quantile(const vec& array, float quantile);
        vec calc_statistic(const vec2& array, Statistic statistic);
        vec calc_quantile(const vec2& array, float quantile=gridpp::MV);
        int num_missing_values(const vec2& iArray);

        int get_lower_index(float iX, const std::vector<float>& iValues);

        int get_upper_index(float iX, const std::vector<float>& iValues);

        float interpolate(float x, const std::vector<float>& iX, const std::vector<float>& iY);
        void not_implemented_error();

        vec2 init_vec2(int Y, int X, float value=gridpp::MV);

        vec calc_even_quantiles(const vec& values, int num);

        bool compatible_size(const Grid& grid, const vec2& v);
        bool compatible_size(const Grid& grid, const vec3& v);

        float* test_array(float* v, int n);
        float test_vec_input(const vec& input);
        int test_ivec_input(const ivec& input);
        float test_vec2_input(const vec2& input);
        float test_vec3_input(const vec3& input);
        vec test_vec_output();
        vec2 test_vec2_output();
        vec3 test_vec3_output();

        ivec test_ivec_output();
        ivec2 test_ivec2_output();
    }

    class Point {
        public:
            Point(float lat, float lon, float elev=gridpp::MV, float laf=gridpp::MV);
            float lat;
            float lon;
            float elev;
            float laf;
    };
    class StructureFunction {
        public:
            StructureFunction(float localization_distance);
            virtual float corr(const Point& p1, const Point& p2) const = 0;
            float localization_distance() const;
            virtual StructureFunction* clone() const = 0;
        protected:
            float barnes_rho(float dist, float length) const;
            float cressman_rho(float dist, float length) const;
            float mLocalizationDistance;
    };
    class BarnesStructure: public StructureFunction {
        public:
            BarnesStructure(float h, float v=0, float w=0, float hmax=gridpp::MV);
            float corr(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
        private:
            float mH;
            float mV;
            float mW;
    };

    class CressmanStructure: public StructureFunction {
        public:
            CressmanStructure(float h, float v=0, float w=0);
            float corr(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
        private:
            float mH;
            float mV;
            float mW;
    };
    class CrossValidation: public StructureFunction {
        public:
            CrossValidation(StructureFunction& structure, float dist=0);
            float corr(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
        private:
            StructureFunction* m_structure;
            float m_dist;
    };

    class Transform {
        public:
            virtual float forward(float value) const = 0;
            virtual float backward(float value) const = 0;
            virtual vec2 forward(const vec2& input) const;
            virtual vec2 backward(const vec2& input) const;
    };
    class Identity : public Transform {
        public:
            float forward(float value) const;
            float backward(float value) const;
    };
    class Log : public Transform {
        public:
            float forward(float value) const;
            float backward(float value) const;
    };
    class BoxCox : public Transform {
        public:
            BoxCox(float threshold);
            float forward(float value) const;
            float backward(float value) const;
        private:
            float mThreshold;
    };

    class KDTree {
        public:
            KDTree(vec lats, vec lons);
            KDTree& operator=(KDTree other);
            KDTree(const KDTree& other);
            KDTree() {};

            int get_nearest_neighbour(float lat, float lon) const;

            ivec get_neighbours(float lat, float lon, float radius) const;

            ivec get_neighbours_with_distance(float lat, float lon, float radius, vec& distances) const;

            int get_num_neighbours(float lat, float lon, float radius) const;

            ivec get_closest_neighbours(float lat, float lon, int num) const;


            static bool convert_coordinates(const vec& lats, const vec& lons, vec& x_coords, vec& y_coords, vec& z_coords);

            static bool convert_coordinates(float lat, float lon, float& x_coord, float& y_coord, float& z_coord);
            static float deg2rad(float deg);
            static float rad2deg(float deg);
            static float calc_distance(float lat1, float lon1, float lat2, float lon2);
            static float calc_distance(float x0, float y0, float z0, float x1, float y1, float z1);
            static float calc_distance_fast(float lat1, float lon1, float lat2, float lon2);
            vec get_lats() const;
            vec get_lons() const;
            int size() const;
        protected:
            typedef boost::geometry::model::point<float, 3, boost::geometry::cs::cartesian> point;
            typedef std::pair<point, unsigned> value;
            typedef boost::geometry::model::box<point> box;
            boost::geometry::index::rtree< value, boost::geometry::index::quadratic<16> > mTree;
            vec mLats;
            vec mLons;
    };

    class Points  {
        public:
            Points();
            Points(vec lats, vec lons, vec elevs=vec(), vec lafs=vec());
            Points& operator=(Points other);
            Points(const Points& other);
            int get_nearest_neighbour(float lat, float lon) const;
            ivec get_neighbours(float lat, float lon, float radius) const;
            ivec get_neighbours_with_distance(float lat, float lon, float radius, vec& distances) const;
            int get_num_neighbours(float lat, float lon, float radius) const;
            ivec get_closest_neighbours(float lat, float lon, int num) const;

            vec get_lats() const;
            vec get_lons() const;
            vec get_elevs() const;
            vec get_lafs() const;
            int size() const;
            ivec get_in_domain_indices(const Grid& grid) const;
            Points get_in_domain(const Grid& grid) const;
        private:
            KDTree mTree;
            vec mLats;
            vec mLons;
            vec mElevs;
            vec mLafs;
    };

    class Grid {
        public:
            Grid();
            Grid(vec2 lats, vec2 lons, vec2 elevs=vec2(), vec2 lafs=vec2());
            ivec get_nearest_neighbour(float lat, float lon) const;
            ivec2 get_neighbours(float lat, float lon, float radius) const;
            ivec2 get_neighbours_with_distance(float lat, float lon, float radius, vec& distances) const;
            int get_num_neighbours(float lat, float lon, float radius) const;
            ivec2 get_closest_neighbours(float lat, float lon, int num) const;

            vec2 get_lats() const;
            vec2 get_lons() const;
            vec2 get_elevs() const;
            vec2 get_lafs() const;
            ivec size() const;
        private:
            KDTree mTree;
            int mX;
            vec2 get_2d(vec input) const;
            ivec get_indices(int index) const;
            ivec2 get_indices(ivec indices) const;
            vec2 mLats;
            vec2 mLons;
            vec2 mElevs;
            vec2 mLafs;
    };
};
#endif