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
#include <exception>

#define GRIDPP_VERSION "0.6.0.dev1"
#define __version__ GRIDPP_VERSION

namespace gridpp { 
    // 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;
    typedef std::vector<ivec2> ivec3;

    // Only use when double is required
    typedef std::vector<double> dvec;
    typedef std::vector<dvec> dvec2; 
    static const float MV = NAN;
    static const float MV_CML = -999;
    static const float pi = 3.14159265;
    static const double radius_earth = 6.378137e6;
    class KDTree;
    class Points;
    class Grid;
    class Point;
    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  
    };

    enum Metric {
        Ets       = 0,  
        Ts        = 1,  
        Kss       = 20, 
        Pc        = 30, 
        Bias      = 40, 
        Hss       = 50, 
    };

    enum CorrectionType {
        Qq        = 0,        
        Multiplicative = 10,  
        Additive  = 20,       
    };

    enum CoordinateType {
        Geodetic = 0,      
        Cartesian = 1,     
    };
 
    vec2 optimal_interpolation(const Grid& bgrid,
            const vec2& background,
            const Points& points,
            const vec& pobs,
            const vec& pratios,
            const vec& pbackground,
            const StructureFunction& structure,
            int max_points);

    vec optimal_interpolation(const Points& bpoints,
            const vec& background,
            const Points& points,
            const vec& pobs,
            const vec& pratios,
            const vec& pbackground,
            const StructureFunction& structure,
            int max_points);

    vec2 optimal_interpolation_full(const Grid& bgrid,
            const vec2& background,
            const vec2& bvariance,
            const Points& points,
            const vec& obs,
            const vec& obs_variance,
            const vec& background_at_points,
            const vec& bvariance_at_points,
            const StructureFunction& structure,
            int max_points,
            vec2& analysis_variance);

    vec optimal_interpolation_full(const Points& bpoints,
            const vec& background,
            const vec& bvariance,
            const Points& points,
            const vec& obs,
            const vec& obs_variance,
            const vec& background_at_points,
            const vec& bvariance_at_points,
            const StructureFunction& structure,
            int max_points,
            vec& analysis_sigmas);

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

    vec2 optimal_interpolation_ensi(const Points& bpoints,
            const vec2& background,
            const Points& points,
            const vec& pobs,
            const vec& psigmas,
            const vec2& pbackground,
            const StructureFunction& structure,
            int max_points);

    vec2 fill(const Grid& igrid, const vec2& input, const Points& points, const vec& radii, float value, bool outside);
    vec2 doping(const Grid& igrid, const vec2& input, const Points& points, const vec& values, int half_width, float max_elev_diff=gridpp::MV);
 
    vec2 neighbourhood(const vec2& input, int halfwidth, Statistic statistic);

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

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

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

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

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

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

    vec2 neighbourhood_quantile_fast(const vec3& input, const vec2& quantile, int halfwidth, const vec& thresholds);

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

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

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

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

    vec2 neighbourhood_ens(const vec3& input, int halfwidth, Statistic statistic);
    vec2 neighbourhood_quantile_ens(const vec3& input, float quantile, int halfwidth);
    vec2 neighbourhood_quantile_ens_fast(const vec3& input, float quantile, int radius, const vec& thresholds); 
    vec quantile_mapping_curve(const vec& ref, const vec& fcst, vec& output_fcst, vec quantiles=vec());

    vec2 metric_optimizer_curve(const vec& ref, const vec& fcst, const vec& thresholds, Metric metric);

    float apply_curve(float fcst, const vec& curve_ref, const vec& curve_fcst, Extrapolation policy_below, Extrapolation policy_above);

    vec apply_curve(const vec& fcst, const vec& curve_ref, const vec& curve_fcst, Extrapolation policy_below, Extrapolation policy_above);

    vec2 apply_curve(const vec2& fcst, const vec& curve_ref, const vec& curve_fcst, Extrapolation policy_below, Extrapolation policy_above);

    vec2 apply_curve(const vec2& fcst, const vec3& curve_ref, const vec3& curve_fcst, Extrapolation policy_below, Extrapolation policy_above);

    vec monotonize_curve(vec curve_ref, vec curve_fcst, vec& output_fcst);

    float get_optimal_threshold(const vec& ref, const vec& fcst, float threshold, Metric metric);

    float calc_score(float a, float b, float c, float d, Metric metric);
    float calc_score(const vec& ref, const vec& fcst, float threshold, Metric metric);
    float calc_score(const vec& ref, const vec& fcst, float threshold, float fthreshold, Metric metric);

    vec2 correction(const Grid& rgrid, const vec2& rvalues, const Points& npoints, const vec& nvalues, float mean_radius, float outer_radius, float inner_radius, int min_num, int max_num, CorrectionType type, ivec2& count);
    // Apply correction based on multiple timesteps
    vec2 correction(const Grid& rgrid, const vec3& rvalues, const Points& npoints, const vec2& nvalues, const vec2& apply_values, float mean_radius, float outer_radius, float inner_radius, int min_num, int max_num, CorrectionType type, ivec2& count);
 
    vec2 nearest(const Grid& igrid, const Grid& ogrid, const vec2& ivalues);
    vec3 nearest(const Grid& igrid, const Grid& ogrid, const vec3& ivalues);

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

    vec nearest(const Points& ipoints, const Points& opoints, const vec& ivalues);
    vec2 nearest(const Points& ipoints, const Points& opoints, const vec2& ivalues);

    vec2 nearest(const Points& ipoints, const Grid& ogrid, const vec& ivalues);
    vec3 nearest(const Points& ipoints, const Grid& ogrid, const vec2& ivalues);

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

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

    vec2 simple_gradient(const Grid& igrid, const Grid& ogrid, const vec2& ivalues, float elev_gradient);
    vec3 simple_gradient(const Grid& igrid, const Grid& ogrid, const vec3& ivalues, float elev_gradient);

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

    vec2 smart(const Grid& igrid, const Grid& ogrid, const vec2& ivalues, int num, const StructureFunction& structure); 
    vec count(const Grid& grid, const Points& points, float radius);

    vec2 count(const Points& points, const Grid& grid, float radius);

    vec distance(const Grid& grid, const Points& points, int num=1);

    vec2 distance(const Grid& igrid, const Grid& ogrid, int num=1);

    vec2 distance(const Points& points, const Grid& grid, int num=1);

    vec2 fill_missing(const vec2& values);

    vec2 gridding(const Grid& grid, const Points& points, const vec& values, float radius, int min_num, Statistic statistic);
 
    float dewpoint(float temperature, float relative_humidity);

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

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

    vec pressure(const vec& ielev, const vec& oelev, const vec& ipressure, const vec& itemperature);

    float qnh(float pressure, float altitude);

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

    float relative_humidity(float temperature, float dewpoint);

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

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

    vec wetbulb(const vec& temperature, const vec& pressure, const vec& relative_humidity);

    float wind_speed(float xwind, float ywind);

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

    float wind_direction(float xwind, float ywind);

    vec wind_direction(const vec& xwind, const vec& ywind);
 
    void set_omp_threads(int num);

    void initialize_omp();
 
    // vec2 calc_gradient(const vec2& values, const vec2& aux, int radius);
    // ivec regression(const vec& x, const vec& y);

    Statistic get_statistic(std::string name);

    std::string version();

    double clock();
    void debug(std::string string);
    void warning(std::string string);
    void error(std::string string);
    void future_deprecation_warning(std::string function, std::string other="");
    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=MV);

    vec2 calc_quantile(const vec3& array, const vec2& quantile);
    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);

    ivec2 init_ivec2(int Y, int X, int value);
    vec2 init_vec2(int Y, int X, float value=MV);

    ivec3 init_ivec3(int Y, int X, int E, int value);
    vec3 init_vec3(int Y, int X, int E, float value=MV);

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

    vec2 calc_gradient(const Grid& grid, const vec2& base, const vec2& values, int radius, int min_num=2, float min_range=0, float default_gradient=0);

    bool compatible_size(const Grid& grid, const vec2& v);
    bool compatible_size(const Grid& grid, const vec3& v);
    bool compatible_size(const Points& points, const vec& v);
    bool compatible_size(const Points& points, const vec2& v);
    bool compatible_size(const vec2& a, const vec3& b);
    bool compatible_size(const vec3& a, const vec3& b);

    bool point_in_rectangle(const Point& A, const Point& B, const Point& C, const Point& D, const Point& m );
 
    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();
    ivec test_ivec_output();
    vec2 test_vec2_output();
    ivec2 test_ivec2_output();
    vec3 test_vec3_output();
    ivec3 test_ivec3_output();

    float test_vec_argout(vec& distances);
    float test_vec2_argout(vec2& distances);

    void test_not_implemented_exception();

    static const float swig_default_value = -1;

    class Point {
        public:
            Point(float lat, float lon, float elev=MV, float laf=MV, CoordinateType type=Geodetic);
            float lat;
            float lon;
            float elev;
            float laf;
            CoordinateType type;
    };
    class StructureFunction {
        public:
            StructureFunction(float localization_distance);
            virtual float corr(const Point& p1, const Point& p2) const = 0;
            virtual float corr_background(const Point& p1, const Point& p2) const;
            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 MultipleStructure: public StructureFunction {
        public:
            MultipleStructure(const StructureFunction& structure_h, const StructureFunction& structure_v, const StructureFunction& structure_w);
            float corr(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
        private:
            StructureFunction* m_structure_h;
            StructureFunction* m_structure_v;
            StructureFunction* m_structure_w;
    };
    class BarnesStructure: public StructureFunction {
        public:
            BarnesStructure(float h, float v=0, float w=0, float hmax=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 VariedStructure: public StructureFunction {
        public:
            VariedStructure(const Points& points, float h_max);
            float corr(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
        private:
            const Points& m_points;
            float m_h_max;
            // StructureFunction* m_structure;
    };
    class CrossValidation: public StructureFunction {
        public:
            CrossValidation(StructureFunction& structure, float dist=MV);
            float corr(const Point& p1, const Point& p2) const;
            float corr_background(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
        private:
            StructureFunction* m_structure;
            float m_dist;
    };

    class Transform {
        public:
            // Note these cannot be pure virtual, otherwise SWIG does not expose
            // the vector functions (with the same name) in python. Therefore, make sure these
            // functions are overloaded in the subclass implementation
            virtual float forward(float value) const;
            virtual float backward(float value) const;

            vec forward(const vec& input) const;
            vec backward(const vec& input) const;
            vec2 forward(const vec2& input) const;
            vec2 backward(const vec2& input) const;
            vec3 forward(const vec3& input) const;
            vec3 backward(const vec3& input) const;
    };
    class Identity : public Transform {
        public:
            // SWIG requires these "using" statements to enable the vectorized versions in the
            // subclasses
            using Transform::forward;
            using Transform::backward;
            float forward(float value) const;
            float backward(float value) const;
    };
    class Log : public Transform {
        public:
            using Transform::forward;
            using Transform::backward;
            float forward(float value) const;
            float backward(float value) const;
    };
    class BoxCox : public Transform {
        public:
            BoxCox(float threshold);
            using Transform::forward;
            using Transform::backward;
            float forward(float value) const;
            float backward(float value) const;
        private:
            float mThreshold;
    };

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

            int get_nearest_neighbour(float lat, float lon, bool include_match=true) const;

            ivec get_neighbours(float lat, float lon, float radius, bool include_match=true) const;

            ivec get_neighbours_with_distance(float lat, float lon, float radius, vec& distances, bool include_match=true) const;

            int get_num_neighbours(float lat, float lon, float radius, bool include_match=true) const;

            ivec get_closest_neighbours(float lat, float lon, int num, bool include_match=true) const;


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

            bool convert_coordinates(float lat, float lon, float& x_coord, float& y_coord, float& z_coord) const;
            static float deg2rad(float deg);
            static float rad2deg(float deg);
            static float calc_distance(float lat1, float lon1, float lat2, float lon2, CoordinateType type=Geodetic);
            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, CoordinateType type=Geodetic);
            static float calc_distance_fast(const Point& p1, const Point& p2);
            vec get_lats() const;
            vec get_lons() const;
            int size() const;
            CoordinateType get_coordinate_type() 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;
            CoordinateType mType;

            struct within_radius {
                public:
                    within_radius(point p, float radius, bool include_match);
                    bool operator()(value const& v) const;
                private:
                    float radius;
                    point p;
                    bool include_match;
            };
            struct is_not_equal {
                public:
                    is_not_equal(point p);
                    bool operator()(value const& v) const;
                private:
                    point p;
            };
    };

    class Points  {
        public:
            Points();
            Points(vec lats, vec lons, vec elevs=vec(), vec lafs=vec(), CoordinateType type=Geodetic);
            Points(KDTree tree, vec elevs=vec(), vec lafs=vec());
            Points& operator=(Points other);
            Points(const Points& other);
            // Returns -1 if there are no neighbours
            int get_nearest_neighbour(float lat, float lon, bool include_match=true) const;
            ivec get_neighbours(float lat, float lon, float radius, bool include_match=true) const;
            ivec get_neighbours_with_distance(float lat, float lon, float radius, vec& distances, bool include_match=true) const;
            int get_num_neighbours(float lat, float lon, float radius, bool include_match=true) const;
            ivec get_closest_neighbours(float lat, float lon, int num, bool include_match=true) 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;
            CoordinateType get_coordinate_type() const;
            Point get_point(int index) const;
            Points subset(const ivec& indices) 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(), CoordinateType type=Geodetic);
            ivec get_nearest_neighbour(float lat, float lon, bool include_match=true) const;
            ivec2 get_neighbours(float lat, float lon, float radius, bool include_match=true) const;
            ivec2 get_neighbours_with_distance(float lat, float lon, float radius, vec& distances, bool include_match=true) const;
            int get_num_neighbours(float lat, float lon, float radius, bool include_match=true) const;
            ivec2 get_closest_neighbours(float lat, float lon, int num, bool include_match=true) const;

            bool get_box(float lat, float lon, int& Y1_out, int& X1_out, int& Y2_out, int& X2_out) const;

            Points to_points() const;

            vec2 get_lats() const;
            vec2 get_lons() const;
            vec2 get_elevs() const;
            vec2 get_lafs() const;
            ivec size() const;
            CoordinateType get_coordinate_type() const;
            Point get_point(int y_index, int x_index) 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;
    };
    class not_implemented_exception: public std::logic_error
    {
        public:
            not_implemented_exception() : std::logic_error("Function not yet implemented") { };
    };
};
#endif