Merge pull request #28 from metno/20231128_sct_dual_for_dataset

added sct_dual method for dataset
metno · Nov 29, 2023 · 74791b2 · 74791b2
2 parents 83bb3bc + a9c5746
commit 74791b2
Show file tree

Hide file tree

Showing 3 changed files with 67 additions and 3 deletions.
diff --git a/include/titanlib.h b/include/titanlib.h
@@ -612,9 +612,27 @@ namespace titanlib {
              */
             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 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 sct_dual(const vec& event_thresholds, 
+                                ConditionType condition, 
+                                int num_min, int num_max,
+                                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,
+                                const ivec& obs_to_check=ivec(), 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));

diff --git a/src/dataset.cpp b/src/dataset.cpp
@@ -51,6 +51,28 @@ void titanlib::Dataset::sct(int num_min, int num_max, float inner_radius, float
     merge(new_flags, indices);
 }
 
+void titanlib::Dataset::sct_dual(const vec& event_thresholds, 
+                                ConditionType condition, 
+                                int num_min, int num_max,
+                                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,
+                                const ivec& obs_to_check,
+                                const ivec& indices){
+    // Only keep flags that are not flagged
+    ivec new_flags = titanlib::sct_dual(get_unflagged_points(), get_unflagged(values), get_unflagged(obs_to_check),
+                                        get_unflagged(event_thresholds), condition,
+                                        num_min, num_max, inner_radius, outer_radius, num_iterations, 
+                                        min_horizontal_scale, max_horizontal_scale,
+                                        kth_closest_obs_horizontal_scale, vertical_scale,
+                                        get_unflagged(test_thresholds), debug);
+    merge(new_flags, indices);
+}
+
 void titanlib::Dataset::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 new_flags = titanlib::buddy_check(get_unflagged_points(), get_unflagged(values), get_unflagged(radius), get_unflagged(num_min), threshold, max_elev_diff, elev_gradient, min_std, num_iterations, get_unflagged(obs_to_check));
     merge(new_flags, indices);

diff --git a/tests/dataset_test.py b/tests/dataset_test.py
@@ -70,6 +70,30 @@ def test_buddy_check(self):
                 max_elev_diff, elev_gradient, min_std, num_iterations)
         np.testing.assert_array_equal(dataset.flags, [0]*8 + [1]*2)
 
+    def test_isolation_check(self):
+        lats = [60, 60, 60, 60]
+        lons = [10, 10.1, 10.2, 12]
+        elevs = [0]*4
+        points = titanlib.Points(lats, lons, elevs)
+        values = [0, 0.1, 1, 2]
+
+        num_min = [1, 1, 1, 1]
+        radii = [100000] * 4
+
+        # Last value should be flagged
+        dataset = titanlib.Dataset(points, values)
+        # Preflag one of the observations
+        dataset.external_check([0, 1, 0, 0])
+        dataset.isolation_check(num_min, radii)
+        np.testing.assert_array_equal(dataset.flags, [0, 1, 0, 1])
+
+        # Check that using a single value for radii also works
+        radii = [100000]
+        dataset = titanlib.Dataset(points, values)
+        dataset.external_check([0, 1, 0, 0])
+        dataset.isolation_check(num_min, radii)
+        np.testing.assert_array_equal(dataset.flags, [0, 1, 0, 1])
+
     def test_get_values(self):
         N = 10
         lats = [60]*N