r.grow.distance: add minimum and maximum distance options (#1453)

* r.grow.distance: mindist must be smaller than maxdist
* r.grow.distance: update COPYRIGHT
 * explain result with both minimum_distance and maximum_distance

Co-authored-by: Markus Neteler <neteler@gmail.com>

raster/r.grow.distance/main.c

@@ -6,10 +6,10 @@
  * AUTHOR(S):    Marjorie Larson - CERL
  *               Glynn Clements
  *
- * PURPOSE:      Generates a raster map layer with contiguous areas 
+ * PURPOSE:      Generates a raster map layer with contiguous areas
  *               grown by one cell.
  *
- * COPYRIGHT:    (C) 2006-2013 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2006-2021 by the GRASS Development Team
  *
  *               This program is free software under the GNU General Public
  *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -124,7 +124,7 @@ int main(int argc, char **argv)
     struct GModule *module;
     struct
     {
-	struct Option *in, *dist, *val, *met;
+	struct Option *in, *dist, *val, *met, *min, *max;
     } opt;
     struct
     {
@@ -140,8 +140,9 @@ int main(int argc, char **argv)
     int row, col;
     struct Colors colors;
     struct History hist;
-    DCELL *out_row;
+    DCELL *out_row, *dist_row_max, *val_row_max;
     double scale = 1.0;
+    double mindist, maxdist;
     int invert;
 
     G_gisinit(argv[0]);
@@ -175,6 +176,18 @@ int main(int argc, char **argv)
     opt.met->options = "euclidean,squared,maximum,manhattan,geodesic";
     opt.met->answer = "euclidean";
 
+    opt.min = G_define_option();
+    opt.min->key = "minimum_distance";
+    opt.min->type = TYPE_DOUBLE;
+    opt.min->required = NO;
+    opt.min->description = _("Minimum distance threshold");
+
+    opt.max = G_define_option();
+    opt.max->key = "maximum_distance";
+    opt.max->type = TYPE_DOUBLE;
+    opt.max->required = NO;
+    opt.max->description = _("Maximum distance threshold");
+
     flag.m = G_define_flag();
     flag.m->key = 'm';
     flag.m->description = _("Output distances in meters instead of map units");
@@ -190,7 +203,43 @@ int main(int argc, char **argv)
     dist_name = opt.dist->answer;
     val_name = opt.val->answer;
 
-    if ((invert = flag.n->answer)) {
+    dist_row_max = val_row_max = NULL;
+    mindist = -1;
+    if (opt.min->answer) {
+	if (sscanf(opt.min->answer, "%lf", &mindist) != 1) {
+	    G_warning(_("Invalid %s value '%s', ignoring."),
+		      opt.min->key, opt.min->answer);
+
+	    mindist = -1;
+	}
+    }
+
+    maxdist = -1;
+    if (opt.max->answer) {
+	if (sscanf(opt.max->answer, "%lf", &maxdist) != 1) {
+	    G_warning(_("Invalid %s value '%s', ignoring."),
+		      opt.max->key, opt.max->answer);
+
+	    maxdist = -1;
+	}
+    }
+
+    if (mindist > 0 && maxdist > 0 && mindist >= maxdist) {
+	/* GTC error because given minimum_distance is not smaller than
+	 * given maximum_distance */
+	G_fatal_error(_("'%s' must be smaller than '%s'."),
+	              opt.min->key, opt.max->key);
+    }
+
+    if (mindist > 0 || maxdist > 0) {
+	dist_row_max = Rast_allocate_d_buf();
+	val_row_max = Rast_allocate_d_buf();
+    }
+
+    invert = flag.n->answer;
+    if (mindist <= 0 && maxdist <= 0 && invert) {
+	/* value output for distance to NULL cells makes sense
+	 * if mindist or maxdist is given */
 	if (!dist_name)
 	    G_fatal_error(_("Distance output is required for distance to NULL cells"));
 	if (val_name) {
@@ -224,7 +273,7 @@ int main(int argc, char **argv)
 	G_fatal_error(_("Unknown metric: '%s'"), opt.met->answer);
 
     if (flag.m->answer) {
-	if (window.proj == PROJECTION_LL && 
+	if (window.proj == PROJECTION_LL &&
 	    strcmp(opt.met->answer, "geodesic") != 0) {
 	    G_fatal_error(_("Output distance in meters for lat/lon is only possible with '%s=%s'"),
 	                  opt.met->key, "geodesic");
@@ -265,7 +314,7 @@ int main(int argc, char **argv)
 
     dist_row = Rast_allocate_d_buf();
 
-    if (dist_name && strcmp(opt.met->answer, "euclidean") == 0)
+    if ((mindist > 0 || maxdist > 0 || dist_name) && strcmp(opt.met->answer, "euclidean") == 0)
 	out_row = Rast_allocate_d_buf();
     else
 	out_row = dist_row;
@@ -349,21 +398,55 @@ int main(int argc, char **argv)
 	for (col = ncols - 1; col >= 0; col--)
 	    check(row, col, 1, 0);
 
-	if (dist_name) {
-	    if (out_row != dist_row)
-		for (col = 0; col < ncols; col++)
-		    out_row[col] = sqrt(dist_row[col]);
+	if (mindist > 0 || maxdist > 0) {
+	    /* do not modify dist_row or new_val_row,
+	     * thus copy */
+	    if (out_row != dist_row) {
+		for (col = 0; col < ncols; col++) {
+		    dist_row_max[col] = sqrt(dist_row[col]);
+		}
+	    }
+	    else {
+		for (col = 0; col < ncols; col++) {
+		    dist_row_max[col] = dist_row[col];
+		}
+	    }
+	    for (col = 0; col < ncols; col++) {
+		if (scale != 1.0)
+		    dist_row_max[col] *= scale;
+
+		val_row_max[col] = new_val_row[col];
+
+		if (mindist > 0 && dist_row_max[col] < mindist) {
+		    Rast_set_d_null_value(&dist_row_max[col], 1);
+		    Rast_set_d_null_value(&val_row_max[col], 1);
+		}
+		if (maxdist > 0 && dist_row_max[col] > maxdist) {
+		    Rast_set_d_null_value(&dist_row_max[col], 1);
+		    Rast_set_d_null_value(&val_row_max[col], 1);
+		}
+	    }
+	    if (dist_name)
+		Rast_put_d_row(dist_fd, dist_row_max);
+	    if (val_name)
+		Rast_put_d_row(val_fd, val_row_max);
+	}
+	else {
+	    if (dist_name) {
+		if (out_row != dist_row)
+		    for (col = 0; col < ncols; col++)
+			out_row[col] = sqrt(dist_row[col]);
 
-	    if (scale != 1.0)
-		for (col = 0; col < ncols; col++)
-		    out_row[col] *= scale;
+		if (scale != 1.0)
+		    for (col = 0; col < ncols; col++)
+			out_row[col] *= scale;
 
-	    Rast_put_d_row(dist_fd, out_row);
+		Rast_put_d_row(dist_fd, out_row);
+	    }
+	    if (val_name)
+		Rast_put_d_row(val_fd, new_val_row);
 	}
 
-	if (val_name)
-	    Rast_put_d_row(val_fd, new_val_row);
-
 	swap_rows();
     }
 

raster/r.grow.distance/r.grow.distance.html

@@ -64,6 +64,21 @@ be used only in latitude-longitude locations. It is recommended
 to use it along with the <em>-m</em> flag in order to output 
 distances in meters instead of map units.
 
+<p>
+If <b>minimum_distance</b> is given, all cells with a distance smaller
+than <b>minimum_distance</b> will be set to NULL.
+
+<p>
+If <b>maximum_distance</b> is given, all cells with a distance larger
+than <b>maximum_distance</b> will be set to NULL. The resultant output
+is equivalent to a buffer.
+
+<p>
+If both <b>minimum_distance</b> and <b>maximum_distance</b> are given,
+the result will be similar to a doughnut, a restricted belt for a
+given distance range. All cells outside this distance range will be set
+to NULL.
+
 <h2>EXAMPLES</h2>
 
 <h3>Distance from the streams network</h3>