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.0-dev3"
#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,
    };

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

    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,
                        vec& scores);

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

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

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

    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);
    vec subset(const vec& array, const ivec& indices);
    Points subset(const Points& input, const ivec& indices);


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

            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());
            void range_check_climatology(int unixtime, const vec& pos, const vec& neg, const ivec& indices=ivec());
            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());
            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, const ivec& indices=ivec());
            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());
            void isolation_check(int num_min, float radius, float vertical_radius, const ivec& indices=ivec());
            void duplicate_check(float radius, float vertical_range=titanlib::MV, const ivec& indices=ivec());
            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());

            vec lats;
            vec lons;
            vec elevs;
            Points points;
            vec values;
            ivec flags;
        private:
            template <class T> T subset(const T& array, const ivec& indices) {
                if(array.size() == 1) {
                    T new_array = array;
                    return new_array;
                }
                else {
                    T new_array(indices.size());
                    for(int i = 0; i < indices.size(); i++) {
                        new_array[i] = array[indices[i]];
                    }
                    return new_array;
                }
            }
            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];
                }
            }
            void merge(const ivec& new_flags, const ivec& indices);
            vec subset(const vec& array, const ivec& indices, ivec& new_indices);
            vec subset(const vec& array, const ivec& indices);
            ivec subset(const ivec& indices);
            Points subset(const Points& input, const ivec& indices);
    };
    class not_implemented_exception: public std::logic_error
    {
        public:
            not_implemented_exception() : std::logic_error("Function not yet implemented") { };
    };

}
#endif