r.buildvrt: new module to build a GRASS virtual raster (VRT)

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@72763 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.buildvrt/Makefile

@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM  = r.buildvrt
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster/r.buildvrt/link.c

@@ -0,0 +1,134 @@
+#include <stdio.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "proto.h"
+
+void make_cell(const char *output, int maptype)
+{
+    FILE *fp;
+
+    fp = G_fopen_new("cell", output);
+    if (!fp)
+	G_fatal_error(_("Unable to create cell/%s file"), output);
+
+    fclose(fp);
+
+    if (maptype == CELL_TYPE)
+	return;
+
+    fp = G_fopen_new("fcell", output);
+    if (!fp)
+	G_fatal_error(_("Unable to create fcell/%s file"), output);
+
+    fclose(fp);
+}
+
+void make_link(const struct input *inputs, int num_inputs,
+               const char *output)
+{
+    int i;
+    FILE *fp;
+
+    fp = G_fopen_new_misc("cell_misc", "vrt", output);
+    if (!fp)
+	G_fatal_error(_("Unable to create cell_misc/%s/vrt file"), output);
+
+    for (i = 0; i < num_inputs; i++) {
+	const struct input *p = &inputs[i];
+	
+	fprintf(fp, "%s@%s\n", p->name, p->mapset);
+    }
+
+    fclose(fp);
+}
+
+void write_fp_format(const char *output, int maptype)
+{
+    struct Key_Value *key_val;
+    const char *type;
+    FILE *fp;
+
+    if (maptype == CELL_TYPE)
+	return;
+
+    key_val = G_create_key_value();
+
+    type = (maptype == FCELL_TYPE)
+	? "float"
+	: "double";
+    G_set_key_value("type", type, key_val);
+
+    G_set_key_value("byte_order", "xdr", key_val);
+
+    fp = G_fopen_new_misc("cell_misc", "f_format", output);
+    if (!fp)
+	G_fatal_error(_("Unable to create cell_misc/%s/f_format file"), output);
+
+    if (G_fwrite_key_value(fp, key_val) < 0)
+	G_fatal_error(_("Error writing cell_misc/%s/f_format file"), output);
+
+    fclose(fp);
+
+    G_free_key_value(key_val);
+}
+
+void write_fp_quant(const char *output)
+{
+    struct Quant quant;
+
+    Rast_quant_init(&quant);
+    Rast_quant_round(&quant);
+
+    Rast_write_quant(output, G_mapset(), &quant);
+}
+
+void create_map(const struct input *inputs, int num_inputs, const char *output,
+		struct Cell_head *cellhd, int maptype, DCELL dmin, DCELL dmax,
+		int have_stats, struct R_stats *ostats, const char *title)
+{
+    struct History history;
+    struct Categories cats;
+    struct Colors colors;
+
+    Rast_put_cellhd(output, cellhd);
+
+    make_cell(output, maptype);
+
+    make_link(inputs, num_inputs, output);
+
+    if (maptype == CELL_TYPE) {
+	struct Range range;
+	range.min = (CELL)dmin;
+	range.max = (CELL)dmax;
+	range.first_time = 0;
+	Rast_write_range(output, &range);
+    }
+    else {
+	struct FPRange fprange;
+	fprange.min = dmin;
+	fprange.max = dmax;
+	fprange.first_time = 0;
+	Rast_write_fp_range(output, &fprange);
+	write_fp_format(output, maptype);
+	write_fp_quant(output);
+    }
+    G_remove_misc("cell_misc", "stats", output);
+    if (have_stats)
+	Rast_write_rstats(output, ostats);
+
+    G_verbose_message(_("Creating support files for %s"), output);
+    Rast_short_history(output, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(output, &history);
+
+    if (Rast_read_colors(inputs[0].name, inputs[0].mapset, &colors) == 1)
+	Rast_write_colors(output, G_mapset(), &colors);
+    Rast_init_cats(NULL, &cats);
+    Rast_write_cats((char *)output, &cats);
+
+    if (title)
+	Rast_put_cell_title(output, title);
+
+    G_done_msg(_("Link to raster map <%s> created."), output);
+}

raster/r.buildvrt/main.c

@@ -0,0 +1,310 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.buildvrt
+ *               
+ * AUTHOR(S):    Markus Metz, based on r.external
+ *
+ * PURPOSE:      Build a VRT (Virtual Raster) that is a mosaic of the
+ *               list of input raster maps.
+ *
+ * COPYRIGHT:    (C) 2018 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 <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+#include "proto.h"
+
+int cmp_wnd(const void *a, const void *b)
+{
+    struct Cell_head *cellhda = &((struct input *) a)->cellhd;
+    struct Cell_head *cellhdb = &((struct input *) b)->cellhd;
+
+    /* sort from descending N to S, then ascending from W to E */
+    if (cellhda->south > cellhdb->south)
+	return -1;
+    if (cellhda->south < cellhdb->south)
+	return 1;
+    if (cellhda->north > cellhdb->north)
+	return -1;
+    if (cellhda->north < cellhdb->north)
+	return 1;
+    if (cellhda->west < cellhdb->west)
+	return -1;
+    if (cellhda->west > cellhdb->west)
+	return 1;
+    if (cellhda->east < cellhdb->east)
+	return -1;
+    if (cellhda->east > cellhdb->east)
+	return 1;
+
+    return 0;
+}
+
+int main(int argc, char *argv[])
+{
+    const char *output;
+    char *title;
+    struct Cell_head cellhd;
+    struct GModule *module;
+    struct {
+	struct Option *input, *file, *output, *title;
+    } parm;
+    int i, j, num_inputs;
+    struct input *inputs = NULL;
+    char nsresstr[1024], ewresstr[1024];
+    int maptype;
+    struct FPRange fprange;
+    DCELL dmin, dmax;
+    int have_stats;
+    struct R_stats rstats, ostats;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("mosaic"));
+    G_add_keyword(_("virtual raster"));
+    module->description =
+	_("Build a VRT (Virtual Raster) from the list of input raster maps.");
+
+    parm.input = G_define_standard_option(G_OPT_R_INPUTS);
+    parm.input->description = _("Name of input raster files");
+    parm.input->required = NO;
+    parm.input->guisection = _("Input");
+
+    parm.file = G_define_standard_option(G_OPT_F_INPUT);
+    parm.file->key = "file";
+    parm.file->description = _("Input file with one raster map name per line");
+    parm.file->required = NO;
+    parm.file->guisection = _("Input");
+
+    parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+    parm.output->guisection = _("Output");
+    
+    parm.title = G_define_option();
+    parm.title->key = "title";
+    parm.title->key_desc = "phrase";
+    parm.title->type = TYPE_STRING;
+    parm.title->required = NO;
+    parm.title->description = _("Title for resultant raster map");
+    parm.title->guisection = _("Output");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    if (parm.input->answer && parm.file->answer)
+        G_fatal_error(_("%s= and %s= are mutually exclusive"),
+			parm.input->key, parm.file->key);
+ 
+    if (!parm.input->answer && !parm.file->answer)
+        G_fatal_error(_("Please specify %s= or %s="),
+			parm.input->key, parm.file->key);
+
+    /* read input maps from file */
+    if (parm.file->answer) {
+	FILE *in;
+	int max_inputs;
+
+	if (strcmp(parm.file->answer, "-") == 0)
+	    in = stdin;
+	else {
+	    in = fopen(parm.file->answer, "r");
+	    if (!in)
+		G_fatal_error(_("Unable to open input file <%s>"), parm.file->answer);
+	}
+    
+	num_inputs = 0;
+	max_inputs = 0;
+
+	for (;;) {
+	    char buf[GNAME_MAX];
+	    char *name;
+	    const char *mapset;
+	    struct input *p;
+
+	    if (!G_getl2(buf, sizeof(buf), in))
+		break;
+
+	    /* Ignore empty lines */
+	    if (!*buf)
+		continue;
+
+	    name = buf;
+	    if ((mapset = G_find_raster(name, "")) == NULL)
+		G_fatal_error(_("Input raster map <%s> not found"), name);
+
+	    Rast_read_fp_range(name, mapset, &fprange);
+	    dmin = fprange.min;
+	    if (Rast_is_d_null_value(&dmin)) {
+		G_verbose_message(_("Input map <%s@%s> is all NULL, skipping"),
+		                  name, mapset);
+		continue;
+	    }
+
+	    if (num_inputs >= max_inputs) {
+		max_inputs += 100;
+		inputs = G_realloc(inputs, max_inputs * sizeof(struct input));
+	    }
+	    p = &inputs[num_inputs++];
+
+	    p->name = G_store(name);
+            p->mapset = G_store(mapset);
+	    p->maptype = Rast_map_type(p->name, p->mapset);
+	    Rast_get_cellhd(p->name, p->mapset, &(p->cellhd));
+	}
+	fclose(in);
+
+	if (num_inputs < 1)
+	    G_fatal_error(_("No raster map name found in input file"));
+        
+	if (num_inputs == 1)
+	    G_fatal_error(_("Only one raster map name found in input file"));
+    }
+    else {
+    	for (i = 0; parm.input->answers[i]; i++)
+	    ;
+    	num_inputs = i;
+
+    	if (num_inputs < 1)
+	    G_fatal_error(_("Raster map not found"));
+
+	if (num_inputs == 1)
+	    G_fatal_error(_("Only one raster map name found"));
+
+    	inputs = G_malloc(num_inputs * sizeof(struct input));
+
+	j = 0;
+    	for (i = 0; i < num_inputs; i++) {
+	    char *name;
+	    const char *mapset;
+	    struct input *p = &inputs[i];
+
+	    name = parm.input->answers[i];
+	    if ((mapset = G_find_raster(name, "")) == NULL)
+		G_fatal_error(_("Input raster map <%s> not found"), name);
+
+	    Rast_read_fp_range(name, mapset, &fprange);
+	    dmin = fprange.min;
+	    if (Rast_is_d_null_value(&dmin)) {
+		G_verbose_message(_("Input map <%s@%s> is all NULL, skipping"),
+		                  name, mapset);
+		continue;
+	    }
+
+	    p = &inputs[j++];
+
+	    p->name = G_store(name);
+            p->mapset = G_store(mapset);
+	    p->maptype = Rast_map_type(p->name, p->mapset);
+	    Rast_get_cellhd(p->name, p->mapset, &(p->cellhd));
+    	}
+	num_inputs = j;
+    }
+
+    qsort(inputs, num_inputs, sizeof(struct input), cmp_wnd);
+
+    /* check resolution and maptype of input maps */
+    cellhd = inputs[0].cellhd;
+    cellhd.compressed = 0;
+    G_format_resolution(cellhd.ns_res, nsresstr, G_projection());
+    G_format_resolution(cellhd.ew_res, ewresstr, G_projection());
+    maptype = inputs[0].maptype;
+
+    Rast_set_d_null_value(&dmin, 1);
+    Rast_set_d_null_value(&dmax, 1);
+    if (Rast_read_fp_range(inputs[0].name, inputs[0].mapset, &fprange) == 1) {
+	dmin = fprange.min;
+	dmax = fprange.max;
+    }
+    Rast_set_d_null_value(&(ostats.sum), 1);
+    Rast_set_d_null_value(&(ostats.sumsq), 1);
+    ostats.count = 0;
+    have_stats = 1;
+    if (Rast_read_rstats(inputs[0].name, inputs[0].mapset, &rstats) == 1) {
+	ostats.sum = rstats.sum;
+	ostats.sumsq = rstats.sumsq;
+	ostats.count = rstats.count;
+    }
+    else
+	have_stats = 0;
+
+    for (i = 1; i < num_inputs; i++) {
+	char tnsresstr[1024], tewresstr[1024];
+	int tmaptype;
+	struct input *p = &inputs[i];
+
+	G_format_resolution(p->cellhd.ns_res, tnsresstr, G_projection());
+	G_format_resolution(p->cellhd.ew_res, tewresstr, G_projection());
+	tmaptype = p->maptype;
+
+	if (tmaptype != maptype)
+	    G_warning(_("Input maptypes are different"));
+	if (strcmp(nsresstr, tnsresstr) != 0)
+	    G_warning(_("Input ns resolutions are different"));
+	if (strcmp(ewresstr, tewresstr) != 0)
+	    G_warning(_("Input ns resolutions are different"));
+
+	if (cellhd.north < p->cellhd.north)
+	    cellhd.north = p->cellhd.north;
+	if (cellhd.south > p->cellhd.south)
+	    cellhd.south = p->cellhd.south;
+	if (cellhd.east < p->cellhd.east)
+	    cellhd.east = p->cellhd.east;
+	if (cellhd.west > p->cellhd.west)
+	    cellhd.west = p->cellhd.west;
+
+	if (Rast_read_fp_range(p->name, p->mapset, &fprange) == 1) {
+	    if (Rast_is_d_null_value(&dmin)) {
+		dmin = fprange.min;
+		dmax = fprange.max;
+	    }
+	    else {
+		if (dmin > fprange.min)
+		    dmin = fprange.min;
+		if (dmax < fprange.max)
+		    dmax = fprange.max;
+	    }
+	}
+	if (have_stats && 
+	    Rast_read_rstats(inputs[0].name, inputs[0].mapset, &rstats) == 1) {
+	    ostats.sum += rstats.sum;
+	    ostats.sumsq += rstats.sumsq;
+	    ostats.count += rstats.count;
+	}
+	else
+	    have_stats = 0;
+    }
+
+    G_adjust_Cell_head(&cellhd, 0, 0);
+
+    if (maptype == CELL_TYPE)
+	cellhd.format = 3;
+    else
+	cellhd.format = -1;
+
+    output = parm.output->answer;
+
+    title = NULL;
+    if (parm.title->answer) {
+	title = G_store(parm.title->answer);
+	G_strip(title);
+    }
+
+    create_map(inputs, num_inputs, output, &cellhd, maptype, dmin, dmax,
+               have_stats, &ostats, title);
+
+    exit(EXIT_SUCCESS);
+}

raster/r.buildvrt/proto.h

@@ -0,0 +1,23 @@
+#include <grass/raster.h>
+
+#undef MIN
+#undef MAX
+#define MIN(a,b)      ((a) < (b) ? (a) : (b))
+#define MAX(a,b)      ((a) > (b) ? (a) : (b))
+
+struct input
+{
+    const char *name;
+    const char *mapset;
+    int maptype;
+    struct Cell_head cellhd;
+};
+
+/* link.c */
+void make_cell(const char *, int);
+void make_link(const struct input *, int, const char *);
+void write_fp_format(const char *, int);
+void write_fp_quant(const char *);
+void create_map(const struct input *, int, const char *,
+		struct Cell_head *, int, DCELL, DCELL,
+		int, struct R_stats *, const char *);

raster/r.buildvrt/r.buildvrt.html

@@ -0,0 +1,53 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.buildvrt</em> builds a virtual raster (VRT) that is a mosaic of 
+the list of input raster maps. The purpose of such a VRT is to provide 
+fast access to small subsets of the VRT, also with multiple simultaneous 
+read requests.
+
+<h2>NOTES</h2>
+
+<em>r.buildvrt</em> creates a list of raster maps that can be 
+located in different mapsets. The ouput is a read-only link to 
+the original raster maps which is only valid if the original raster 
+maps remain at the originally indicated mapset. A VRT can also be built 
+from raster maps registered with <em>r.external</em>.
+
+<p>
+Reading the whole VRT is slower than reading the equivalent single 
+raster map. Only reading small parts of the VRT provides a performance 
+benefit.
+
+
+<h2>EXAMPLES</h2>
+
+<h3>VRT from a DEM in the North Carolina sample dataset</h3>
+
+<div class="code"><pre>
+# set the region
+g.region raster=elev_state_500m -p
+# higher resolution
+g.region res=50 -p
+# resample the DEM to 50 meter
+r.resamp.interp input=elev_state_500m output=elev_state_50m method=bilinear
+# create tiles
+r.tile input=elev_state_50m output=elev_state_50m_tile_ width=1000 height=1000 overlap=0
+# dump list of tiles to a file
+g.list type=raster pattern=elev_state_50m_tile_* > tilelist
+# build a vrt
+r.buildvrt file=tilelist output=elev_state_50m_vrt
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.tile.html">r.tile</a>,
+<a href="r.patch.html">r.patch</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Markus Metz
+
+<p>
+<i>Last changed: $Date$</i>