r.accumulate: Support power-2 directions (#1050)

* r.accumulate: Support power-2 directions * Add a missing break
OSGeo · Apr 5, 2024 · a148e72 · a148e72
1 parent 5c2ee29
commit a148e72
Showing 1 changed file with 27 additions and 8 deletions.
diff --git a/src/raster/r.accumulate/main.c b/src/raster/r.accumulate/main.c
@@ -29,6 +29,7 @@
 #define DIR_UNKNOWN 0
 #define DIR_DEG     1
 #define DIR_DEG45   2
+#define DIR_POW2    3
 
 int main(int argc, char *argv[])
 {
@@ -321,26 +322,35 @@ int main(int argc, char *argv[])
         G_fatal_error(_("Unable to read range file"));
     Rast_get_range_min_max(&dir_range, &dir_min, &dir_max);
     if (dir_max <= 0)
-        G_fatal_error(_("Invalid directions map <%s>"), dir_name);
+        G_fatal_error(_("Invalid direction map <%s>"), dir_name);
 
     dir_format = DIR_UNKNOWN;
     if (strcmp(opt.format->answer, "degree") == 0) {
         if (dir_max > 360)
-            G_fatal_error(_("Directional degrees can not be > 360"));
+            G_fatal_error(_("Directional degrees cannot be > 360"));
         dir_format = DIR_DEG;
     }
     else if (strcmp(opt.format->answer, "45degree") == 0) {
         if (dir_max > 8)
-            G_fatal_error(
-                _("Directional degrees divided by 45 can not be > 8"));
+            G_fatal_error(_("Directional degrees divided by 45 cannot be > 8"));
         dir_format = DIR_DEG45;
     }
+    else if (strcmp(opt.format->answer, "power2") == 0) {
+        if (dir_max > 128)
+            G_fatal_error(_("Powers of 2 cannot be > 128"));
+        dir_format = DIR_POW2;
+    }
     else if (strcmp(opt.format->answer, "auto") == 0) {
         if (dir_max <= 8) {
             dir_format = DIR_DEG45;
             G_important_message(_("Input direction format assumed to be "
                                   "degrees CCW from East divided by 45"));
         }
+        else if (dir_max <= 128) {
+            dir_format = DIR_POW2;
+            G_important_message(_("Input direction format assumed to be "
+                                  "powers of 2 CW from East"));
+        }
         else if (dir_max <= 360) {
             dir_format = DIR_DEG;
             G_important_message(_(
@@ -352,7 +362,7 @@ int main(int argc, char *argv[])
                 dir_name);
     }
     if (dir_format == DIR_UNKNOWN)
-        G_fatal_error(_("Invalid directions format '%s'"), opt.format->answer);
+        G_fatal_error(_("Invalid direction format '%s'"), opt.format->answer);
     /* end of r.path */
 
     /* read outlet coordinates and IDs */
@@ -476,19 +486,28 @@ int main(int argc, char *argv[])
         G_percent(row, nrows, 1);
         dir_buf.c[row] = Rast_allocate_c_buf();
         Rast_get_c_row(dir_fd, dir_buf.c[row], row);
-        if (dir_format == DIR_DEG) {
+        switch (dir_format) {
+        case DIR_DEG:
             for (col = 0; col < ncols; col++) {
                 CELL dir = abs(dir_buf.c[row][col] / 45.0);
 
                 dir_buf.c[row][col] = dir >= NE && dir <= E ? dir : 0;
             }
-        }
-        else {
+            break;
+        case DIR_DEG45:
             for (col = 0; col < ncols; col++) {
                 CELL dir = abs(dir_buf.c[row][col]);
 
                 dir_buf.c[row][col] = dir >= NE && dir <= E ? dir : 0;
             }
+            break;
+        default:
+            for (col = 0; col < ncols; col++) {
+                CELL dir = abs(8 - log2(dir_buf.c[row][col]));
+
+                dir_buf.c[row][col] = dir >= NE && dir <= E ? dir : 0;
+            }
+            break;
         }
     }
     G_percent(1, 1, 1);