titanlib.h

Go to the documentation of this file.

#ifndef TITANLIB_H
#define TITANLIB_H
#include <iostream>
#include <vector>
#include <assert.h>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#ifdef _OPENMP
    #include <omp.h>
#endif

#define TITANLIB_VERSION "0.3.3"
#define __version__ TITANLIB_VERSION

namespace titanlib { 
    // Preferred vector types
    typedef std::vector<int> ivec;
    typedef std::vector<float> vec;
    typedef std::vector<double> dvec;
    typedef std::vector<vec> vec2; 
    static const float MV = NAN;
    static const float MV_CML = -999;
    static const float pi = 3.14159265;
    static const double radius_earth = 6.378137e6;
    enum CoordinateType {
        Geodetic = 0,      
        Cartesian = 1,     
    };

    enum BackgroundType {
        VerticalProfile = 0,
        VerticalProfileTheilSen = 1,
        MeanOuterCircle = 2,
        MedianOuterCircle = 3,
        External = 4,
    };
 
    enum ConditionType {
        Eq = 0,
        Gt = 1,
        Geq = 2,
        Lt = 3,
        Leq = 4,
    };

    class Points;
 
    ivec sct(const Points& points,
            const vec& values,
            int num_min,
            int num_max,
            float inner_radius,
            float outer_radius,
            int num_iterations,
            int num_min_prof,
            float min_elev_diff,
            float min_horizontal_scale,
            float vertical_scale,
            const vec& pos,
            const vec& neg,
            const vec& eps2,
            vec& prob_gross_error,
            vec& rep,
            const ivec& obs_to_check=ivec());

    ivec sct_resistant( const Points& points,
                        const vec& values,
                        const ivec& obs_to_check,
                        const vec& background_values,
                        BackgroundType background_elab_type,
                        int num_min_outer,
                        int num_max_outer,
                        float inner_radius,
                        float outer_radius,
                        int num_iterations,
                        int num_min_prof,
                        float min_elev_diff,
                        float min_horizontal_scale,
                        float max_horizontal_scale,
                        int kth_closest_obs_horizontal_scale,
                        float vertical_scale,
                        const vec& value_mina,
                        const vec& value_maxa,
                        const vec& value_minv,
                        const vec& value_maxv,
                        const vec& eps2,
                        const vec& tpos,
                        const vec& tneg,
                        bool debug,
                        bool basic,
                        vec& scores);

    ivec sct_dual( const Points& points,
                   const vec& values,
                   const ivec& obs_to_check,
                   const vec& event_thresholds,
                   ConditionType condition,
                   int num_min_outer,
                   int num_max_outer,
                   float inner_radius,
                   float outer_radius,
                   int num_iterations,
                   float min_horizontal_scale,
                   float max_horizontal_scale,
                   int kth_closest_obs_horizontal_scale,
                   float vertical_scale,
                   const vec& test_thresholds,
                   bool debug);

    ivec fgt( const Points& points,
              const vec& values,
              const ivec& obs_to_check,
              const vec& background_values,
              const vec& background_uncertainties,
              BackgroundType background_elab_type,
              int num_min_outer,
              int num_max_outer,
              float inner_radius,
              float outer_radius,
              int num_iterations,
              int num_min_prof,
              float min_elev_diff,
              const vec& value_mina,
              const vec& value_maxa,
              const vec& value_minv,
              const vec& value_maxv,
              const vec& tpos,
              const vec& tneg,
              bool debug,
              bool basic,
              vec& scores);

    ivec range_check(const vec& values,
            const vec& min,
            const vec& max);

    ivec range_check_climatology(const Points& points,
            const vec& values,
            int unixtime,
            const vec& pos,
            const vec& neg);

    ivec buddy_check(const Points& points,
            const vec& values,
            const vec& radius,
            const ivec& num_min,
            float threshold,
            float max_elev_diff,
            float elev_gradient,
            float min_std,
            int num_iterations,
            const ivec& obs_to_check = ivec());

    ivec buddy_event_check(const Points& points,
            const vec& values,
            const vec& radius,
            const ivec& num_min,
            float event_threshold,
            float threshold,
            float max_elev_diff,
            float elev_gradient,
            int num_iterations,
            const ivec& obs_to_check = ivec());

    ivec isolation_check(const Points& points,
            int num_min,
            float radius,
            float vertical_radius=MV);

    ivec isolation_check(const Points& points,
            const ivec& num_min,
            const vec& radius,
            const vec& vertical_radius=vec());

    ivec duplicate_check(const Points& points, float radius, float vertical_range=titanlib::MV);

    ivec metadata_check(const Points& points, bool check_lat=true, bool check_lon=true, bool check_elev=true, bool check_laf=true);
 
    vec lag_reduction_filter(const vec& times, const vec& values, float a=0.5, float b=0.5, float k1=0.25, float k2=0.25, int n=10);
 
    void set_omp_threads(int num);

    int get_omp_threads();

    void initialize_omp();
 
    std::string version();

    double clock();
    bool is_valid(float value);

    vec compute_vertical_profile(const vec& elevs, const vec& oelevs, const vec& values, int num_min_prof, double min_elev_diff, bool debug=false);

    vec compute_vertical_profile_Theil_Sen(const vec& elevs, const vec& oelevs, const vec& values, int num_min_prof, double min_elev_diff, bool debug);

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

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

    vec interpolate_to_points(const vec2& input_lats, const vec2& input_lons, const vec2& input_values, const vec& output_lats, const vec& output_lons);
    float deg2rad(float deg);
    float calc_distance(float lat1, float lon1, float lat2, float lon2);
    float calc_distance(float x0, float y0, float z0, float x1, float y1, float z1);

    float compute_quantile(double quantile, const vec& array);
    float find_k_closest(const vec& array, int k);
    template <class T> T subset(const T& array, const ivec& indices) {
        if(array.size() == 0)
            return array;
        if(indices.size() == 0) {
            return array;
        }

        T new_array(indices.size());
        for(int i = 0; i < indices.size(); i++) {
            new_array[i] = array[indices[i]];
        }
        return new_array;
    }
    Points subset(const Points& input, const ivec& indices);
    template <class T> void unsubset(const T& array, T& orig_array, const ivec& indices) {
        orig_array.clear();
        orig_array.resize(indices.size());
        assert(array.size() == indices.size());
        for(int i = 0; i < indices.size(); i++) {
            orig_array[indices[i]] = array[i];
        }
    }


    vec background(const vec& elevs, const vec& values, int num_min_prof, float min_elev_diff, float value_minp, float value_maxp, BackgroundType background_type, const vec& external_background_values, const ivec& indices_global_outer, bool debug);

    bool invert_matrix (const boost::numeric::ublas::matrix<float>& input, boost::numeric::ublas::matrix<float>& inverse);

    bool set_indices( const ivec& indices_global_outer_guess, const ivec& obs_test, const ivec& dqcflags, const vec& dist_outer_guess, float inner_radius, int test_just_this, ivec& indices_global_outer, ivec& indices_global_test, ivec& indices_outer_inner, ivec& indices_outer_test, ivec& indices_inner_test);
    template<class T1, class T2> struct sort_pair_first {
        bool operator()(const std::pair<T1,T2>&left, const std::pair<T1,T2>&right) {
            return left.first < right.first;
        };
    }; 
    float* test_array(float* v, int n);

    void test_not_implemented_exception();
    // ivec nearest_neighbours(const vec& lats, const vec& lons, float radius, float lat, float lon);

    // bool prioritize(const vec& values, const ivec& priority, float distance, ivec& flags);

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


            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 operator()(value const& v) const;
                private:
                    float radius;
                    point p;
            };
            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;
            CoordinateType get_coordinate_type() const;
        private:
            KDTree mTree;
            vec mLats;
            vec mLons;
            vec mElevs;
            vec mLafs;
    };

    class Dataset {
        public:
            Dataset(Points points, vec ivalues);
            void range_check(const vec& min, const vec& max, const ivec& indices=ivec(1, -1));
            void range_check_climatology(int unixtime, const vec& pos, const vec& neg, const ivec& indices=ivec(1, -1));
            void sct(int num_min, int num_max, float inner_radius, float outer_radius, int num_iterations, int num_min_prof, float min_elev_diff, float min_horizontal_scale, float vertical_scale, const vec& t2pos, const vec& t2neg, const vec& eps2, vec& prob_gross_error, vec& rep, const ivec& indices=ivec(1, -1));
            void buddy_check(const vec& radius, const ivec& num_min, float threshold, float max_elev_diff, float elev_gradient, float min_std, int num_iterations, const ivec& obs_to_check=ivec(), const ivec& indices=ivec(1, -1));
            void buddy_event_check(const vec& radius, const ivec& num_min, float event_threshold, float threshold, float max_elev_diff, float elev_gradient, int num_iterations, const ivec& obs_to_check=ivec(), const ivec& indices=ivec(1, -1));
            void isolation_check(int num_min, float radius, float vertical_radius=MV, const ivec& indices=ivec(1, -1));
            void isolation_check(const ivec& num_min, const vec& radius, const vec& vertical_radius=vec(), const ivec& indices=ivec(1, -1));
            void duplicate_check(float radius, float vertical_range=titanlib::MV, const ivec& indices=ivec(1, -1));
            void dem_check(const vec& dem, float max_elev_diff);
            void external_check(const ivec& flags);
            void metadata_check(bool check_lat=true, bool check_lon=true, bool check_elev=true, bool check_laf=true, const ivec& indices=ivec(1, -1));

            Points get_points() const;
            vec get_values() const;
            ivec get_flags() const;
            void set_values(vec ivalues);
            void set_flags(ivec ivalues);
            void set_points(Points ipoints);

            Points points;
            vec values;
            ivec flags;
        private:
            /* Get a subset of array where flags=0 and for indices. This function is different than
             * titanlib::subset in that it takes into account flags.
             * @param array Vector of values
             * @param indices Vector of indices
             * @return Vector of values
            */
            template <class T> T subset_valid(const T& array, const ivec& indices=ivec(1, -1)) {
                if(array.size() != 1 && array.size() != flags.size()) {
                    std::stringstream ss;
                    ss << "Array (" << array.size() << ") must be either size 1 or the same as dataset size (" << flags.size() << ")";
                    throw std::invalid_argument(ss.str());
                }
                if(indices.size() == 0)
                    return T();

                ivec indices0 = indices;
                if(indices.size() == 1 && indices[0] == -1) {
                    // Use all indices
                    indices0.clear();
                    indices0.resize(flags.size());
                    for(int i = 0; i < flags.size(); i++)
                        indices0[i] = i;
                }
                T new_array;
                new_array.reserve(indices0.size());
                for(int i = 0; i < indices0.size(); i++) {
                    int index = indices0[i];
                    if(flags[index] == 0) {
                        if(array.size() == 1)
                            // Array is broadcast
                            new_array.push_back(array[0]);
                        else {
                            assert(index < array.size());
                            new_array.push_back(array[index]);
                        }
                    }
                }
                return new_array;
            }
            /*  Same as T subset_valid, except for Points */
            Points subset_valid(const Points& input, const ivec& indices=ivec(1, -1));

            ivec get_unflagged_indices();

            // Create a vector with results, but only for locations that have not been flagged
            // To broadcase the values to point, use a singleton array
            // An empty array will return an empty array
            template <class T> T get_unflagged(const T& array) {
                if(array.size() == 0)
                    return array;

                if(array.size() > 1 && array.size() < flags.size()) {
                    throw std::invalid_argument("Array is shorter than flags");
                }

                ivec indices = get_unflagged_indices();
                T output(indices.size());
                for(int i = 0; i < indices.size(); i++) {
                    if(array.size() > 1) {
                        int index = indices[i];
                        assert(index < array.size());
                        output[i] = array[index];
                    }
                    else
                        output[i] = array[0];
                }
                return output;
            }
            Points get_unflagged_points();

            /* Same as T subset_valid, except the subsetted indices are returned */
            ivec subset_valid(const ivec& indices);

            /* Merge the flags into the objects flags. If existing flags are 1 but new_flags are
             * 0, the existing flags are not updated.
             * @param new_flags Vector of new flags
             * @param indices The location indices that flags are valid for
            */
            void merge_simple(const ivec& new_flags, ivec indices);

            /* Merge new flags into the object's flags. If existing flags are 1 but new_flags are
             * 0, the existing flags are not updated.
             * @param new_flags Vector of new flags
             * @param indices The location indices that flags are valid for
            */
            void merge(const ivec& new_flags, ivec indices);
    };
    class not_implemented_exception: public std::logic_error
    {
        public:
            not_implemented_exception() : std::logic_error("Function not yet implemented") { };
    };

}
#endif