Added average_le_quantile,average_ge_quantile operators to libstats a…

…nd r.neighbors module. These are very simple averages for filtered occurrences by quantile.
fpl · Dec 23, 2022 · 6b83795 · 6b83795
1 parent 2b6ab28
commit 6b83795
Show file tree

Hide file tree

Showing 3 changed files with 187 additions and 1 deletion.
diff --git a/include/defs/stats.h b/include/defs/stats.h
@@ -29,6 +29,8 @@ extern stat_func c_perc90;
 extern stat_func c_quant;
 extern stat_func c_skew;
 extern stat_func c_kurt;
+extern stat_func c_ave_ge_quant;
+extern stat_func c_ave_le_quant;
 
 extern stat_func_w w_ave;
 extern stat_func_w w_count;
@@ -49,6 +51,8 @@ extern stat_func_w w_sum;
 extern stat_func_w w_var;
 extern stat_func_w w_skew;
 extern stat_func_w w_kurt;
+extern stat_func_w w_ave_ge_quant;
+extern stat_func_w w_ave_le_quant;
 
 extern int sort_cell(DCELL *, int);
 extern int sort_cell_w(DCELL(*)[2], int);

diff --git a/lib/stats/c_percentile.c b/lib/stats/c_percentile.c
@@ -26,6 +26,92 @@ void c_quant(DCELL * result, DCELL * values, int n, const void *closure)
 	: values[i0] * (i1 - k) + values[i1] * (k - i0);
 }
 
+void c_ave_ge_quant(DCELL * result, DCELL * values, int n, const void *closure)
+{
+    double quant = *(const double *)closure;
+    double k;
+    int i0, i1;
+    DCELL q;
+    DCELL sum;
+    int count;
+    int i;
+
+    n = sort_cell(values, n);
+
+    if (n < 1) {
+	Rast_set_d_null_value(result, 1);
+	return;
+    }
+
+    k = n * quant;
+    i0 = (int)floor(k);
+    i1 = (int)ceil(k);
+
+    q = (i0 == i1)
+	? values[i0]
+	: values[i0] * (i1 - k) + values[i1] * (k - i0);
+
+    sum = 0.0;
+    count = 0;
+    for (i = 0; i < n; i++) {
+        if (Rast_is_d_null_value(&values[i]))
+            continue;
+
+        if (values[i] >= q) {
+            sum += values[i];
+            count++;
+        }
+    }
+
+    if (count == 0)
+	Rast_set_d_null_value(result, 1);
+    else
+	*result = sum / count;
+}
+
+void c_ave_le_quant(DCELL * result, DCELL * values, int n, const void *closure)
+{
+    double quant = *(const double *)closure;
+    double k;
+    int i0, i1;
+    DCELL q;
+    DCELL sum;
+    int count;
+    int i;
+
+    n = sort_cell(values, n);
+
+    if (n < 1) {
+	Rast_set_d_null_value(result, 1);
+	return;
+    }
+
+    k = n * quant;
+    i0 = (int)floor(k);
+    i1 = (int)ceil(k);
+
+    q = (i0 == i1)
+	? values[i0]
+	: values[i0] * (i1 - k) + values[i1] * (k - i0);
+
+    sum = 0.0;
+    count = 0;
+    for (i = 0; i < n; i++) {
+        if (Rast_is_d_null_value(&values[i]))
+            continue;
+
+        if (values[i] <= q) {
+            sum += values[i];
+            count++;
+        }
+    }
+
+    if (count == 0)
+	Rast_set_d_null_value(result, 1);
+    else
+	*result = sum / count;
+}
+
 void c_quart1(DCELL * result, DCELL * values, int n, const void *closure)
 {
     static const double q = 0.25;
@@ -72,6 +158,100 @@ void w_quant(DCELL * result, DCELL(*values)[2], int n, const void *closure)
     *result = values[i][0];
 }
 
+void w_ave_ge_quant(DCELL * result, DCELL(*values)[2], int n, const void *closure)
+{
+    double quant = *(const double *)closure;
+    DCELL total;
+    int i, count;
+    DCELL k;
+    DCELL q;
+    DCELL sum;
+
+    n = sort_cell_w(values, n);
+
+    if (n < 1) {
+	Rast_set_d_null_value(result, 1);
+	return;
+    }
+
+    total = 0.0;
+    for (i = 0; i < n; i++)
+	total += values[i][1];
+
+    k = 0.0;
+    for (i = 0; i < n; i++) {
+	k += values[i][1];
+	if (k >= total * quant)
+	    break;
+    }
+
+    q = values[i][0];
+
+    sum = 0.0;
+    count = 0;
+    for (i = 0; i < n; i++) {
+        if (Rast_is_d_null_value(&values[i][0]))
+            continue;
+
+        if (values[i][0] >= q) {
+            sum += values[i][0];
+            count++;
+        }
+    }
+
+    if (count == 0)
+	Rast_set_d_null_value(result, 1);
+    else
+	*result = sum / count;
+}
+
+void w_ave_le_quant(DCELL * result, DCELL(*values)[2], int n, const void *closure)
+{
+    double quant = *(const double *)closure;
+    DCELL total;
+    int i, count;
+    DCELL k;
+    DCELL q;
+    DCELL sum;
+
+    n = sort_cell_w(values, n);
+
+    if (n < 1) {
+	Rast_set_d_null_value(result, 1);
+	return;
+    }
+
+    total = 0.0;
+    for (i = 0; i < n; i++)
+	total += values[i][1];
+
+    k = 0.0;
+    for (i = 0; i < n; i++) {
+	k += values[i][1];
+	if (k >= total * quant)
+	    break;
+    }
+
+    q = values[i][0];
+
+    sum = 0.0;
+    count = 0;
+    for (i = 0; i < n; i++) {
+        if (Rast_is_d_null_value(&values[i][0]))
+            continue;
+
+        if (values[i][0] <= q) {
+            sum += values[i][0];
+            count++;
+        }
+    }
+
+    if (count == 0)
+	Rast_set_d_null_value(result, 1);
+    else
+	*result = sum / count;
+}
+
 void w_quart1(DCELL * result, DCELL(*values)[2], int n, const void *closure)
 {
     static const double q = 0.25;

diff --git a/raster/r.neighbors/main.c b/raster/r.neighbors/main.c
@@ -72,6 +72,8 @@ static struct menu menu[] = {
     {c_quart3, w_quart3, NO_CATS, 1, 0, T_FLOAT, "quart3", "third quartile"},
     {c_perc90, w_perc90, NO_CATS, 1, 0, T_FLOAT, "perc90", "ninetieth percentile"},
     {c_quant, w_quant, NO_CATS, 1, 0, T_FLOAT, "quantile", "arbitrary quantile"},
+    {c_ave_ge_quant, w_ave_ge_quant, NO_CATS, 1, 0, T_FLOAT, "average_ge_quantile", "average of arbitrary quantile"},
+    {c_ave_le_quant, w_ave_le_quant, NO_CATS, 1, 0, T_FLOAT, "average_le_quantile", "average of arbitrary quantile"},
     {0, 0, 0, 0, 0, 0, 0, 0}
 };
 
@@ -236,7 +238,7 @@ int main(int argc, char *argv[])
     parm.quantile->type = TYPE_DOUBLE;
     parm.quantile->required = NO;
     parm.quantile->multiple = YES;
-    parm.quantile->description = _("Quantile to calculate for method=quantile");
+    parm.quantile->description = _("Quantile to calculate for method=quantile,average_le_quantile,average_ge_quantile");
     parm.quantile->options = "0.0-1.0";
 
     flag.align = G_define_flag();