[GRASS-SVN] r49960 - in grass/trunk/raster: . r.series.interpol
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Dec 29 08:15:50 EST 2011
Author: huhabla
Date: 2011-12-29 05:15:50 -0800 (Thu, 29 Dec 2011)
New Revision: 49960
Added:
grass/trunk/raster/r.series.interpol/
grass/trunk/raster/r.series.interpol/Makefile
grass/trunk/raster/r.series.interpol/main.c
grass/trunk/raster/r.series.interpol/r.series.interpol.html
grass/trunk/raster/r.series.interpol/test.r.series.interpol.sh
grass/trunk/raster/r.series.interpol/test_1_prec_2.ref
grass/trunk/raster/r.series.interpol/test_1_prec_3.ref
grass/trunk/raster/r.series.interpol/test_1_prec_4.ref
grass/trunk/raster/r.series.interpol/test_2_prec_2.ref
grass/trunk/raster/r.series.interpol/test_2_prec_3.ref
grass/trunk/raster/r.series.interpol/test_2_prec_4.ref
Log:
New module for raster map interpolation located temporal or spatial in between existing raster maps
Added: grass/trunk/raster/r.series.interpol/Makefile
===================================================================
--- grass/trunk/raster/r.series.interpol/Makefile (rev 0)
+++ grass/trunk/raster/r.series.interpol/Makefile 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.series.interpol
+
+LIBES = $(STATSLIB) $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(STATSDEP) $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass/trunk/raster/r.series.interpol/main.c
===================================================================
--- grass/trunk/raster/r.series.interpol/main.c (rev 0)
+++ grass/trunk/raster/r.series.interpol/main.c 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,337 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.series.interpol
+ * AUTHOR(S): Soeren Gebbert <soerengebbert googlemail.com>
+ * Code is based on r.series from Glynn Clements <glynn gclements.plus.com>
+ *
+ * PURPOSE: Interpolate raster maps located (temporal or spatial) in between input raster maps at specific sampling positions
+ * COPYRIGHT: (C) 2011-2012 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 <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <math.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include <grass/stats.h>
+
+#define LINEAR_INTERPOLATION 1
+#define QUADRATIC_INTERPOLATION 2
+#define CUBIC_INTERPOLATION 3
+
+struct map_store
+{
+ const char *name;
+ double pos;
+ DCELL *buf;
+ int fd;
+};
+
+double linear_position[2] = {0.0, 1.0};
+double quadratic_position[3] = {0.0, 0.5, 1.0};
+double cubic_position[4] = {0.0, 1.0/3.0, 2.0/3.0, 1.0};
+
+static void linear_interpolation(struct map_store *inputs, struct map_store *out, int ncols);
+static int start_interpolation(struct map_store *inputs, int num_inputs, int interpolation_method, struct map_store *out);
+
+int main(int argc, char *argv[])
+{
+ struct GModule *module;
+ struct
+ {
+ struct Option *input, *file, *output, *sampoint, *method;
+ } parm;
+ int i;
+ int num_outputs;
+ int num_inputs;
+ int num_sampoint;
+ struct map_store *inputs = NULL;
+ struct map_store *outputs = NULL;
+ int interpol_method = LINEAR_INTERPOLATION;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("series"));
+ G_add_keyword(_("interpolation"));
+ module->description =
+ _("Interpolate raster maps located (temporal or spatial) "
+ "in between input raster maps at specific sampling positions.");
+
+ parm.input = G_define_standard_option(G_OPT_R_INPUTS);
+
+ parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.output->multiple = YES;
+ parm.output->required = NO;
+
+ parm.sampoint = G_define_option();
+ parm.sampoint->key = "sampoint";
+ parm.sampoint->type = TYPE_DOUBLE;
+ parm.sampoint->required = NO;
+ parm.sampoint->description = _("Sampling point for each input map,"
+ " the point must in between the interval (0;1)");
+ parm.sampoint->multiple = YES;
+
+ parm.file = G_define_standard_option(G_OPT_F_INPUT);
+ parm.file->key = "file";
+ parm.file->description = _("Input file with output a raster map name and sample point per line,"
+ " field separator between name and sample point is |");
+ parm.file->required = NO;
+
+ parm.method = G_define_option();
+ parm.method->key = "method";
+ parm.method->type = TYPE_STRING;
+ parm.method->required = NO;
+ parm.method->options = "linear";
+ parm.method->answer = "linear";
+ parm.method->description = _("Interpolation method, currently only linear interpolation is supported");
+ parm.method->multiple = NO;
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ if (parm.output->answer && parm.file->answer)
+ G_fatal_error(_("input= and file= are mutually exclusive"));
+
+ if (!parm.input->answer && !parm.file->answer)
+ G_fatal_error(_("Please specify input= or file="));
+
+ if (parm.output->answer && !parm.sampoint->answer)
+ G_fatal_error(_("Please specify input= and sampoint="));
+
+ if(G_strncasecmp(parm.method->answer, "linear", 6))
+ interpol_method = LINEAR_INTERPOLATION;
+
+ /* process the input maps */
+ for (i = 0; parm.input->answers[i]; i++)
+ ;
+ num_inputs = i;
+
+ if (num_inputs < 1)
+ G_fatal_error(_("No input raster map not found"));
+
+ if(interpol_method == LINEAR_INTERPOLATION)
+ if (num_inputs != 2)
+ G_fatal_error(_("You need to specify two input maps for linear interpolation"));
+
+ if(interpol_method == QUADRATIC_INTERPOLATION)
+ if (num_inputs != 3)
+ G_fatal_error(_("You need to specify three input maps for quadratic interpolation"));
+
+ if(interpol_method == CUBIC_INTERPOLATION)
+ if (num_inputs != 4)
+ G_fatal_error(_("You need to specify 4 input maps for linear interpolation"));
+
+
+ inputs = G_calloc(num_inputs, sizeof(struct map_store));
+
+ for (i = 0; i < num_inputs; i++) {
+ struct map_store *in = &inputs[i];
+
+ in->name = parm.input->answers[i];
+
+ if(interpol_method == LINEAR_INTERPOLATION)
+ in->pos = linear_position[i];
+ if(interpol_method == QUADRATIC_INTERPOLATION)
+ in->pos = quadratic_position[i];
+ if(interpol_method == CUBIC_INTERPOLATION)
+ in->pos = cubic_position[i];
+
+ in->buf = Rast_allocate_d_buf();
+ in->fd = Rast_open_old(in->name, "");
+ G_verbose_message(_("Reading input raster map <%s> at sample point %g..."), in->name, in->pos);
+ }
+
+ /* process the output maps from the file */
+ if (parm.file->answer) {
+ FILE *in;
+ int max_outputs;
+ double left = 0.0;
+ double right = 1.0;
+
+ if(interpol_method == LINEAR_INTERPOLATION) {
+ left = linear_position[0];
+ right = linear_position[1];
+ }
+
+ in = fopen(parm.file->answer, "r");
+ if (!in)
+ G_fatal_error(_("Unable to open input file <%s>"), parm.file->answer);
+
+ num_outputs = 0;
+ max_outputs = 0;
+
+ for (;;) {
+ char buf[GNAME_MAX + 50]; /* Name and position */
+ char tok_buf[GNAME_MAX + 50]; /* Name and position */
+ char *name;
+ int ntokens;
+ char **tokens;
+ struct map_store *p;
+ double pos = -1;
+
+ if (!G_getl2(buf, sizeof(buf), in))
+ break;
+
+ strcpy(tok_buf, buf);
+ tokens = G_tokenize(tok_buf, "|");
+ ntokens = G_number_of_tokens(tokens);
+
+ if(ntokens > 1) {
+ name = G_chop(tokens[0]);
+ pos = atof(G_chop(tokens[1]));
+ if(pos < left || pos > right)
+ G_fatal_error(_("Wrong sampling point for output map <%s> in file <%s> "
+ "near line %i, sampling point must be in between (%g:%g) not: %g"),
+ name, parm.file->answer, num_outputs, left, right, pos);
+ } else {
+ name = G_chop(buf);
+ }
+
+ /* Ignore empty lines */
+ if (!*name)
+ continue;
+
+ if(pos == -1)
+ G_fatal_error(_("Missing sampling point for output map <%s>"
+ " in file <%s> near line %i"),
+ name, parm.file->answer, num_outputs);
+
+
+ if (num_outputs >= max_outputs) {
+ max_outputs += 100;
+ outputs = G_realloc(outputs, max_outputs * sizeof(struct map_store));
+ }
+ p = &outputs[num_outputs++];
+
+ p->name = G_store(name);
+ p->pos = pos;
+ p->fd = -1;
+ p->buf = NULL;
+ }
+
+ if (num_outputs < 1)
+ G_fatal_error(_("No raster map name found in file <%s>"), parm.file->answer);
+
+ fclose(in);
+ }
+ else {
+ for (i = 0; parm.output->answers[i]; i++)
+ ;
+ num_outputs = i;
+
+ if (num_outputs < 1)
+ G_fatal_error(_("No output raster map not found"));
+
+ for (i = 0; parm.sampoint->answers[i]; i++)
+ ;
+ num_sampoint = i;
+
+ if (num_sampoint != num_outputs)
+ G_fatal_error(_("input= and sampoint= must have the same number of values"));
+
+ outputs = G_malloc(num_outputs * sizeof(struct map_store));
+
+ for (i = 0; i < num_outputs; i++) {
+ struct map_store *p = &outputs[i];
+
+ p->name = parm.output->answers[i];
+ p->pos = (DCELL)atof(parm.sampoint->answers[i]);
+ p->fd = -1;
+ p->buf = NULL;
+ }
+ }
+
+ /* Start the interpolation for each output map */
+ for (i = 0; i < num_outputs; i++) {
+ G_verbose_message(_("Processing output raster map <%s> at sample point %g..."), outputs[i].name, outputs[i].pos);
+ start_interpolation(inputs, num_inputs, interpol_method, &outputs[i]);
+ }
+
+ /* Close input maps */
+ for (i = 0; i < num_inputs; i++)
+ Rast_close(inputs[i].fd);
+
+ exit(EXIT_SUCCESS);
+}
+
+int start_interpolation(struct map_store *inputs, int num_inputs, int interpol_method, struct map_store *out) {
+
+ struct History history;
+ int row, i;
+ int nrows = Rast_window_rows();
+ int ncols = Rast_window_cols();
+
+ out->fd = Rast_open_new(out->name, DCELL_TYPE);
+ out->buf = Rast_allocate_d_buf();
+
+ /* process the data */
+ G_verbose_message(_("Percent complete..."));
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+
+ for (i = 0; i < num_inputs; i++)
+ Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
+
+ if(interpol_method == LINEAR_INTERPOLATION)
+ linear_interpolation(inputs, out, ncols);
+ /* Add new interpolation methods here */
+
+ Rast_put_d_row(out->fd, out->buf);
+ }
+
+ G_percent(row, nrows, 2);
+
+ Rast_close(out->fd);
+
+ Rast_short_history(out->name, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(out->name, &history);
+
+ G_free(out->buf);
+
+ return 1;
+}
+
+/* linear function v = (1 - pos/dist) * u1 + pos/dist * u2
+ *
+ * v -> The value of the output map
+ * pos -> The normalized position of the output map (0;1)
+ * u1 -> The value of the left input map
+ * u2 -> The value of the right input map
+ * dist -> The distance between the position of u1 and u2
+ *
+ * */
+void linear_interpolation(struct map_store *inputs, struct map_store *out, int ncols)
+{
+ DCELL v;
+ DCELL u1;
+ DCELL u2;
+ DCELL dist;
+ int col;
+
+ for (col = 0; col < ncols; col++) {
+ u1 = inputs[0].buf[col];
+ u2 = inputs[1].buf[col];
+ dist = fabs(inputs[1].pos - inputs[0].pos);
+
+
+ if(Rast_is_d_null_value(&u1) || Rast_is_d_null_value(&u2)) {
+ Rast_set_d_null_value(&v, 1);
+ } else {
+ v = (1 - (out->pos/dist)) * u1 + (out->pos/dist) * u2;
+ }
+ out->buf[col] = v;
+ }
+ return;
+}
Added: grass/trunk/raster/r.series.interpol/r.series.interpol.html
===================================================================
--- grass/trunk/raster/r.series.interpol/r.series.interpol.html (rev 0)
+++ grass/trunk/raster/r.series.interpol/r.series.interpol.html 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,60 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.series.interpol</em>
+interpolates raster maps which are located temporal or spatial in between existing input raster maps.
+The interpolation is performed at specific sampling positions. The sampling position for each output map must be specified.
+Sampling positions are always located in the normalized intervall (0;1).
+The following interpolation methods are supported.
+<ul>
+ <li>linear: Linear interpolation. Two input maps are required. The first input map is located at position 0.0, the second input map at position 1.0.
+ </ul>
+
+<h2>EXAMPLES</h2>
+<p>Interpolate linear three new maps at 3 sampling positions in between the interval (0;1)
+
+First prepare the input maps:
+<br>
+<div class="code"><pre>
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
+
+r.mapcalc --o expr="prec_1 = 100"
+r.mapcalc --o expr="prec_5 = 500"
+</pre></div>
+
+<p>Interpolate
+
+<div class="code"><pre>
+r.series.interpol --o --v input=prec_1,prec_5 output=prec_2,prec_3,prec_4 sampoint=0.25,0.5,0.75 method=linear
+</pre></div>
+
+<p>Interpolate using the file option:
+First prepare the input file:
+<br>
+<div class="code"><pre>
+TMP_FILE=`g.tempfile pid=1`
+cat > "${TMP_FILE}" << EOF
+prec_2|0.25
+prec_3|0.5
+prec_4|0.75
+EOF
+</pre></div>
+
+<p>Interpolate
+
+<div class="code"><pre>
+r.series.interpol --o --v input=prec_1,prec_5 file="${TMP_FILE}" method=linear
+</pre></div>
+
+The created maps will have the values 200, 300 and 400.
+
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.series.html">r.series</a></em>
+<em><a href="g.region.html">g.region</a></em>
+
+<h2>AUTHOR</h2>
+
+Sören Gebbert
+
+<p><i>Last changed: $Date: 2011-12-28 15:26:22 +0100 (Mi, 28 Dez 2011) $</i>
Added: grass/trunk/raster/r.series.interpol/test.r.series.interpol.sh
===================================================================
--- grass/trunk/raster/r.series.interpol/test.r.series.interpol.sh (rev 0)
+++ grass/trunk/raster/r.series.interpol/test.r.series.interpol.sh 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,39 @@
+# Test r.series.interpol
+# We need to set a specific region in the
+# @preprocess step of this test. We generate
+# raster maps with r.mapcalc
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
+
+r.mapcalc --o expr="prec_1 = 100"
+r.mapcalc --o expr="prec_5 = 500"
+
+TMP_FILE=`g.tempfile pid=1`
+# We create an input files containing empty lines. However, r.series.interpol should process the
+# valid raster maps listed in the files.
+
+cat > "${TMP_FILE}" << EOF
+
+prec_2|0.25
+
+
+prec_3|0.5
+
+prec_4|0.75
+
+
+EOF
+
+# The first @test with map input and @precision=3
+r.series.interpol --o --v input=prec_1,prec_5 output=prec_2,prec_3,prec_4 sampoint=0.25,0.5,0.75 method=linear
+
+#r.out.ascii --o input=prec_2 output=test_1_prec_2.ref dp=3
+#r.out.ascii --o input=prec_3 output=test_1_prec_3.ref dp=3
+#r.out.ascii --o input=prec_4 output=test_1_prec_4.ref dp=3
+
+# The second @test with file input and @precision=3
+r.series.interpol --o --v input=prec_1,prec_5 file="${TMP_FILE}" method=linear
+
+#r.out.ascii --o input=prec_2 output=test_2_prec_2.ref dp=3
+#r.out.ascii --o input=prec_3 output=test_2_prec_3.ref dp=3
+#r.out.ascii --o input=prec_4 output=test_2_prec_4.ref dp=3
Added: grass/trunk/raster/r.series.interpol/test_1_prec_2.ref
===================================================================
--- grass/trunk/raster/r.series.interpol/test_1_prec_2.ref (rev 0)
+++ grass/trunk/raster/r.series.interpol/test_1_prec_2.ref 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
Added: grass/trunk/raster/r.series.interpol/test_1_prec_3.ref
===================================================================
--- grass/trunk/raster/r.series.interpol/test_1_prec_3.ref (rev 0)
+++ grass/trunk/raster/r.series.interpol/test_1_prec_3.ref 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
Added: grass/trunk/raster/r.series.interpol/test_1_prec_4.ref
===================================================================
--- grass/trunk/raster/r.series.interpol/test_1_prec_4.ref (rev 0)
+++ grass/trunk/raster/r.series.interpol/test_1_prec_4.ref 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
Added: grass/trunk/raster/r.series.interpol/test_2_prec_2.ref
===================================================================
--- grass/trunk/raster/r.series.interpol/test_2_prec_2.ref (rev 0)
+++ grass/trunk/raster/r.series.interpol/test_2_prec_2.ref 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
+200 200 200 200 200 200 200 200 200 200 200 200
Added: grass/trunk/raster/r.series.interpol/test_2_prec_3.ref
===================================================================
--- grass/trunk/raster/r.series.interpol/test_2_prec_3.ref (rev 0)
+++ grass/trunk/raster/r.series.interpol/test_2_prec_3.ref 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
+300 300 300 300 300 300 300 300 300 300 300 300
Added: grass/trunk/raster/r.series.interpol/test_2_prec_4.ref
===================================================================
--- grass/trunk/raster/r.series.interpol/test_2_prec_4.ref (rev 0)
+++ grass/trunk/raster/r.series.interpol/test_2_prec_4.ref 2011-12-29 13:15:50 UTC (rev 49960)
@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
+400 400 400 400 400 400 400 400 400 400 400 400
More information about the grass-commit
mailing list