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>
#include <boost/math/distributions/gamma.hpp>
#include <boost/math/distributions/normal.hpp>
#ifdef _OPENMP
    #include <omp.h>
#endif
#include <exception>

#define GRIDPP_VERSION "0.7.0"
#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 these 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;
    static const float lapse_rate=0.0065;
    static const float standard_surface_temperature = 288.15;
    static const float gravit = 9.80665;
    static const float molar_mass = 0.0289644;
    static const float gas_constant_mol = 8.31447;
    static const float gas_constant_si = 287.05;
    class KDTree;
    class Points;
    class Grid;
    class Point;
    class Nearest;
    class StructureFunction;
    class Transform;

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

    enum Statistic {
        Mean      = 0,  
        Min       = 10, 
        Median    = 20, 
        Max       = 30, 
        Quantile  = 40, 
        Std       = 50, 
        Variance  = 60, 
        Sum       = 70, 
        Count     = 80, 
        RandomChoice = 90, 
        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,     
    };

    enum GradientType {
        MinMax = 0,
        LinearRegression = 10,
    };

    enum Downscaler {
        Nearest = 0,      
        Bilinear = 1,     
    };

    enum ComparisonOperator {
        Lt    = 0,         
        Leq   = 10,        
        Gt    = 20,        
        Geq   = 30,        
    };
 
    vec2 optimal_interpolation(const Grid& bgrid,
            const vec2& background,
            const Points& obs_points,
            const vec& obs,
            const vec& variance_ratios,
            const vec& background_at_points,
            const StructureFunction& structure,
            int max_points,
            bool allow_extrapolation=true);

    vec optimal_interpolation(const Points& bpoints,
            const vec& background,
            const Points& obs_points,
            const vec& obs,
            const vec& variance_ratios,
            const vec& background_at_points,
            const StructureFunction& structure,
            int max_points,
            bool allow_extrapolation=true);

    vec2 optimal_interpolation_full(const Grid& bgrid,
            const vec2& background,
            const vec2& bvariance,
            const Points& obs_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,
            bool allow_extrapolation=true);

    vec optimal_interpolation_full(const Points& bpoints,
            const vec& background,
            const vec& bvariance,
            const Points& obs_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_variance,
            bool allow_extrapolation=true);

    vec3 optimal_interpolation_ensi(const Grid& bgrid,
            const vec3& background,
            const Points& obs_points,
            const vec& obs,
            const vec& obs_standard_deviations,
            const vec2& background_at_points,
            const StructureFunction& structure,
            int max_points,
            bool allow_extrapolation=true);

    vec2 optimal_interpolation_ensi(const Points& bpoints,
            const vec2& background,
            const Points& obs_points,
            const vec& obs,
            const vec& obs_standard_deviations,
            const vec2& background_at_points,
            const StructureFunction& structure,
            int max_points,
            bool allow_extrapolation=true);

    vec2 local_distribution_correction(const Grid& bgrid,
            const vec2& background,
            const Points& obs_points,
            const vec& obs,
            const vec& background_at_points,
            const StructureFunction& structure,
            float min_quantile,
            float max_quantile,
            int min_points=0);

    vec2 local_distribution_correction(const Grid& bgrid,
            const vec2& background,
            const Points& obs_points,
            const vec2& obs,
            const vec2& background_at_points,
            const StructureFunction& structure,
            float min_quantile,
            float max_quantile,
            int min_points=0);

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

    vec2 doping_square(const Grid& igrid,
            const vec2& background,
            const Points& points,
            const vec& observations,
            const ivec& halfwidths,
            float max_elev_diff=gridpp::MV);

    vec2 doping_circle(const Grid& igrid,
            const vec2& background,
            const Points& points,
            const vec& observations,
            const vec& radii,
            float max_elev_diff=gridpp::MV);
 
    vec gamma_inv(const vec& levels, const vec& shape, const vec& scale);
 
    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 calc_gradient(const vec2& base, const vec2& values, GradientType gradient_type, int halfwidth, int min_num=2, float min_range=gridpp::MV, float default_gradient=0);

    vec2 neighbourhood_search(const vec2& array, const vec2& search_array, int halfwidth, float search_target_min, float search_target_max, float search_delta, const ivec2& apply_array=ivec2());

    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());

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

    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& curve_ref, const vec& curve_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 downscaling(const Grid& igrid, const Grid& ogrid, const vec2& ivalues, Downscaler downscaler);

    vec3 downscaling(const Grid& igrid, const Grid& ogrid, const vec3& ivalues, Downscaler downscaler);

    vec downscaling(const Grid& igrid, const Points& opoints, const vec2& ivalues, Downscaler downscaler);

    vec2 downscaling(const Grid& igrid, const Points& opoints, const vec3& ivalues, Downscaler downscaler);

    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 downscale_probability(const Grid& igrid, const Grid& ogrid, const vec3& ivalues, const vec2& threshold, const ComparisonOperator& comparison_operator);

    vec2 mask_threshold_downscale_consensus(const Grid& igrid, const Grid& ogrid, const vec3& ivalues_true, const vec3& ivalues_false, const vec3& theshold_values, const vec2& threshold, const ComparisonOperator& comparison_operator, const Statistic& statistic);

    vec2 mask_threshold_downscale_quantile(const Grid& igrid, const Grid& ogrid, const vec3& ivalues_true, const vec3& ivalues_false, const vec3& theshold_values, const vec2& threshold, const ComparisonOperator& comparison_operator, const float quantile_level);

    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, Downscaler downscaler=Nearest);

    vec3 simple_gradient(const Grid& igrid, const Grid& ogrid, const vec3& ivalues, float elev_gradient, Downscaler downscaler=Nearest);

    vec simple_gradient(const Grid& igrid, const Points& opoints, const vec2& ivalues, float elev_gradient, Downscaler downscaler=Nearest);

    vec2 simple_gradient(const Grid& igrid, const Points& opoints, const vec3& ivalues, float elev_gradient, Downscaler downscaler=Nearest);

    vec2 full_gradient(const Grid& igrid, const Grid& ogrid, const vec2& ivalues,  const vec2& elev_gradient, const vec2& laf_gradient=vec2(), Downscaler downscaler=Nearest);

    vec3 full_gradient(const Grid& igrid, const Grid& ogrid, const vec3& ivalues, const vec3& elev_gradient, const vec3& laf_gradient, Downscaler downscaler=Nearest);

    vec full_gradient(const Grid& igrid, const Points& opoints, const vec2& ivalues, const vec2& elev_gradient, const vec2& laf_gradient, Downscaler downscaler=Nearest);

    vec2 full_gradient(const Grid& igrid, const Points& opoints, const vec3& ivalues, const vec3& elev_gradient, const vec3& laf_gradient, Downscaler downscaler=Nearest);

    /* Elevation and land area fraction downscaling with debug output fields. For debugging only,
    */
    vec3 full_gradient_debug(const Grid& igrid, const Grid& ogrid, const vec2& ivalues,  const vec2& elev_gradient, const vec2& laf_gradient=vec2(), Downscaler downscaler=Nearest);

    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 Grid& igrid, const Grid& ogrid, float radius);

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

    vec count(const Points& ipoints, const Points& opoints, 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);

    vec distance(const Points& ipoints, const Points& opoint, 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);

    vec gridding(const Points& opoints, const Points& ipoints, const vec& values, float radius, int min_num, Statistic statistic);

    vec2 gridding_nearest(const Grid& grid, const Points& points, const vec& values, int min_num, gridpp::Statistic statistic);

    vec gridding_nearest(const Points& opoints, const Points& ipoints, const vec& values, int min_num, gridpp::Statistic statistic);

    vec2 neighbourhood_score(const Grid& grid, const Points& points, const vec2& fcst, const vec& ref, int half_width, gridpp::Metric metric, float threshold);
 
    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 sea_level_pressure(float ps, float altitude, float temperature, float rh=gridpp::MV, float dewpoint=gridpp::MV);

    vec sea_level_pressure(const vec& ps, const vec& altitude, const vec& temperature, const vec& rh, const vec& dewpoint);

    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);

    int get_omp_threads();

    void initialize_omp();
 
    // ivec regression(const vec& x, const vec& y);

    void set_debug_level(int level);

    static int _debug_level = 0;

    int get_debug_level();

    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_level);

    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_levels);

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

    bool convert_coordinates(float lat, float lon, CoordinateType type, float& x_coord, float& y_coord, float& z_coord);

    bool is_valid_lat(float lat, CoordinateType type);

    bool is_valid_lon(float lon, CoordinateType type);

    int num_missing_values(const vec2& iArray);

    int get_lower_index(float iX, const vec& iValues);

    int get_upper_index(float iX, const vec& iValues);

    float interpolate(float x, const vec& iX, const vec& iY);

    vec interpolate(const vec& x, const vec& iX, const vec& 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 window(const vec2& array, int length, gridpp::Statistic statistic, bool before=false, bool keep_missing=false, bool missing_edges=true);

    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 vec2& b);

    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);

            Point(float lat, float lon, float elev, float laf, CoordinateType type, float x, float y, float z);
            float lat;
            float lon;
            float elev;
            float laf;
            CoordinateType type;
            float x;
            float y;
            float z;
    };

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

            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;

            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_straight_distance(float x0, float y0, float z0, float x1, float y1, float z1);

            static float calc_distance(const Point& p1, const Point& p2);

            static float calc_straight_distance(const Point& p1, const Point& p2);

            static float calc_distance_fast(float lat1, float lon1, float lat2, float lon2, CoordinateType type=Geodetic);

            vec get_lats() const;
            vec get_lons() const;
            int size() const;
            CoordinateType get_coordinate_type() const;
            const vec& get_x() const;
            const vec& get_y() const;
            const vec& get_z() 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;
            vec mX;
            vec mY;
            vec mZ;
            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);

            /* Get the index of the nearest neighbour
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param include_match Return matching points. If false, don't use this point.
             * @returns Index of nearest neighbour (-1 if there are no neighbours).
            */
            int get_nearest_neighbour(float lat, float lon, bool include_match=true) const;

            /* Get the indices of all neighbours within a search radius
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param radius Neighbourhood radius [m]
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 1D vector of indices of neighbours within radius
            */
            ivec get_neighbours(float lat, float lon, float radius, bool include_match=true) const;

            /* Get the indices and distances to all neighbours within a search radius
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param radius Neighbourhood radius [m]
             * @param distances 1D output vector of distances to each point [m]
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 1D vector of indices of neighbours within radius
            */
            ivec get_neighbours_with_distance(float lat, float lon, float radius, vec& distances, bool include_match=true) const;

            /* Get the number of neighbouring points within a search radius
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param radius Neighbourhood radius [m]
             * @param include_match Return matching points. If false, don't use this point.
             * @returns the number of neighbours within radius
            */
            int get_num_neighbours(float lat, float lon, float radius, bool include_match=true) const;

            /* Get the N nearest neighbouring points
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param num Number of nearest neighbours
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 1D vector of indices of neighbours
            */
            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);

            /* Get the index of the nearest neighbour
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 1D vector of (x, y) index (empty vector if no neighbours)
            */
            ivec get_nearest_neighbour(float lat, float lon, bool include_match=true) const;

            /* Get the indices of all neighbours within a search radius
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param radius Neighbourhood radius [m]
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 2D vector (Point, 2) of x,y-indices of neighbours within radius
            */
            ivec2 get_neighbours(float lat, float lon, float radius, bool include_match=true) const;

            /* Get the indices and distances to all neighbours within a search radius
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param radius Neighbourhood radius [m]
             * @param distances 1D output vector of distances to each point [m]
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 2D vector (Point, 2) of x,y-indices of neighbours within radius
            */
            ivec2 get_neighbours_with_distance(float lat, float lon, float radius, vec& distances, bool include_match=true) const;

            /* Get the number of neighbouring points within a search radius
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param radius Neighbourhood radius [m]
             * @param include_match Return matching points. If false, don't use this point.
             * @returns the number of neighbours within radius
            */
            int get_num_neighbours(float lat, float lon, float radius, bool include_match=true) const;

            /* Get the N nearest neighbouring points
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param num Number of nearest neighbours
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 2D vector (Point, 2) of x,y-indices of neighbours
            */
            ivec2 get_closest_neighbours(float lat, float lon, int num, bool include_match=true) const;

            /* Get the N nearest neighbouring points
             * @param lat Lookup latitude (or y-axis value)
             * @param lon Lookup longitude (or x-axis value)
             * @param Y1_out Number of nearest neighbours
             * @param include_match Return matching points. If false, don't use this point.
             * @returns 2D vector (Point, 2) of x,y-indices of neighbours
            */
            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") { };
    };

    class StructureFunction {
        public:
            StructureFunction(float localization_distance=0);
            virtual float corr(const Point& p1, const Point& p2) const = 0;

            virtual float corr_background(const Point& p1, const Point& p2) const;

            virtual float localization_distance(const Point& p) const;
            virtual StructureFunction* clone() const = 0;
            static const float default_min_rho;
        protected:
            float barnes_rho(float dist, float length) const;

            float cressman_rho(float dist, float length) const;
            float m_localization_distance;
    };
    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;
            float localization_distance(const Point& p) 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);

            BarnesStructure(Grid grid, vec2 h, vec2 v, vec2 w, float min_rho=StructureFunction::default_min_rho);
            float corr(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
            float localization_distance(const Point& p) const;
        private:
            Grid m_grid;
            vec2 mH;
            vec2 mV;
            vec2 mW;
            float m_min_rho;
            bool m_is_spatial;
    };

    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=MV);
            float corr(const Point& p1, const Point& p2) const;
            float corr_background(const Point& p1, const Point& p2) const;
            StructureFunction* clone() const;
            float localization_distance(const Point& p) 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 Gamma : public Transform {
        public:
            Gamma(float shape, float scale, float tolerance=0.01);
            using Transform::forward;
            using Transform::backward;
            float forward(float value) const;
            float backward(float value) const;
        private:
            float m_tolerance;
            boost::math::gamma_distribution<> m_gamma_dist;
            boost::math::normal m_norm_dist;
    };
};
#endif