1
2 /****************************************************************************
3 *
4 * MODULE: r.series.interp
5 * AUTHOR(S): Soeren Gebbert <soerengebbert googlemail.com>
6 * Code is based on r.series from Glynn Clements <glynn gclements.plus.com>
7 *
8 * PURPOSE: Interpolate raster maps located (temporal or spatial) in between input raster maps at specific sampling positions
9 * COPYRIGHT: (C) 2011-2012 by the GRASS Development Team
10 *
11 * This program is free software under the GNU General Public
12 * License (>=v2). Read the infile COPYING that comes with GRASS
13 * for details.
14 *
15 *****************************************************************************/
16 #include <string.h>
17 #include <stdlib.h>
18 #include <unistd.h>
19 #include <math.h>
20
21 #include <grass/gis.h>
22 #include <grass/raster.h>
23 #include <grass/glocale.h>
24 #include <grass/stats.h>
25
26 #define LINEAR_INTERPOLATION 1
27 #define SPLINE_INTERPOLATION 2
28
29 struct map_store
30 {
31 const char *name;
32 double pos;
33 DCELL *buf;
34 int fd;
35 int has_run;
36 };
37
38 void selection_sort(struct map_store **array, int num);
39 static struct map_store *get_parameter_input(const char *type, char **map_names, char **positions, char *file, int *number_of_maps);
40 static void linear_interpolation(struct map_store **inp, int num_inputs, struct map_store **outp, int num_outputs);
41 static void interpolate_row_linear(struct map_store *left, struct map_store *right, struct map_store *out, int ncols);
42 static void start_interpolation(struct map_store *inputs, int num_inputs, struct map_store *outputs, int num_outputs, int interpol_method);
43
44 /* *************************************************************** */
45 /* *************************************************************** */
46 /* *************************************************************** */
47
main(int argc,char * argv[])48 int main(int argc, char *argv[])
49 {
50 struct GModule *module;
51 struct
52 {
53 struct Option *input, *datapos, *infile, *output, *samplingpos, *outfile, *method;
54 } parm;
55 int num_outputs;
56 int num_inputs;
57 struct map_store *inputs = NULL;
58 struct map_store *outputs = NULL;
59 int interpol_method = LINEAR_INTERPOLATION;
60
61 G_gisinit(argv[0]);
62
63 module = G_define_module();
64 G_add_keyword(_("raster"));
65 G_add_keyword(_("series"));
66 G_add_keyword(_("interpolation"));
67 /* TODO: re-phrase the description */
68 module->description =
69 _("Interpolates raster maps located (temporal or spatial) "
70 "in between input raster maps at specific sampling positions.");
71
72 parm.input = G_define_standard_option(G_OPT_R_INPUTS);
73 parm.input->required = NO;
74
75 parm.datapos = G_define_option();
76 parm.datapos->key = "datapos";
77 parm.datapos->type = TYPE_DOUBLE;
78 parm.datapos->required = NO;
79 parm.datapos->description = _("Data point position for each input map");
80 parm.datapos->multiple = YES;
81
82 parm.infile = G_define_standard_option(G_OPT_F_INPUT);
83 parm.infile->key = "infile";
84 parm.infile->description = _("Input file with one input raster map name and data point position per line,"
85 " field separator between name and sample point is |");
86 parm.infile->required = NO;
87
88 parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
89 parm.output->multiple = YES;
90 parm.output->required = NO;
91
92 parm.samplingpos = G_define_option();
93 parm.samplingpos->key = "samplingpos";
94 parm.samplingpos->type = TYPE_DOUBLE;
95 parm.samplingpos->required = NO;
96 parm.samplingpos->multiple = YES;
97 parm.samplingpos->description = _("Sampling point position for each output map");
98
99 parm.outfile = G_define_standard_option(G_OPT_F_INPUT);
100 parm.outfile->key = "outfile";
101 parm.outfile->description = _("Input file with one output raster map name and sample point position per line,"
102 " field separator between name and sample point is |");
103 parm.outfile->required = NO;
104
105 parm.method = G_define_option();
106 parm.method->key = "method";
107 parm.method->type = TYPE_STRING;
108 parm.method->required = NO;
109 parm.method->options = "linear";
110 parm.method->answer = "linear";
111 parm.method->description = _("Interpolation method, currently only linear interpolation is supported");
112 parm.method->multiple = NO;
113
114 if (G_parser(argc, argv))
115 exit(EXIT_FAILURE);
116
117 if (parm.output->answer && parm.outfile->answer)
118 G_fatal_error(_("%s= and %s= are mutually exclusive"),
119 parm.output->key, parm.outfile->key);
120
121 if (parm.samplingpos->answer && parm.outfile->answer)
122 G_fatal_error(_("%s= and %s= are mutually exclusive"),
123 parm.samplingpos->key, parm.outfile->key);
124
125 if (!parm.output->answer && !parm.outfile->answer)
126 G_fatal_error(_("Please specify %s= or %s="),
127 parm.output->key, parm.outfile->key);
128
129 if (parm.output->answer && !parm.samplingpos->answer)
130 G_fatal_error(_("Please specify %s= and %s="),
131 parm.output->key, parm.samplingpos->key);
132
133 if (parm.input->answer && parm.infile->answer)
134 G_fatal_error(_("%s= and %s= are mutually exclusive"),
135 parm.input->key, parm.infile->key);
136
137 if (parm.datapos->answer && parm.infile->answer)
138 G_fatal_error(_("%s= and %s= are mutually exclusive"),
139 parm.datapos->key, parm.infile->key);
140
141 if (!parm.input->answer && !parm.infile->answer)
142 G_fatal_error(_("Please specify %s= or %s="),
143 parm.input->key, parm.infile->key);
144
145 if (parm.input->answer && !parm.datapos->answer)
146 G_fatal_error(_("Please specify %s= and %s="),
147 parm.input->key, parm.datapos->key);
148
149 if(G_strncasecmp(parm.method->answer, "linear", 6) == 0)
150 interpol_method = LINEAR_INTERPOLATION;
151
152 if(G_strncasecmp(parm.method->answer, "spline", 6) == 0)
153 interpol_method = SPLINE_INTERPOLATION;
154
155 inputs = get_parameter_input("input", parm.input->answers, parm.datapos->answers, parm.infile->answer, &num_inputs);
156 outputs = get_parameter_input("output", parm.output->answers, parm.samplingpos->answers, parm.outfile->answer, &num_outputs);
157
158 start_interpolation(inputs, num_inputs, outputs, num_outputs, interpol_method);
159
160 exit(EXIT_SUCCESS);
161 }
162
163 /* *************************************************************** */
164 /* *************************************************************** */
165 /* *************************************************************** */
166
get_parameter_input(const char * type,char ** map_names,char ** positions,char * file,int * number_of_maps)167 struct map_store *get_parameter_input(const char *type, char **map_names, char **positions, char *file, int *number_of_maps)
168 {
169 struct map_store *maps = NULL;
170 int max_maps;
171 int num_maps;
172 int num_points;
173
174 /* process the output maps from the infile */
175 if (file) {
176 FILE *in;
177
178 in = fopen(file, "r");
179 if (!in)
180 G_fatal_error(_("Unable to open %s file <%s>"), type, file);
181
182 num_maps = 0;
183 max_maps = 0;
184
185 for (;;) {
186 char buf[GNAME_MAX + 50]; /* Name and position */
187 char tok_buf[GNAME_MAX + 50]; /* Name and position */
188 char *name;
189 int ntokens;
190 char **tokens;
191 struct map_store *p;
192 double pos = -1;
193
194 if (!G_getl2(buf, sizeof(buf), in))
195 break;
196
197 strcpy(tok_buf, buf);
198 tokens = G_tokenize(tok_buf, "|");
199 ntokens = G_number_of_tokens(tokens);
200
201 if(ntokens > 1) {
202 name = G_chop(tokens[0]);
203 pos = atof(G_chop(tokens[1]));
204 } else {
205 name = G_chop(buf);
206 }
207
208 /* Ignore empty lines */
209 if (!*name)
210 continue;
211
212 if(pos == -1)
213 G_fatal_error(_("Missing point position for %s map <%s>"
214 " in file <%s> near line %i"),
215 type, name, file, num_maps);
216
217 if (num_maps >= max_maps) {
218 max_maps += 100;
219 maps = G_realloc(maps, max_maps * sizeof(struct map_store));
220 }
221 p = &maps[num_maps++];
222
223 p->name = G_store(name);
224 p->pos = pos;
225 p->fd = -1;
226 p->buf = NULL;
227 p->has_run = 0;
228 G_verbose_message(_("Preparing %s map <%s> at position %g"), type, p->name, p->pos);
229 }
230
231 if (num_maps < 1)
232 G_fatal_error(_("No raster map name found in %s file <%s>"), type, file);
233
234 fclose(in);
235 }
236 else {
237 int i;
238 for (i = 0; map_names[i]; i++)
239 ;
240 num_maps = i;
241
242 if (num_maps < 1)
243 G_fatal_error(_("No %s raster map not found"), type);
244
245 for (i = 0; positions[i]; i++)
246 ;
247 num_points = i;
248
249 if (num_points != num_maps)
250 G_fatal_error(_("The number of %s maps and %s point positions must be equal"), type, type);
251
252 maps = G_malloc(num_maps * sizeof(struct map_store));
253
254 for (i = 0; i < num_maps; i++) {
255 struct map_store *p = &maps[i];
256
257 p->name = map_names[i];
258 p->pos = (DCELL)atof(positions[i]);
259 p->fd = -1;
260 p->buf = NULL;
261 p->has_run = 0;
262 G_verbose_message(_("Preparing %s map <%s> at position %g"), type, p->name, p->pos);
263 }
264 }
265 *number_of_maps = num_maps;
266 return maps;
267 }
268
269 /* *************************************************************** */
270 /* *************************************************************** */
271 /* *************************************************************** */
272
start_interpolation(struct map_store * inputs,int num_inputs,struct map_store * outputs,int num_outputs,int interpol_method)273 void start_interpolation(struct map_store *inputs, int num_inputs, struct map_store *outputs, int num_outputs, int interpol_method)
274 {
275 int i;
276 struct map_store **inp = (struct map_store**) G_malloc(num_inputs * sizeof(struct map_store*));
277 struct map_store **outp = (struct map_store**) G_malloc(num_outputs * sizeof(struct map_store*));
278
279 G_verbose_message(_("Start interpolation run with %i input maps and %i output maps"),
280 num_inputs, num_outputs);
281
282 for(i = 0; i < num_inputs; i++)
283 inp[i] = &inputs[i];
284 for(i = 0; i < num_outputs; i++)
285 outp[i] = &outputs[i];
286
287 /* Sort input and output pointer by their point position
288 * using brute force. :)
289 */
290
291 selection_sort(inp, num_inputs);
292 selection_sort(outp, num_outputs);
293
294 if(interpol_method == LINEAR_INTERPOLATION)
295 linear_interpolation(inp, num_inputs, outp, num_outputs);
296
297 for(i = 0; i < num_outputs; i++) {
298 if(outp[i]->has_run == 0)
299 G_warning(_("map <%s> at position %g was not interpolated. Check the interpolation interval."), outp[i]->name, outp[i]->pos);
300 }
301 }
302
303 /* *************************************************************** */
304 /* *************************************************************** */
305 /* *************************************************************** */
306
linear_interpolation(struct map_store ** inp,int num_inputs,struct map_store ** outp,int num_outputs)307 void linear_interpolation(struct map_store **inp, int num_inputs,
308 struct map_store **outp, int num_outputs)
309 {
310 struct map_store *left;
311 struct map_store *right;
312 struct History history;
313 int interval, l, row;
314 int nrows = Rast_window_rows();
315 int ncols = Rast_window_cols();
316 int start = 0;
317
318 if(num_inputs < 2)
319 G_fatal_error(_("At least 2 input maps are required for linear interpolation"));
320
321 /* Interpolate for each interval */
322 for(interval = 0; interval < num_inputs - 1; interval++) {
323
324 left = inp[interval];
325 right = inp[interval + 1];
326 left->fd = Rast_open_old(left->name, "");
327 right->fd = Rast_open_old(right->name, "");
328 left->buf = Rast_allocate_d_buf();
329 right->buf = Rast_allocate_d_buf();
330
331 for(l = start; l < num_outputs; l++) {
332 /* Check if the map is in the interval and process it */
333 if(outp[l]->pos >= left->pos && outp[l]->pos <= right->pos) {
334 outp[l]->fd = Rast_open_new(outp[l]->name, DCELL_TYPE);
335 outp[l]->buf = Rast_allocate_d_buf();
336
337 G_verbose_message(_("Interpolate map <%s> at position %g in interval (%g;%g)"),
338 outp[l]->name, outp[l]->pos, left->pos, right->pos);
339 G_verbose_message(_("Percent complete..."));
340
341 for (row = 0; row < nrows; row++) {
342 G_percent(row, nrows, 2);
343
344 Rast_get_d_row(left->fd, left->buf, row);
345 Rast_get_d_row(right->fd, right->buf, row);
346
347 interpolate_row_linear(left, right, outp[l], ncols);
348 Rast_put_d_row(outp[l]->fd, outp[l]->buf);
349 }
350
351 G_percent(row, nrows, 2);
352
353 Rast_close(outp[l]->fd);
354 Rast_short_history(outp[l]->name, "raster", &history);
355 Rast_command_history(&history);
356 Rast_write_history(outp[l]->name, &history);
357
358 G_free(outp[l]->buf);
359 outp[l]->has_run = 1;
360
361 start = l;
362 }
363 }
364 Rast_close(left->fd);
365 G_free(left->buf);
366 Rast_close(right->fd);
367 G_free(right->buf);
368 }
369 }
370
371 /* *************************************************************** */
372 /* *************************************************************** */
373 /* *************************************************************** */
374 /* linear function v = (1 - pos/dist) * u1 + pos/dist * u2
375 *
376 * v -> The value of the output map
377 * pos -> The normalized position of the output map
378 * u1 -> The value of the left input map
379 * u2 -> The value of the right input map
380 * dist -> The distance between the position of u1 and u2
381 *
382 * */
interpolate_row_linear(struct map_store * left,struct map_store * right,struct map_store * out,int ncols)383 void interpolate_row_linear(struct map_store *left, struct map_store *right, struct map_store *out, int ncols)
384 {
385 DCELL v;
386 DCELL u1;
387 DCELL u2;
388 DCELL dist;
389 int col;
390
391 for (col = 0; col < ncols; col++) {
392 u1 = left->buf[col];
393 u2 = right->buf[col];
394 dist = fabs(right->pos - left->pos);
395
396 if(Rast_is_d_null_value(&u1) || Rast_is_d_null_value(&u2)) {
397 Rast_set_d_null_value(&v, 1);
398 } else {
399 v = (1 - ((out->pos - left->pos)/dist)) * u1 + ((out->pos - left->pos)/dist) * u2;
400 }
401 out->buf[col] = v;
402 }
403 return;
404 }
405
406 /* *************************************************************** */
407 /* *************************************************************** */
408 /* *************************************************************** */
409
selection_sort(struct map_store ** array,int num)410 void selection_sort(struct map_store **array, int num)
411 {
412 int i, j, b;
413 struct map_store *min;
414
415 for (i = 0; i < num - 1; i++) {
416 b = i;
417 min = array[b];
418
419 for (j = i + 1; j < num; j++) {
420 if (array[j]->pos < min->pos) {
421 b = j;
422 min = array[b];
423 }
424 }
425 array[b] = array[i];
426 array[i] = min;
427 }
428 }
429
430