/*-------------------------------------------------------------------- * $Id: nearneighbor.c,v 1.5.4.2 2002/02/27 17:41:10 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; version 2 of the License. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * Based on a specified grid size, nearneighbor reads an xyz file and * determines the nearest points to each node in sectors. The default * looks for the nearest point for each quadrant. The points must also * be within a maximum search-radius from the node. For the nodes that * have a full set of nearest neighbors, a weighted average value is * computed. New feature is full support for boundary conditions so * that geographic or periodic conditions are explicitly dealt with * in the sense that a data point may wrap around to serve as a * constraint on the other side of the periodic boundary. * * Author: Paul Wessel * Date: 14-JUL-2000 * Version: 3.4.1 */ #include "gmt.h" struct NODE { /* Structure with point id and distance pairs for all sectors */ float *distance; /* Distance of nearest datapoint to this node per sector */ int *datum; /* Point id of this data point */ }; struct POINT { /* Structure with input data constraints */ float x, y, z, w; }; struct NODE *add_new_node(int n); void assign_node (struct NODE **node, int n_sector, int sector, double distance, int id); main (int argc, char **argv) { int i, j, k, ij, i0, j0, *di, dj, n_sectors = 4, sector, n, n_alloc = 5 * GMT_CHUNK, n_fields, nx_2; int n_set, n_almost, n_none, n_files = 0, n_args, fno, one_or_zero, n_expected_fields, pad[4], distance_flag = 0; int max_di, actual_max_di, ii, jj, n_req, x_wrap, y_wrap; BOOLEAN go, error = FALSE, done = FALSE, nofile = TRUE; BOOLEAN set_empty = FALSE, weighted = FALSE, wrap_180, replicate_x, replicate_y; double radius = 0.0, weight, weight_sum, *x0, *y0, dx, dy, delta, distance, factor; double *in, *shrink, km_pr_deg, x_left, x_right, y_top, y_bottom, offset, xinc2, yinc2, idx, idy; double half_y_width, y_width, half_x_width, x_width, three_over_radius; float empty = 0.0, *grd; char *outfile = CNULL, line[BUFSIZ]; FILE *fp = NULL; struct GRD_HEADER header; struct GMT_EDGEINFO edgeinfo; struct NODE **grid_node; struct POINT *point; argc = GMT_begin (argc, argv); GMT_boundcond_init (&edgeinfo); GMT_grd_init (&header, argc, argv, FALSE); pad[0] = pad[1] = pad[2] = pad[3] = 0; for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch (argv[i][1]) { /* Common parameters */ case 'H': case 'R': case 'V': case ':': case '\0': error += GMT_get_common_args (argv[i], &header.x_min, &header.x_max, &header.y_min, &header.y_max); break; /* Supplemental parameters */ case 'b': /* Input triplets [quadruplets] are binary, not ascii */ error += GMT_io_selection (&argv[i][2]); break; case 'I': GMT_getinc (&argv[i][2], &header.x_inc, &header.y_inc); break; case 'E': if (!argv[i][2]) { fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: Must specify value or NaN\n", GMT_program); error++; } else empty = (argv[i][2] == 'N' || argv[i][2] == 'n') ? GMT_f_NaN : (float)atof (&argv[i][2]); set_empty = TRUE; break; case 'F': header.node_offset = TRUE; break; case 'G': outfile = &argv[i][2]; break; case 'L': if (argv[i][2]) { error += GMT_boundcond_parse (&edgeinfo, &argv[i][2]); if (edgeinfo.gn) GMT_geographic_in = TRUE; } else GMT_geographic_in = TRUE; break; case 'N': n_sectors = atoi (&argv[i][2]); break; case 'S': GMT_getinc (&argv[i][2], &radius, &radius); if (argv[i][strlen(argv[i])-1] == 'k') distance_flag = 1; if (argv[i][strlen(argv[i])-1] == 'K') distance_flag = 2; break; case 'W': weighted = TRUE; break; default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else n_files++; } if (argc == 1 || GMT_quick) { fprintf (stderr, "nearneighbor %s - A \"Nearest neighbor\" gridding algorithm\n\n", GMT_VERSION); fprintf(stderr, "usage: nearneighbor [xyzfile(s)] -G -I[m|c][/[m|c]]\n"); fprintf(stderr, "\t-N -R -S[m|c|k|K] [-E] [-F]\n"); fprintf(stderr, "\t[-H ] [-L[]] [-V ] [-W] [-:] [-bi[s][]]\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf(stderr, "\t-G name of output grid.\n"); fprintf(stderr, "\t-I sets the grid spacing for the grid. Append m for minutes, c for seconds.\n"); fprintf(stderr, "\t-N sets number of sectors. Default is quadrant search [4].\n"); GMT_explain_option ('R'); fprintf(stderr, "\t-S sets search radius in -R, -I units; append m or c for minutes or seconds.\n"); fprintf(stderr, "\t Append k for km (implies -R,-I in degrees), use flat Earth approximation.\n"); fprintf(stderr, "\t Append K for km (implies -R,-I in degrees), use great circle distances.\n"); fprintf(stderr, "\n\tOPTIONS:\n"); fprintf(stderr, "\t-E value to use for empty nodes [Default is NaN].\n"); fprintf(stderr, "\t-F Force pixel registration [Default is gridline registration].\n"); GMT_explain_option ('H'); fprintf(stderr, "\t-L sets boundary conditions. can be either\n"); fprintf(stderr, "\t g for geographic boundary conditions, or one or both of\n"); fprintf(stderr, "\t x for periodic boundary conditions on x\n"); fprintf(stderr, "\t y for periodic boundary conditions on y\n"); fprintf (stderr, "\t If no flag is given, x is assumed to be longitude and periodic in 360.\n"); GMT_explain_option ('V'); fprintf(stderr, "\t-W input file has observation weights in 4th column.\n"); GMT_explain_option (':'); GMT_explain_option ('i'); GMT_explain_option ('n'); fprintf(stderr, "\t Default is 3 (or 4 if -W is set) columns\n"); GMT_explain_option ('.'); exit (EXIT_FAILURE); } if (!project_info.region_supplied) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify -R option\n", GMT_program); error++; } if (!outfile) { fprintf (stderr, "%s: GMT SYNTAX ERROR option -G: Must specify output file\n", GMT_program); error++; } if (n_sectors <= 0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Must specify a positive number of sectors\n", GMT_program); error++; } if (radius <= 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: Must specify a positive search radius\n", GMT_program); error++; } if (header.x_inc <= 0.0 || header.y_inc <= 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -I option. Must specify positive increment(s)\n", GMT_program); error++; } if (GMT_io.binary[0] && gmtdefs.io_header) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have header -H\n", GMT_program); error++; } n_req = (weighted) ? 4 : 3; if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = n_req; if (GMT_io.binary[0] && n_req > GMT_io.ncol[0]) { fprintf (stderr, "%s: GMT SYNTAX ERROR. binary input data must have at least %d columns\n", GMT_program, n_req); error++; } if (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ if (GMT_io.binary[0] && gmtdefs.verbose) { char *type[2] = {"double", "single"}; fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]); } GMT_grd_RI_verify (&header, 1); n_expected_fields = (GMT_io.binary[0]) ? GMT_io.ncol[0] : 3 + weighted; if (n_files > 0) nofile = FALSE; else n_files = 1; n_args = (argc > 1) ? argc : 2; if (header.node_offset) { one_or_zero = 0; offset = 0.0; xinc2 = 0.5 * header.x_inc; yinc2 = 0.5 * header.y_inc; } else { one_or_zero = 1; offset = 0.5; xinc2 = yinc2 = 0.0; } idx = 1.0 / header.x_inc; idy = 1.0 / header.y_inc; header.nx = irint ( (header.x_max - header.x_min) * idx) + one_or_zero; header.ny = irint ( (header.y_max - header.y_min) * idy) + one_or_zero; GMT_boundcond_param_prep (&header, &edgeinfo); if (gmtdefs.verbose) fprintf (stderr, "%s: Grid dimensions are nx = %d, ny = %d\n", GMT_program, header.nx, header.ny); grid_node = (struct NODE **) GMT_memory (VNULL, (size_t)(header.nx * header.ny), sizeof (struct NODE *), GMT_program); point = (struct POINT *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (struct POINT), GMT_program); di = (int *) GMT_memory (VNULL, (size_t)header.ny, sizeof (int), GMT_program); shrink = (double *) GMT_memory (VNULL, (size_t)header.ny, sizeof (double), GMT_program); x0 = (double *) GMT_memory (VNULL, (size_t)header.nx, sizeof (double), GMT_program); y0 = (double *) GMT_memory (VNULL, (size_t)header.ny, sizeof (double), GMT_program); for (i = 0; i < header.nx; i++) x0[i] = header.x_min + i * header.x_inc + xinc2; for (j = 0; j < header.ny; j++) y0[j] = header.y_max - j * header.y_inc - yinc2; if (distance_flag) { /* Input data is geographical */ km_pr_deg = 0.001 * 2.0 * M_PI * gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius / 360.0; max_di = (int) (ceil (header.nx / 2.0) + 0.1); actual_max_di = 0; for (j = 0; j < header.ny; j++) { shrink[j] = cosd (y0[j]); di[j] = (fabs (y0[j]) == 90.0) ? max_di : (int)(ceil (radius / (km_pr_deg * header.x_inc * shrink[j])) + 0.1); if (di[j] > max_di) di[j] = max_di; if (di[j] > actual_max_di) actual_max_di = di[j]; } dj = (int) (ceil (radius / (km_pr_deg * header.y_inc)) + 0.1); } else { /* Plain Cartesian data */ max_di = (int) (ceil (radius * idx) + 0.1); for (j = 0; j < header.ny; j++) di[j] = max_di; dj = (int) (ceil (radius * idy) + 0.1); actual_max_di = max_di; } factor = n_sectors / (2.0 * M_PI); x_left = header.x_min - actual_max_di * header.x_inc; x_right = header.x_max + actual_max_di * header.x_inc; y_top = header.y_max + dj * header.y_inc; y_bottom = header.y_min - dj * header.y_inc; x_width = header.x_max - header.x_min; y_width = header.y_max - header.y_min; half_x_width = 0.5 * x_width; half_y_width = 0.5 * y_width; nx_2 = edgeinfo.nxp / 2; n = 0; replicate_x = (edgeinfo.nxp && !header.node_offset); /* Gridline registration has duplicate column */ replicate_y = (edgeinfo.nyp && !header.node_offset); /* Gridline registration has duplicate row */ x_wrap = header.nx - 1; /* Add to node index to go to right column */ y_wrap = (header.ny - 1) * header.nx; /* Add to node index to go to bottom row */ for (fno = 1; !done && fno < n_args; fno++) { /* Loop over input files, if any */ if (!nofile && argv[fno][0] == '-') continue; if (nofile) { /* Just read standard input */ fp = GMT_stdin; done = TRUE; #ifdef SET_IO_MODE GMT_setmode (0); #endif } else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) { fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]); continue; } if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]); if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp); while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) { /* Not yet EOF */ if (GMT_io.status & GMT_IO_MISMATCH) { fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, n); exit (EXIT_FAILURE); } if (in[0] < x_left || in[0] > x_right) continue; if (in[1] < y_bottom || in[1] > y_top) continue; point[n].x = (float)in[0]; point[n].y = (float)in[1]; point[n].z = (float)in[2]; if (weighted) point[n].w = (float)in[3]; /* Find indeces of the node closest to this data point */ i0 = (int)floor (((in[0] - header.x_min) * idx) + offset); j0 = (int)floor (((header.y_max - in[1]) * idy) + offset); /* Loop over all nodes within radius of this node */ for (j = j0 - dj; j <= (j0 + dj); j++) { jj = j; if (GMT_y_out_of_bounds (&jj, &header, &edgeinfo, &wrap_180)) continue; /* Outside y-range */ for (i = i0 - di[jj]; i <= (i0 + di[jj]); i++) { ii = i; if (GMT_x_out_of_bounds (&ii, &header, &edgeinfo, wrap_180)) continue; /* Outside x-range */ /* Here, (ii,jj) is index of a node (k) inside the grid */ k = jj * header.nx + ii; dx = in[0] - x0[ii]; dy = in[1] - y0[jj]; /* Check for wrap-around in x or y. This should only occur if the search radius is larger than 1/2 the grid width/height so that the shortest distance is going through the periodic boundary. For longitudes the dx obviously cannot exceed 180 (half_x_width) since we could then go the other direction instead. */ if (edgeinfo.nxp && fabs (dx) > half_x_width) dx -= copysign (x_width, dx); if (edgeinfo.nyp && fabs (dy) > half_y_width) dy -= copysign (y_width, dy); switch (distance_flag) { /* Take different action depending on how we want distances calculated */ case 0: /* Cartesian distance */ distance = hypot (dx, dy); break; case 1: /* Flat Earth Approximation */ distance = km_pr_deg * hypot (dx * shrink[jj], dy); break; case 2: /* Full spherical calculation */ distance = km_pr_deg * GMT_great_circle_dist (x0[ii], y0[jj], in[0], in[1]); break; default: break; } if (distance > radius) continue; /* Data constraint is too far from this node */ /* OK, this point should constrain this node. Calculate which sector and assign the value */ sector = ((int)((d_atan2 (dy, dx) + M_PI) * factor)) % n_sectors; assign_node (&grid_node[k], n_sectors, sector, distance, n); /* With periodic, gridline-registrered grids there are duplicate rows and/or columns so we may have to assign the point to more than one node. The next section deals with this situation. */ if (replicate_x) { /* Must check if we have to replicate a column */ if (ii == 0) /* Must replicate left to right column */ assign_node (&grid_node[k+x_wrap], n_sectors, sector, distance, n); else if (ii == edgeinfo.nxp) /* Must replicate right to left column */ assign_node (&grid_node[k-x_wrap], n_sectors, sector, distance, n); } if (replicate_y) { /* Must check if we have to replicate a row */ if (jj == 0) /* Must replicate top to bottom row */ assign_node (&grid_node[k+y_wrap], n_sectors, sector, distance, n); else if (jj == edgeinfo.nyp) /* Must replicate bottom to top row */ assign_node (&grid_node[k-y_wrap], n_sectors, sector, distance, n); } } } n++; if (n == n_alloc) { n_alloc += GMT_CHUNK; point = (struct POINT *) GMT_memory ((void *)point, (size_t)n_alloc, sizeof (struct POINT), GMT_program); } } if (fp != GMT_stdin) GMT_fclose (fp); } point = (struct POINT *) GMT_memory ((void *)point, (size_t)n, sizeof (struct POINT), GMT_program); grd = (float *) GMT_memory (VNULL, (size_t)(header.nx * header.ny), sizeof (float), GMT_program); /* Compute weighted averages based on the nearest neighbors */ n_set = n_almost = n_none = 0; if (!set_empty) empty = GMT_f_NaN; three_over_radius = 3.0 / radius; for (j = ij = 0; j < header.ny; j++) { for (i = 0; i < header.nx; i++, ij++) { if (!grid_node[ij]) { /* No nearest neighbors, set to empty and goto next node */ n_none++; grd[ij] = empty; continue; } for (k = 0, go = TRUE; go && k < n_sectors; k++) if (grid_node[ij]->datum[k] < 0) go = FALSE; if (!go) { /* Not full set of neighbors in all sectors, set to empty and goto next node */ n_almost++; grd[ij] = empty; continue; } /* OK, here we have enough data and need to calculate the weighted value */ n_set++; weight_sum = grd[ij] = 0.0; /* Initialize sums */ for (k = 0; k < n_sectors; k++) { delta = three_over_radius * grid_node[ij]->distance[k]; weight = 1.0 / (1.0 + delta * delta); /* This is distance weight */ if (weighted) weight *= point[grid_node[ij]->datum[k]].w; /* This is observation weight */ grd[ij] += (float)(weight * point[grid_node[ij]->datum[k]].z); weight_sum += weight; } grd[ij] /= (float)weight_sum; } } if (GMT_write_grd (outfile, &header, grd, 0.0, 0.0, 0.0, 0.0, pad, FALSE)) { fprintf (stderr, "%s: Error writing file %s\n", GMT_program, outfile); exit (EXIT_FAILURE); } if (gmtdefs.verbose) { sprintf (line, "%s)\n", gmtdefs.d_format); fprintf (stderr, "%s: %d nodes were assigned an average value\n", GMT_program, n_set); fprintf (stderr, "%s: %d nodes failed sector criteria and %d nodes had no neighbor points (all set to ", GMT_program, n_almost, n_none); (GMT_is_dnan (empty)) ? fprintf (stderr, "NaN)\n") : fprintf (stderr, line, empty); } GMT_free ((void *)grd); GMT_free ((void *)point); GMT_free ((void *)grid_node); GMT_free ((void *)shrink); GMT_free ((void *)di); GMT_free ((void *)x0); GMT_free ((void *)y0); GMT_end (argc, argv); } struct NODE *add_new_node(int n) { struct NODE *new; new = (struct NODE *) GMT_memory (VNULL, (size_t)1, sizeof (struct NODE), GMT_program); new->distance = (float *) GMT_memory (VNULL, (size_t)n, sizeof (float), GMT_program); new->datum = (int *) GMT_memory (VNULL, (size_t)n, sizeof (int), GMT_program); while (n > 0) new->datum[--n] = -1; return (new); } void assign_node (struct NODE **node, int n_sector, int sector, double distance, int id) { /* Allocates node space if not already used and updates the value if closer to node */ if (!(*node)) *node = add_new_node (n_sector); if ((*node)->datum[sector] == -1 || (*node)->distance[sector] > distance) { (*node)->distance[sector] = (float)distance; (*node)->datum[sector] = id; } }