[GRASS-user] distance from extent to origin cost surfaces
Dylan Beaudette
dylan.beaudette at gmail.com
Tue Apr 1 12:15:46 EDT 2008
On Mon, Mar 31, 2008 at 4:57 AM, Glynn Clements
<glynn at gclements.plus.com> wrote:
>
> Dylan Beaudette wrote:
>
> > > > This is slightly incorrect -- as we are growing the cost surface to
> > > > compute the perimeter, when we should be shrinking it. However it
> > > > appears that r.grow cannot use a negative radius.
> > >
> > > create a negative mask, grow that, then invert again:
> > >
> > >
> > > r.mapcalc "inv_mask = if(isnull(mask), 1, null())"
> > >
> > > r.grow in=inv_mask out=inv_mask.grown
> > >
> > > r.mapcalc "mask_shrunk = if(isnull(inv_mask.grown), 1, null())"
> >
>
> > I haven't looked at the source for r.grow -- does anyone know how much effort
> > it would take so that it could accept a negative "growth radius" ?
>
> It would make more sense to add a specific "shrink" flag.
>
> I have attached a modified (and untested) version of
> raster/r.grow2/main.c which implements such a flag.
>
> It would be possible to do so with less code duplication (the shrink
> and grow cases are structurally very similar), but would might make
> the code harder to understand. I'm not sure which is preferable.
>
> --
> Glynn Clements <glynn at gclements.plus.com>
>
>
> /****************************************************************************
> *
> * MODULE: r.grow2
> *
> * AUTHOR(S): Marjorie Larson - CERL
> * Glynn Clements
> *
> * PURPOSE: Generates a raster map layer with contiguous areas
> * grown by one cell.
> *
> * COPYRIGHT: (C) 2006 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
> * for details.
> *
> ***************************************************************************/
>
> #include <stdio.h>
> #include <stdlib.h>
> #include <string.h>
> #include <grass/gis.h>
> #include <grass/glocale.h>
>
>
> #define MAX(a, b) ((a) > (b) ? (a) : (b))
>
> static int size;
> static int count;
> static int (*neighbors)[2];
>
> typedef int metric_fn(int, int);
>
>
> static int distance_euclidian_squared(int dx, int dy)
> {
> return dx * dx + dy * dy;
> }
>
> static int distance_maximum(int dx, int dy)
> {
> return MAX(abs(dx), abs(dy));
> }
>
> static int distance_manhattan(int dx, int dy)
> {
> return abs(dx) + abs(dy);
> }
>
> static void setup_neighbors(double radius, int limit, metric_fn *dist)
> {
> int i, dx, dy;
> int n;
>
> size = (int) radius;
>
> n = size * 2 + 1;
>
> neighbors = G_malloc(n * n * 2 * sizeof(int));
>
> count = 0;
>
> for (i = 1; i <= limit; i++)
> {
> for (dy = -size; dy <= size; dy++)
> {
> for (dx = -size; dx <= size; dx++)
> {
> if ((*dist)(dx, dy) != i)
> continue;
>
> neighbors[count][0] = dx;
> neighbors[count][1] = dy;
> count++;
> }
> }
> }
> }
>
> static void setup_neighbors_euclidian(double radius)
> {
> int r2 = (int) (radius * radius);
>
> setup_neighbors(radius, r2, distance_euclidian_squared);
> }
>
> static void setup_neighbors_maximum(double radius)
> {
> setup_neighbors(radius, (int) radius, distance_maximum);
> }
>
> static void setup_neighbors_manhattan(double radius)
> {
> setup_neighbors(radius, (int) radius, distance_manhattan);
> }
>
> int main(int argc, char **argv)
> {
> struct GModule *module;
> struct {
> struct Option *in, *out, *rad, *met, *old, *new;
> } opt;
> struct {
> struct Flag *q, *s;
> } flag;
> struct Colors colr;
> struct Categories cats;
> int verbose;
> int colrfile;
> char *in_name;
> char *out_name;
> double radius;
> int oldval = 0;
> int newval = 0;
> char *mapset;
> RASTER_MAP_TYPE type;
> int in_fd;
> int out_fd;
> DCELL **in_rows;
> DCELL *out_row;
> int nrows, row;
> int ncols, col;
> int shrink;
>
> G_gisinit(argv[0]);
>
> module = G_define_module();
> module->keywords = _("raster");
> module->description =
> _("Generates a raster map layer "
> "with contiguous areas grown by one cell.");
>
> opt.in = G_define_standard_option(G_OPT_R_INPUT);
>
> opt.out = G_define_standard_option(G_OPT_R_OUTPUT);
>
> opt.rad = G_define_option();
> opt.rad->key = "radius";
> opt.rad->type = TYPE_DOUBLE;
> opt.rad->required = NO;
> opt.rad->description = _("Radius of buffer in raster cells");
> opt.rad->answer = "1.01";
>
> opt.met = G_define_option();
> opt.met->key = "metric";
> opt.met->type = TYPE_STRING;
> opt.met->required = NO;
> opt.met->description = _("Metric");
> opt.met->options = "euclidian,maximum,manhattan";
> opt.met->answer = "euclidian";
>
> opt.old = G_define_option();
> opt.old->key = "old";
> opt.old->type = TYPE_INTEGER;
> opt.old->required = NO;
> opt.old->description = _("Value to write for input cells which are non-NULL (-1 => NULL)");
>
> opt.new = G_define_option();
> opt.new->key = "new";
> opt.new->type = TYPE_INTEGER;
> opt.new->required = NO;
> opt.new->description = _("Value to write for \"grown\" cells");
>
> flag.q = G_define_flag();
> flag.q->key = 'q';
> flag.q->description = _("Quiet");
>
> flag.s = G_define_flag();
> flag.s->key = 's';
> flag.s->description = _("Shrink (grow null areas)");
>
> if (G_parser(argc, argv))
> exit(EXIT_FAILURE);
>
> in_name = opt.in->answer;
> out_name = opt.out->answer;
>
> radius = atof(opt.rad->answer);
>
> if (opt.old->answer)
> oldval = atoi(opt.old->answer);
>
> if (opt.new->answer)
> newval = atoi(opt.new->answer);
>
> verbose = !flag.q->answer;
>
> shrink = flag.s->answer ? 1 : 0;
>
> mapset = G_find_cell(in_name, "");
> if (!mapset)
> G_fatal_error(_("Raster map <%s> not found"), in_name);
>
> nrows = G_window_rows();
> ncols = G_window_cols();
>
> if (strcmp(opt.met->answer, "euclidian") == 0)
> setup_neighbors_euclidian(radius);
> else if (strcmp(opt.met->answer, "maximum") == 0)
> setup_neighbors_maximum(radius);
> else if (strcmp(opt.met->answer, "manhattan") == 0)
> setup_neighbors_manhattan(radius);
> else
> G_fatal_error(_("Unknown metric: [%s]."), opt.met->answer);
>
> in_fd = G_open_cell_old(in_name, mapset);
> if (in_fd < 0)
> G_fatal_error(_("Unable to open raster map <%s>"), in_name);
>
> type = G_get_raster_map_type(in_fd);
>
> out_fd = G_open_raster_new(out_name, type);
> if (out_fd < 0)
> G_fatal_error(_("Unable to create raster map <%s>"), out_name);
>
> if (G_read_cats(in_name, mapset, &cats) == -1)
> {
> G_warning(_("Error reading category file for <%s>"), in_name);
> G_init_cats(0, "", &cats);
> }
>
> if (G_read_colors(in_name, mapset, &colr) == -1)
> {
> G_warning(_("Error in reading color file for <%s>"), in_name);
> colrfile = 0;
> }
> else
> colrfile = 1;
>
> if (opt.old->answer && oldval >= 0)
> G_set_cat(oldval, "original cells", &cats);
>
> if (opt.new->answer)
> G_set_cat(newval, "grown cells", &cats);
>
> in_rows = G_malloc((size * 2 + 1) * sizeof(DCELL *));
>
> for (row = 0; row <= size * 2; row++)
> in_rows[row] = G_allocate_d_raster_buf();
>
> out_row = G_allocate_d_raster_buf();
>
> for (row = 0; row < size; row++)
> G_get_d_raster_row(in_fd, in_rows[size + row], row);
>
> for (row = 0; row < nrows; row++)
> {
> DCELL *tmp;
> int i;
>
> if (row + size < nrows)
> G_get_d_raster_row(in_fd, in_rows[size * 2], row + size);
>
> for (col = 0; col < ncols; col++)
> {
> DCELL *c = &in_rows[size][col];
>
> if (shrink)
> {
> if (G_is_d_null_value(c))
> {
> G_set_d_null_value(&out_row[col], 1);
> continue;
> }
>
> for (i = 0; i < count; i++)
> {
> int dx = neighbors[i][0];
> int dy = neighbors[i][1];
> int x = col + dx;
> int y = row + dy;
>
> if (x < 0 || x >= ncols || y < 0 || y >= nrows)
> continue;
>
> c = &in_rows[size + dy][x];
>
> if (G_is_d_null_value(c))
> {
> G_set_d_null_value(&out_row[col], 1);
> break;
> }
> }
>
> if (i == count)
> out_row[col] = *c;
> }
> else
> {
> if (!G_is_d_null_value(c))
> {
> if (opt.old->answer)
> {
> if (oldval < 0)
> G_set_d_null_value(&out_row[col], 1);
> else
> out_row[col] = oldval;
> }
> else
> out_row[col] = *c;
>
> continue;
> }
>
> for (i = 0; i < count; i++)
> {
> int dx = neighbors[i][0];
> int dy = neighbors[i][1];
> int x = col + dx;
> int y = row + dy;
>
> if (x < 0 || x >= ncols || y < 0 || y >= nrows)
> continue;
>
> c = &in_rows[size + dy][x];
>
> if (!G_is_d_null_value(c))
> {
> out_row[col] = opt.new->answer
> ? newval
> : *c;
> break;
> }
> }
>
> if (i == count)
> G_set_d_null_value(&out_row[col], 1);
> }
> }
>
> G_put_d_raster_row(out_fd, out_row);
>
> if (verbose)
> G_percent(row, nrows, 2);
>
> tmp = in_rows[0];
> for (i = 0; i < size * 2; i++)
> in_rows[i] = in_rows[i + 1];
> in_rows[size * 2] = tmp;
> }
>
> if (verbose)
> G_percent(row, nrows, 2);
>
> G_close_cell(in_fd);
> G_close_cell(out_fd);
>
> if (G_write_cats(out_name, &cats) == -1)
> G_warning(_("Error writing category file for <%s>"), out_name);
>
> if (colrfile)
> if (G_write_colors(out_name, G_mapset(), &colr) == -1)
> G_warning(_("Error writing color file for <%s>"), out_name);
>
> return EXIT_SUCCESS;
> }
>
>
>
Thanks for suggesting this Glynn. You are quite correct in that there
is a delicate balance between transparency and efficiency. Before I
make a suggestion, I need a little help understanding this little bit
from the original code:
if (!G_is_d_null_value(c))
{
out_row[col] = opt.new->answer
? newval
: *c;
break;
}
what exactly is this stanza doing, and why does it contain newval
*and* opt.new->answer ??
Either way, I think that the ability to grow null regions would be a
very useful option-- similar to mask "eroding" functions in R.
If it would help I can submit a feature request for this.
Cheers,
Dylan
More information about the grass-user
mailing list