source: grass/trunk/vector/v.surf.idw/main.c

Last change on this file was 69811, checked in by neteler, 8 years ago

Numerous typos fixed (identified with tools/fix_typos.sh)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-csrc
File size: 15.4 KB
Line 
1
2/****************************************************************************
3 *
4 * MODULE: v.surf.idw
5 * AUTHOR(S): Michael Shapiro, U.S. Army Construction Engineering Research Laboratory
6 * Improved algorithm (indexes points according to cell and ignores
7 * points outside current region) by Paul Kelly
8 * further: Radim Blazek <radim.blazek gmail.com>, Huidae Cho <grass4u gmail.com>,
9 * Glynn Clements <glynn gclements.plus.com>, Markus Neteler <neteler itc.it>
10 * OGR support by Martin Landa <landa.martin gmail.com>
11 * PURPOSE: Surface interpolation from vector point data by Inverse
12 * Distance Squared Weighting
13 * COPYRIGHT: (C) 2003-2010 by the GRASS Development Team
14 *
15 * This program is free software under the GNU General
16 * Public License (>=v2). Read the file COPYING that
17 * comes with GRASS for details.
18 *
19 *****************************************************************************/
20#include <stdlib.h>
21#include <math.h>
22#include <grass/gis.h>
23#include <grass/raster.h>
24#include <grass/glocale.h>
25#include "proto.h"
26
27int search_points = 12;
28
29long npoints = 0;
30long **npoints_currcell;
31int nsearch;
32static int i;
33
34struct Point
35{
36 double north, east;
37 double z;
38};
39struct list_Point
40{
41 double north, east;
42 double z;
43 double dist;
44};
45struct Point ***points;
46struct Point *noidxpoints = NULL;
47struct list_Point *list;
48static struct Cell_head window;
49
50int main(int argc, char *argv[])
51{
52 int fd, maskfd;
53 CELL *mask;
54 DCELL *dcell;
55 struct GModule *module;
56 struct History history;
57 int row, col;
58 int searchrow, searchcolumn, pointsfound;
59 int *shortlistrows = NULL, *shortlistcolumns = NULL;
60 long ncells = 0;
61 double north, east;
62 double dist;
63 double sum1, sum2, interp_value;
64 int n;
65 double p;
66 struct
67 {
68 struct Option *input, *npoints, *power, *output, *dfield, *col;
69 } parm;
70 struct
71 {
72 struct Flag *noindex;
73 } flag;
74 struct cell_list
75 {
76 int row, column;
77 struct cell_list *next;
78 };
79 struct cell_list **search_list = NULL, **search_list_start = NULL;
80 int max_radius, radius;
81 int searchallpoints = 0;
82 char *tmpstr1, *tmpstr2;
83
84 G_gisinit(argv[0]);
85
86 module = G_define_module();
87 G_add_keyword(_("vector"));
88 G_add_keyword(_("surface"));
89 G_add_keyword(_("interpolation"));
90 G_add_keyword(_("IDW"));
91 module->description =
92 _("Provides surface interpolation from vector point data by Inverse "
93 "Distance Squared Weighting.");
94
95 parm.input = G_define_standard_option(G_OPT_V_INPUT);
96
97 parm.dfield = G_define_standard_option(G_OPT_V_FIELD);
98
99 parm.col = G_define_standard_option(G_OPT_DB_COLUMN);
100 parm.col->required = NO;
101 parm.col->label = _("Name of attribute column with values to interpolate");
102 parm.col->description = _("If not given and input is 2D vector map then category values are used. "
103 "If input is 3D vector map then z-coordinates are used.");
104 parm.col->guisection = _("Values");
105
106 parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
107
108 parm.npoints = G_define_option();
109 parm.npoints->key = "npoints";
110 parm.npoints->key_desc = "count";
111 parm.npoints->type = TYPE_INTEGER;
112 parm.npoints->required = NO;
113 parm.npoints->description = _("Number of interpolation points");
114 parm.npoints->answer = "12";
115 parm.npoints->guisection = _("Settings");
116
117 parm.power = G_define_option();
118 parm.power->key = "power";
119 parm.power->type = TYPE_DOUBLE;
120 parm.power->answer = "2.0";
121 parm.power->label = _("Power parameter");
122 parm.power->description =
123 _("Greater values assign greater influence to closer points");
124 parm.power->guisection = _("Settings");
125
126 flag.noindex = G_define_flag();
127 flag.noindex->key = 'n';
128 flag.noindex->label = _("Don't index points by raster cell");
129 flag.noindex->description = _("Slower but uses"
130 " less memory and includes points from outside region"
131 " in the interpolation");
132 flag.noindex->guisection = _("Settings");
133
134 if (G_parser(argc, argv))
135 exit(EXIT_FAILURE);
136
137 if (sscanf(parm.npoints->answer, "%d", &search_points) != 1 ||
138 search_points < 1)
139 G_fatal_error(_("Illegal number (%s) of interpolation points"),
140 parm.npoints->answer);
141
142 list =
143 (struct list_Point *) G_calloc((size_t) search_points,
144 sizeof(struct list_Point));
145
146 p = atof(parm.power->answer);
147
148 /* get the window, dimension arrays */
149 G_get_window(&window);
150
151 if (!flag.noindex->answer) {
152 npoints_currcell = (long **)G_malloc(window.rows * sizeof(long *));
153 points =
154 (struct Point ***)G_malloc(window.rows * sizeof(struct Point **));
155
156
157 for (row = 0; row < window.rows; row++) {
158 npoints_currcell[row] =
159 (long *)G_malloc(window.cols * sizeof(long));
160 points[row] =
161 (struct Point **)G_malloc(window.cols *
162 sizeof(struct Point *));
163
164 for (col = 0; col < window.cols; col++) {
165 npoints_currcell[row][col] = 0;
166 points[row][col] = NULL;
167 }
168 }
169 }
170
171 /* read the elevation points from the input sites file */
172 read_sites(parm.input->answer, parm.dfield->answer,
173 parm.col->answer, flag.noindex->answer);
174
175 if (npoints == 0)
176 G_fatal_error(_("No points found"));
177 nsearch = npoints < search_points ? npoints : search_points;
178
179 if (!flag.noindex->answer) {
180 /* Arbitrary point to switch between searching algorithms. Could do
181 * with refinement PK */
182 if ((window.rows * window.cols) / npoints > 400) {
183 /* Using old algorithm.... */
184 searchallpoints = 1;
185 ncells = 0;
186
187 /* Make an array to contain the row and column indices that have
188 * sites in them; later will just search through all these. */
189 for (searchrow = 0; searchrow < window.rows; searchrow++)
190 for (searchcolumn = 0; searchcolumn < window.cols;
191 searchcolumn++)
192 if (npoints_currcell[searchrow][searchcolumn] > 0) {
193 shortlistrows = (int *)G_realloc(shortlistrows,
194 (1 +
195 ncells) *
196 sizeof(int));
197 shortlistcolumns =
198 (int *)G_realloc(shortlistcolumns,
199 (1 + ncells) * sizeof(int));
200 shortlistrows[ncells] = searchrow;
201 shortlistcolumns[ncells] = searchcolumn;
202 ncells++;
203 }
204 }
205 else {
206 /* Fill look-up table of row and column offsets for
207 * doing a circular region growing search looking for sites */
208 /* Use units of column width */
209 max_radius = (int)(0.5 + sqrt(window.cols * window.cols +
210 (window.rows * window.ns_res /
211 window.ew_res) * (window.rows *
212 window.ns_res /
213 window.ew_res)));
214
215 search_list =
216 (struct cell_list **)G_malloc(max_radius *
217 sizeof(struct cell_list *));
218 search_list_start =
219 (struct cell_list **)G_malloc(max_radius *
220 sizeof(struct cell_list *));
221
222 for (radius = 0; radius < max_radius; radius++)
223 search_list[radius] = NULL;
224
225 for (row = 0; row < window.rows; row++)
226 for (col = 0; col < window.cols; col++) {
227 radius = (int)sqrt(col * col +
228 (row * window.ns_res / window.ew_res) *
229 (row * window.ns_res / window.ew_res));
230 if (search_list[radius] == NULL)
231 search_list[radius] =
232 search_list_start[radius] =
233 G_malloc(sizeof(struct cell_list));
234 else
235 search_list[radius] =
236 search_list[radius]->next =
237 G_malloc(sizeof(struct cell_list));
238
239 search_list[radius]->row = row;
240 search_list[radius]->column = col;
241 search_list[radius]->next = NULL;
242 }
243 }
244 }
245
246 /* allocate buffers, etc. */
247
248 dcell = Rast_allocate_d_buf();
249
250 if ((maskfd = Rast_maskfd()) >= 0)
251 mask = Rast_allocate_c_buf();
252 else
253 mask = NULL;
254
255
256 fd = Rast_open_new(parm.output->answer, DCELL_TYPE);
257
258 /* GTC Count of window rows */
259 G_asprintf(&tmpstr1, n_("%d row", "%d rows", window.rows), window.rows);
260 /* GTC Count of window columns */
261 G_asprintf(&tmpstr2, n_("%d column", "%d columns", window.cols), window.cols);
262 /* GTC First argument is map name, second - message about number of rows, third - columns. */
263 G_important_message(_("Interpolating raster map <%s> (%s, %s)..."),
264 parm.output->answer, tmpstr1, tmpstr2);
265 G_free(tmpstr1);
266 G_free(tmpstr2);
267
268 north = window.north + window.ns_res / 2.0;
269 for (row = 0; row < window.rows; row++) {
270 G_percent(row, window.rows, 1);
271
272 if (mask)
273 Rast_get_c_row(maskfd, mask, row);
274
275 north -= window.ns_res;
276 east = window.west - window.ew_res / 2.0;
277 for (col = 0; col < window.cols; col++) {
278 east += window.ew_res;
279 /* don't interpolate outside of the mask */
280 if (mask && mask[col] == 0) {
281 Rast_set_d_null_value(&dcell[col], 1);
282 continue;
283 }
284
285 /* If current cell contains more than nsearch points just average
286 * all the points in this cell and don't look in any others */
287
288 if (!(flag.noindex->answer) && npoints_currcell[row][col] >= nsearch) {
289 sum1 = 0.0;
290 for (i = 0; i < npoints_currcell[row][col]; i++)
291 sum1 += points[row][col][i].z;
292
293 interp_value = sum1 / npoints_currcell[row][col];
294 }
295 else {
296 if (flag.noindex->answer)
297 calculate_distances_noindex(north, east);
298 else {
299 pointsfound = 0;
300 i = 0;
301
302 if (searchallpoints == 1) {
303 /* If there aren't many sites just check them all to find
304 * the nearest */
305 for (n = 0; n < ncells; n++)
306 calculate_distances(shortlistrows[n],
307 shortlistcolumns[n], north,
308 east, &pointsfound);
309 }
310 else {
311 radius = 0;
312 while (pointsfound < nsearch) {
313 /* Keep widening the search window until we find
314 * enough points */
315 search_list[radius] = search_list_start[radius];
316 while (search_list[radius] != NULL) {
317 /* Always */
318 if (row <
319 (window.rows - search_list[radius]->row)
320 && col <
321 (window.cols -
322 search_list[radius]->column)) {
323 searchrow =
324 row + search_list[radius]->row;
325 searchcolumn =
326 col + search_list[radius]->column;
327 calculate_distances(searchrow,
328 searchcolumn, north,
329 east, &pointsfound);
330 }
331
332 /* Only if at least one offset is not 0 */
333 if ((search_list[radius]->row > 0 ||
334 search_list[radius]->column > 0) &&
335 row >= search_list[radius]->row &&
336 col >= search_list[radius]->column) {
337 searchrow =
338 row - search_list[radius]->row;
339 searchcolumn =
340 col - search_list[radius]->column;
341 calculate_distances(searchrow,
342 searchcolumn, north,
343 east, &pointsfound);
344 }
345
346 /* Only if both offsets are not 0 */
347 if (search_list[radius]->row > 0 &&
348 search_list[radius]->column > 0) {
349 if (row <
350 (window.rows -
351 search_list[radius]->row) &&
352 col >= search_list[radius]->column) {
353 searchrow =
354 row + search_list[radius]->row;
355 searchcolumn =
356 col - search_list[radius]->column;
357 calculate_distances(searchrow,
358 searchcolumn,
359 north, east,
360 &pointsfound);
361 }
362 if (row >= search_list[radius]->row &&
363 col <
364 (window.cols -
365 search_list[radius]->column)) {
366 searchrow =
367 row - search_list[radius]->row;
368 searchcolumn =
369 col + search_list[radius]->column;
370 calculate_distances(searchrow,
371 searchcolumn,
372 north, east,
373 &pointsfound);
374 }
375 }
376
377 search_list[radius] =
378 search_list[radius]->next;
379 }
380 radius++;
381 }
382 }
383 }
384
385 /* interpolate */
386 sum1 = 0.0;
387 sum2 = 0.0;
388 for (n = 0; n < nsearch; n++) {
389 if ((dist = sqrt(list[n].dist))) {
390 sum1 += list[n].z / pow(dist, p);
391 sum2 += 1.0 / pow(dist, p);
392 }
393 else {
394 /* If one site is dead on the centre of the cell, ignore
395 * all the other sites and just use this value.
396 * (Unlikely when using floating point numbers?) */
397 sum1 = list[n].z;
398 sum2 = 1.0;
399 break;
400 }
401 }
402 interp_value = sum1 / sum2;
403 }
404 dcell[col] = (DCELL) interp_value;
405 }
406 Rast_put_d_row(fd, dcell);
407 }
408 G_percent(1, 1, 1);
409
410 Rast_close(fd);
411
412 /* writing history file */
413 Rast_short_history(parm.output->answer, "raster", &history);
414 Rast_command_history(&history);
415 Rast_write_history(parm.output->answer, &history);
416
417 G_done_msg(" ");
418
419 exit(EXIT_SUCCESS);
420}
421
422void newpoint(double z, double east, double north, int noindex)
423{
424 int row, column;
425
426 row = (int)((window.north - north) / window.ns_res);
427 column = (int)((east - window.west) / window.ew_res);
428
429 if (!noindex) {
430 if (row < 0 || row >= window.rows || column < 0 ||
431 column >= window.cols) ;
432 else { /* Ignore sites outside current region as can't be indexed */
433
434 points[row][column] =
435 (struct Point *)G_realloc(points[row][column],
436 (1 +
437 npoints_currcell[row][column]) *
438 sizeof(struct Point));
439 points[row][column][npoints_currcell[row][column]].north = north;
440 points[row][column][npoints_currcell[row][column]].east = east;
441 points[row][column][npoints_currcell[row][column]].z = z;
442 npoints_currcell[row][column]++;
443 npoints++;
444 }
445 }
446 else {
447 noidxpoints = (struct Point *)G_realloc(noidxpoints,
448 (1 +
449 npoints) *
450 sizeof(struct Point));
451 noidxpoints[npoints].north = north;
452 noidxpoints[npoints].east = east;
453 noidxpoints[npoints].z = z;
454 npoints++;
455 }
456}
457
458void calculate_distances(int row, int column, double north,
459 double east, int *pointsfound)
460{
461 int j, n;
462 static int max;
463 double dx, dy, dist;
464 static double maxdist;
465
466 /* Check distances and find the points to use in interpolation */
467 for (j = 0; j < npoints_currcell[row][column]; j++) {
468 /* fill list with first nsearch points */
469 if (i < nsearch) {
470 dy = points[row][column][j].north - north;
471 dx = points[row][column][j].east - east;
472 list[i].dist = dy * dy + dx * dx;
473 list[i].z = points[row][column][j].z;
474 i++;
475
476 /* find the maximum distance */
477 if (i == nsearch) {
478 maxdist = list[max = 0].dist;
479 for (n = 1; n < nsearch; n++) {
480 if (maxdist < list[n].dist)
481 maxdist = list[max = n].dist;
482 }
483 }
484 }
485 else {
486
487 /* go through rest of the points now */
488 dy = points[row][column][j].north - north;
489 dx = points[row][column][j].east - east;
490 dist = dy * dy + dx * dx;
491
492 if (dist < maxdist) {
493 /* replace the largest dist */
494 list[max].z = points[row][column][j].z;
495 list[max].dist = dist;
496 maxdist = list[max = 0].dist;
497 for (n = 1; n < nsearch; n++) {
498 if (maxdist < list[n].dist)
499 maxdist = list[max = n].dist;
500 }
501 }
502
503 }
504 }
505 *pointsfound += npoints_currcell[row][column];
506}
507
508void calculate_distances_noindex(double north, double east)
509{
510 int n, max;
511 double dx, dy, dist;
512 double maxdist;
513
514 /* fill list with first nsearch points */
515 for (i = 0; i < nsearch; i++) {
516 dy = noidxpoints[i].north - north;
517 dx = noidxpoints[i].east - east;
518 list[i].dist = dy * dy + dx * dx;
519 list[i].z = noidxpoints[i].z;
520 }
521 /* find the maximum distance */
522 maxdist = list[max = 0].dist;
523 for (n = 1; n < nsearch; n++) {
524 if (maxdist < list[n].dist)
525 maxdist = list[max = n].dist;
526 }
527 /* go through rest of the points now */
528 for (; i < npoints; i++) {
529 dy = noidxpoints[i].north - north;
530 dx = noidxpoints[i].east - east;
531 dist = dy * dy + dx * dx;
532
533 if (dist < maxdist) {
534 /* replace the largest dist */
535 list[max].z = noidxpoints[i].z;
536 list[max].dist = dist;
537 maxdist = list[max = 0].dist;
538 for (n = 1; n < nsearch; n++) {
539 if (maxdist < list[n].dist)
540 maxdist = list[max = n].dist;
541 }
542 }
543 }
544
545}
Note: See TracBrowser for help on using the repository browser.