File: grinterp.cc

package info (click to toggle)
gri 2.4.2-1
  • links: PTS
  • area: main
  • in suites: potato
  • size: 4,540 kB
  • ctags: 1,966
  • sloc: cpp: 32,542; lisp: 3,243; perl: 806; makefile: 548; sh: 253
file content (180 lines) | stat: -rw-r--r-- 4,886 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
#include	<string>
#include	<math.h>
#include	<stdio.h>
#include        "private.hh"
#include        "gr.hh"
extern double   _grMissingValue;/* defined in gr.c */
int             number_good_xyz(double x[], double y[], double f[], int n);
#define root_of_2  1.414213562

int
number_good_xyz(double *x, double *y, double *f, int n)
{
    int             i, good = 0;
    for (i = 0; i < n; i++)
	if (!gr_missingx(x[i]) && !gr_missingy(y[i]) && !gr_missing(f[i]))
	    good++;
    return good;
}

/*
 * gr_grid1: objective analysis of 2d field
 * 
 * DESCRIPTION: Given f[i] defined at points (x[i],y[i]), where 0<=i<n, this
 * routine uses objective analysis to get a predicted value fOut. The points
 * f[i] used to calculate fOut come from a specified search rectangle xr wide
 * by yr high.
 *
 * Also does boxcar gridding if 'method' is 0.
 * 
 * Third method exists to .... better clean this doc up!
 *
 * If there are `neighbors' or more points in a rectangular neighborhood 2*xr
 * wide and 2*yr high, then only the data inside that region are used, with x
 * and y distances weighted by xr and yr. If there are fewer than `neighbors'
 * points in the neighborhood, then the search rectangle is enlarged in area
 * [xr and yr are multiplied by sqrt(2)] and the process starts over. This
 * enlargement may be repeated up to 'enlargements' times, until either
 * `neighbors' data are located in the search region, or all the data are
 * used.  If `enlargements' is set to a negative value, the enlargement
 * process will be continued until all the data are used, if necessary.
 * 
 * When all relevant points are located, the output is calculated using a
 * weighting formula for the data within the rectangle:
 *
 * fOut = sum(f[i]*w[i])/sum(w[i]
 *
 * Various weightings 'w' are provided, each as a function of
 * the scaled distance 'r = sqrt([(x1-x2)/xr]^2 + [(y1-y2)/yr]^2)'.
 *
 * method = 0: w = 1 (boxcar smoothing)
 *
 * method = 1: w = (2 - r^2) / (2 + r^2) (Levy + Brown method "1")
 *
 * method = 2: w = [(2 - r^2) / (2 + r^2)]^2 (Levy + Brown method "2")
 *
 * method = 3: w = [(2 - r^2) / (2 + r^2)]^4 (Levy + Brown method "3")
 *
 * REFERENCE: Levy & Brown [1986 Journal of
 * Geophysical Research 91C4: 5153-5158].  Note that the present method 
 * extends for unequal weightings in the (x,y) direction, so that nonisotropic
 * data can be handled.
 * 
 * RETURN VALUE: The number of points finally used in the sum, or 0 if there was
 * an error (if n<=0).
 */
int
gr_grid1(double *x, double *y, double *f,
	 unsigned int n,
	 double x0, double y0,
	 double xr, double yr,
	 int method,
	 unsigned int neighbors,
	 int enlargements,
	 double *fOut)
{
    double          dx, dy;
    int             enlarge = 0;
    unsigned int    i;
    unsigned int    nn, num_in_rect;
    double          d2;	/* squared distance */
    double          sumw;	/* sum of weights */
    double          sumfw;	/* sum of weighted values */
    double          w;		/* weight of f[i] */
    /* Check for obvious errors */
    if (n <= 0 || neighbors == 0) {
	*fOut = _grMissingValue;
	return 0;
    }
    nn = number_good_xyz(x, y, f, n);
    if (neighbors > nn)
	neighbors = nn;
    /*
     * Search the rectangle, increasing its size if necessary.
     */
    do {
	num_in_rect = 0;
	sumw = sumfw = 0.0;
	switch (method) {
	case 0:
	    for (i = 0; i < n; i++) {
		dx = x0 - x[i];
		if (-xr <= dx && dx <= xr) {
		    dy = y0 - y[i];
		    if (-yr <= dy && dy <= yr) {
			sumw += 1.0;
			sumfw += f[i];
			num_in_rect++;
		    }
		}
	    }
	    break;
	case 1:
	    for (i = 0; i < n; i++) {
		dx = GRI_ABS(x0 - x[i]);
		if (dx <= xr) {
		    dy = GRI_ABS(y0 - y[i]);
		    if (dy <= yr) {
			dx /= xr;
			dy /= yr;
			d2 = dx * dx + dy * dy;	/* note 0<=d2<=2 */
			w = (2.0 - d2) / (2.0 + d2);
			sumw += w;
			sumfw += f[i] * w;
			num_in_rect++;
		    }
		}
	    }
	    break;
	case 2:
	    for (i = 0; i < n; i++) {
		dx = GRI_ABS(x0 - x[i]);
		if (dx <= xr) {
		    dy = GRI_ABS(y0 - y[i]);
		    if (dy <= yr) {
			dx /= xr;
			dy /= yr;
			d2 = dx * dx + dy * dy;	/* note 0<=d2<=2 */
			w = (2.0 - d2) / (2.0 + d2);
			w *= w;
			sumw += w;
			sumfw += f[i] * w;
			num_in_rect++;
		    }
		}
	    }
	    break;
	case 3:
	    for (i = 0; i < n; i++) {
		dx = GRI_ABS(x0 - x[i]);
		if (dx <= xr) {
		    dy = GRI_ABS(y0 - y[i]);
		    if (dy <= yr) {
			dx /= xr;
			dy /= yr;
			d2 = dx * dx + dy * dy;	/* note 0<=d2<=2 */
			w = (2.0 - d2) / (2.0 + d2);
			w *= w;
			w *= w;
			sumw += w;
			sumfw += f[i] * w;
			num_in_rect++;
		    }
		}
	    }
	    break;
	default:
	    return 0;		/* unknown method! */
	}
	xr *= root_of_2;
	yr *= root_of_2;
    } while ((++enlarge <= enlargements || enlargements < 0)
	     && num_in_rect < neighbors);
    if (num_in_rect > 0) {
	*fOut = sumfw / sumw;
	return num_in_rect;
    } else {
	*fOut = _grMissingValue;
	return 0;
    }
}