-
Notifications
You must be signed in to change notification settings - Fork 23
/
CalcWeights.c
executable file
·299 lines (260 loc) · 8.68 KB
/
CalcWeights.c
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
/*
* SUMMARY: CalcWeights.c - Calculate interpolation weights
* USAGE: Part of DHSVM
*
* AUTHOR: Bart Nijssen
* ORG: University of Washington, Department of Civil Engineering
* E-MAIL: nijssen@u.washington.edu
* ORIG-DATE: Apr-96
* DESCRIPTION: calculate the interpolation weights for the meteorological
* inputs for each point in the modeling area. The number of
* stations is variable.
* DESCRIP-END.
* FUNCTIONS: CalcWeights()
* COMMENTS:
* $Id: CalcWeights.c,v 1.5 2003/10/28 20:02:41 colleen Exp $
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "constants.h"
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "functions.h"
/*****************************************************************************
Function name: CalcWeights()
Purpose : Calculate interpolation weights to interpolate the
meteorological data from the various stations for every
individual pixel
Required :
METLOCATION *Station - Location of meteorological stations
int NStats - Number of meteorological stations
int NX - Number of pixels in East - West direction
int NY - Number of pixels in North - South direction
uchar ** BasinMask - BasinMask
uchar ***WeightArray - 3D array with interpolation weights
Returns : void
Modifies :
The values stored at the addresses pointed to by WeightArray (i.e. it
calculates the weights and stores them)
Comments :
*****************************************************************************/
void CalcWeights(METLOCATION * Station, int NStats, int NX, int NY,
uchar ** BasinMask, uchar **** WeightArray,
OPTIONSTRUCT * Options)
{
double *Weights; /* Array with weights for all stations */
double *Distance; /* Array with distances to all stations */
double *InvDist2; /* Array with inverse distance squared */
double Denominator; /* Sum of 1/Distance^2 */
double mindistance;
double avgdistance;
double tempdistance;
double cr, crt;
int totalweight;
int y; /* Counter for rows */
int x; /* Counter for columns */
int i, j; /* Counter for stations */
int CurrentStation; /* Station at current location (if any) */
int *stationid; /* index array for sorted list of station distances */
int *stat;
int tempid;
int closest;
int crstat;
COORD Loc; /* Location of current point */
/* Allocate memory for a 3 dimensional array */
if (DEBUG)
printf("Calculating interpolation weights for %d stations\n", NStats);
if (!((*WeightArray) = (uchar ***) calloc(NY, sizeof(uchar **))))
ReportError("CalcWeights()", 1);
for (y = 0; y < NY; y++)
if (!((*WeightArray)[y] = (uchar **) calloc(NX, sizeof(uchar *))))
ReportError("CalcWeights()", 1);
for (y = 0; y < NY; y++)
for (x = 0; x < NX; x++)
if (!((*WeightArray)[y][x] = (uchar *) calloc(NStats, sizeof(uchar))))
ReportError("CalcWeights()", 1);
/* Allocate memory for the array that will contain weights, and the array for
the distances to each of the towers, and the inverse distance squared */
if (!(Weights = (double *) calloc(NStats, sizeof(double))))
ReportError("CalcWeights()", 1);
if (!(Distance = (double *) calloc(NStats, sizeof(double))))
ReportError("CalcWeights()", 1);
if (!(InvDist2 = (double *) calloc(NStats, sizeof(double))))
ReportError("CalcWeights()", 1);
if (!(stationid = (int *) calloc(NStats, sizeof(int))))
ReportError("CalcWeights()", 1);
if (!(stat = (int *) calloc(NStats + 1, sizeof(int))))
ReportError("CalcWeights()", 1);
/* Calculate the weights for each location that is inside the basin mask */
/* note stations themselves can be outside the mask */
/* this first scheme is an inverse distance squared scheme */
if (Options->Interpolation == INVDIST) {
for (y = 0; y < NY; y++) {
Loc.N = y;
for (x = 0; x < NX; x++) {
Loc.E = x;
if (INBASIN(BasinMask[y][x])) {
if (IsStationLocation(&Loc, NStats, Station, &CurrentStation)) {
for (i = 0; i < NStats; i++) {
if (i == CurrentStation)
(*WeightArray)[y][x][i] = MAXUCHAR;
else
(*WeightArray)[y][x][i] = 0;
}
}
else {
for (i = 0, Denominator = 0; i < NStats; i++) {
Distance[i] = CalcDistance(&(Station[i].Loc), &Loc);
InvDist2[i] = 1 / (Distance[i] * Distance[i]);
Denominator += InvDist2[i];
}
for (i = 0; i < NStats; i++) {
(*WeightArray)[y][x][i] =
(uchar) Round(InvDist2[i] / Denominator * MAXUCHAR);
}
}
}
else {
for (i = 0; i < NStats; i++)
(*WeightArray)[y][x][i] = 0;
}
}
}
}
if (Options->Interpolation == NEAREST) {
/* this next scheme is a nearest station */
printf("Number of stations is %d \n", NStats);
for (y = 0; y < NY; y++) {
Loc.N = y;
for (x = 0; x < NX; x++) {
Loc.E = x;
if (INBASIN(BasinMask[y][x])) { /*we are inside the basin mask */
/* find the distance to nearest station */
mindistance = DHSVM_HUGE;
avgdistance = 0.0;
for (i = 0; i < NStats; i++) {
Distance[i] = CalcDistance(&(Station[i].Loc), &Loc);
avgdistance += Distance[i] / ((double) NStats);
if (Distance[i] < mindistance) {
mindistance = Distance[i];
closest = i;
}
}
/* got closest station */
for (i = 0; i < NStats; i++) {
if (i == closest)
(*WeightArray)[y][x][i] = MAXUCHAR;
else
(*WeightArray)[y][x][i] = 0;
}
} /* done in basin mask */
else {
for (i = 0; i < NStats; i++)
(*WeightArray)[y][x][i] = 0;
}
}
}
}
if (Options->Interpolation == VARCRESS) {
/* this next scheme is a variable radius cressman */
/* find the distance to the nearest station */
/* make a decision based on the maximum allowable radius, cr */
/* and the distance to the closest station */
/* while limiting the number of interpolation stations to three */
cr = (double) Options->CressRadius;
if (cr < 2)
ReportError("CalcWeights.c", 42);
crstat = Options->CressStations;
if (crstat < 2)
ReportError("CalcWeights.c", 42);
for (y = 0; y < NY; y++) {
Loc.N = y;
for (x = 0; x < NX; x++) {
Loc.E = x;
if (INBASIN(BasinMask[y][x])) { /*we are inside the basin mask */
/* find the distance to nearest station */
for (i = 0; i < NStats; i++) {
Distance[i] = CalcDistance(&(Station[i].Loc), &Loc);
stationid[i] = i;
}
/* got distances for each station */
/* now sort the list by distance */
for (i = 0; i < NStats; i++) {
for (j = 0; j < NStats; j++) {
if (Distance[j] > Distance[i]) {
tempdistance = Distance[i];
tempid = stationid[i];
Distance[i] = Distance[j];
stationid[i] = stationid[j];
Distance[j] = tempdistance;
stationid[j] = tempid;
}
}
}
crt = Distance[0] * 2.0;
if (crt < 1.0)
crt = 1.0;
for (i = 0, Denominator = 0; i < NStats; i++) {
if (i < crstat && Distance[i] < crt) {
InvDist2[i] =
(crt * crt - Distance[i] * Distance[i]) /
(crt * crt + Distance[i] * Distance[i]);
Denominator += InvDist2[i];
}
else
InvDist2[i] = 0.0;
}
for (i = 0; i < NStats; i++)
(*WeightArray)[y][x][stationid[i]] =
(uchar) Round(InvDist2[i] / Denominator * MAXUCHAR);
/*at this point all weights have been assigned to one or more stations */
}
}
}
}
/*check that all weights add up to MAXUCHAR */
/* and output some stats on the interpolation field */
for (i = 0; i <= NStats; i++)
stat[i] = 0;
for (i = 0; i < NStats; i++)
stationid[i] = 0;
printf("\nChecking interpolation weights\n");
printf("Sum should be 255 for all pixels \n");
printf("Some error is expected due to roundoff \n");
printf("Errors greater than +/- 2 Percent are: \n");
for (y = 0; y < NY; y++) {
for (x = 0; x < NX; x++) {
if (INBASIN(BasinMask[y][x])) { /*we are inside the basin mask */
tempid = 0;
totalweight = 0;
for (i = 0; i < NStats; i++) {
totalweight += (int) (*WeightArray)[y][x][i];
if ((*WeightArray)[y][x][i] > 0) {
tempid += 1;
stationid[i] = 1;
}
}
if (totalweight < 250 || totalweight > 260)
printf("error in interpolation weight at pixel y %d x %d : %d \n", y,
x, totalweight);
stat[tempid] += 1;
}
}
}
for (i = 0; i <= NStats; i++)
if (stat[i] > 0)
printf("%d pixels are linked to %d met stations \n", stat[i],
i);
for (i = 0; i < NStats; i++)
if (stationid[i] > 0)
printf("%s station used in interpolation \n", Station[i].Name);
/* Free memory */
free(Weights);
free(Distance);
free(InvDist2);
free(stationid);
free(stat);
}