[GRASSLIST:6487] Re: use of c_minx and c_maxx in r.series

[Glynn Clements]

vincent at ecovla.nl wrote:

> > > What you describe will give me an output raster file which is in fact
> > > an index referring to my input series of raster files, right? Now what
> > > I want to do is use this result, the index raster, to get values from
> > > comparable but different raster series, e.g. different bands of the
> > > same time series. My problem is not how to define the series or have
> > > them as input using g.mlist, I know how to do that. My problem is that
> > > I want to do e.g. r.series.fromindex input=$series_red
> > > index=$output_of_min_x.
> > 
> > x=1
> > r.mapcalc if(index==$x)
> > 
> > for each step then r.patch all the parts together?
> 
> Hmm it's a solution, but not really what I hoped for: something more
> efficient, and easier to use. Maybe I'll try to create an
> r.series.fromindex...

Currently, there isn't any better approach than the one which Hamish
suggested.

It wouldn't be a lot of work to make a copy of r.series (you only need
main.c) then modify it to do what you need. I've attached a draft
version (untested).

-- 
Glynn Clements <glynn at gclements.plus.com>

-------------- next part --------------
#include <string.h>
#include <stdlib.h>
#include <unistd.h>

#include "gis.h"
#include "glocale.h"

struct input
{
	char *name, *mapset;
	int fd;
	DCELL *buf;
};

int main(int argc, char *argv[])
{
	struct GModule *module;
	struct {
		struct Option *index, *input, *output;
	} parm;
	struct {
		struct Flag *quiet;
	} flag;
	int verbose;
	int i;

	char *idx_name, *idx_mapset;
	int idx_fd;
	CELL *idx_buf;

	int num_inputs;
	struct input *inputs;

	char *out_name;
	int out_fd;
	DCELL *out_buf;

	int nrows, ncols;
	int row, col;

	G_gisinit(argv[0]);

	module = G_define_module();
	module->description =
		_("Selects each cell from one of the input maps "
		  "according to the cell value from the index map.");

	parm.index = G_define_option() ;
	parm.index->key        = "index" ;
	parm.index->type       = TYPE_STRING ;
	parm.index->required   = YES ;
	parm.index->multiple   = NO ;
	parm.index->gisprompt  = "old,cell,raster" ;
	parm.index->description= _("Names of existing raster files") ;

	parm.input = G_define_option() ;
	parm.input->key        = "input" ;
	parm.input->type       = TYPE_STRING ;
	parm.input->required   = YES ;
	parm.input->multiple   = YES ;
	parm.input->gisprompt  = "old,cell,raster" ;
	parm.input->description= _("Names of existing raster files") ;

	parm.output = G_define_option() ;
	parm.output->key        = "output" ;
	parm.output->type       = TYPE_STRING ;
	parm.output->required   = YES ;
	parm.output->gisprompt  = "any,cell,raster" ;
	parm.output->description= _("Name of the new raster file") ;

	flag.quiet = G_define_flag();
	flag.quiet->key = 'q';
	flag.quiet->description = _("Run quietly");

	if (G_parser(argc,argv))
		exit(1);

	verbose = !flag.quiet->answer;

	/* process the index map */

	idx_name = parm.index->answer;
	idx_mapset = G_find_cell2(idx_name,"");
	if (!idx_mapset)
		G_fatal_error("raster file <%s> not found", idx_name);
	idx_fd = G_open_cell_old(idx_name, idx_mapset);
	if (idx_fd < 0)
		G_fatal_error("unable to open index map <%s> in mapset <%s>",
			      idx_name, idx_mapset);
	idx_buf = G_allocate_c_raster_buf();

	/* process the input maps */
	for (i = 0; parm.input->answers[i]; i++)
		;
	num_inputs = i;
	inputs = G_malloc(num_inputs * sizeof(struct input));

	for (i = 0; i < num_inputs; i++)
	{
		struct input *p = &inputs[i];
		p->name = parm.input->answers[i];
		p->mapset = G_find_cell2(p->name,"");
		if (!p->mapset)
			G_fatal_error("raster file <%s> not found", p->name);
		p->fd = G_open_cell_old(p->name, p->mapset);
		if (p->fd < 0)
			G_fatal_error("unable to open input map <%s> in mapset <%s>",
				      p->name, p->mapset);
		p->buf = G_allocate_d_raster_buf();
	}

	/* process the output map */
	out_name = parm.output->answer;

	out_fd = G_open_raster_new(out_name, DCELL_TYPE);
	if (out_fd < 0)
		G_fatal_error("unable to create output map <%s>", out_name);

	out_buf = G_allocate_d_raster_buf();

	/* initialise variables */

	nrows = G_window_rows();
	ncols = G_window_cols();

	/* process the data */
	if (verbose)
		fprintf(stderr, "Percent complete ... ");

	for (row = 0; row < nrows; row++)
	{
		if (verbose)
			G_percent(row, nrows, 2);

		G_get_c_raster_row(idx_fd, idx_buf, row);

		for (i = 0; i < num_inputs; i++)
			G_get_d_raster_row(inputs[i].fd, inputs[i].buf, row);

		for (col = 0; col < ncols; col++)
		{
			CELL idx = idx_buf[col];

			if (G_is_c_null_value(&idx) || idx < 0 || idx >= num_inputs)
				G_set_d_null_value(&out_buf[col], 1);
			else
				memcpy(&out_buf[col], &inputs[idx].buf[col], sizeof(DCELL));
		}

		G_put_d_raster_row(out_fd, out_buf);
	}

	if (verbose)
		G_percent(row, nrows, 2);

	/* close maps */
	G_close_cell(out_fd);

	for (i = 0; i < num_inputs; i++)
		G_close_cell(inputs[i].fd);

	G_close_cell(idx_fd);

	exit(0);
}