145 #define VGROUPCLASS "PlanetaryGrid"
148 #define REGISTRATION CENTER
153 #define MAX_NORTH 89.0
154 #define MAX_SOUTH -89.0
155 #define MAX_WEST -179.0
156 #define MAX_EAST 179.0
158 #define VALIDMIN -99.1
161 int main(
int argc,
char *argv[]) {
167 float32 minval = 1200.0;
168 float32 maxval = -100.0;
170 float32 loZ_WIND = -30.0;
171 float32 hiZ_WIND = 30.0;
172 float32 loM_WIND = -30.0;
173 float32 hiM_WIND = 30.0;
174 float32 loPRESS = 850.0;
175 float32 hiPRESS = 1100.0;
176 float32 loREL_HUM = 5.0;
177 float32 hiREL_HUM = 100.0;
178 float32 loP_WATER = 0.;
179 float32 hiP_WATER = 200.;
197 int32 shape_nrt[2], shape_clm[2];
211 int32 sdsid, sdfid, fid, gridid, geomid;
213 float32 stdlimit1, stdlimit2, stdlimit3;
219 float32 *float_SDSdataRT;
220 float32 *float_SDSdataCLMAVG;
221 float32 *float_SDSdataCLMSTD;
222 float32 *float_SDSdataQC;
223 float32 *c_avg, *c_std;
229 int resize_2d(
float *,
int,
int,
int,
int,
float *);
230 int rd_size(
char *,
int *,
int *);
239 strcpy(infileRT, argv[1]);
240 strcpy(infileCLM, argv[2]);
242 strcpy(monthstr, argv[4]);
244 if (!strcmp(
lowcase(argv[5]),
"ozone"))
245 pexit(
"OZONE parameter type not supported by this program.");
247 if (!strcmp(
lowcase(argv[5]),
"z_wind")) {
250 }
else if (!strcmp(
lowcase(argv[5]),
"m_wind")) {
253 }
else if (!strcmp(
lowcase(argv[5]),
"press")) {
256 }
else if (!strcmp(
lowcase(argv[5]),
"rel_hum")) {
259 }
else if (!strcmp(
lowcase(argv[5]),
"p_water")) {
264 "parameter must be 'z_wind', 'm_wind', 'press', 'rel_hum' or 'p_water'.");
266 strcpy(vgroupname, argv[5]);
268 sscanf(argv[6],
"%f", &stdlimit1);
269 sscanf(argv[7],
"%f", &stdlimit2);
270 sscanf(argv[8],
"%f", &stdlimit3);
271 sscanf(argv[9],
"%d", &maxmissing);
272 if (argc >= 11) sscanf(argv[10],
"%f", &loval);
273 if (argc >= 12) sscanf(argv[11],
"%f", &hival);
274 if (argc >= 13) sscanf(argv[12],
"%d", &outmax);
276 strncpy(gridflag, argv[14], 1);
279 if (argc >= 15) ifileflag = atoi(argv[15]);
288 strcpy(sdsname, vgroupname);
289 strcat(sdsname,
"_mean");
291 strcpy(vgname,
"Geophysical Data");
292 strcpy(sdsname, vgroupname);
297 if (
rd_size(infileRT, (
int*) &shape_nrt[1], (
int*) &shape_nrt[0]) != 0) {
298 pexit(
"metqc: Unable to read the data field size");
308 array_size = shape_nrt[0] * shape_nrt[1];
310 if ((float_SDSdataRT =
311 (float32 *) malloc(
sizeof (float32) * array_size)) ==
NULL)
312 pexit(
"malloc float_SDSdataRT");
314 if ((float_SDSdataCLMAVG =
315 (float32 *) malloc(
sizeof (float32) * array_size)) ==
NULL)
316 pexit(
"malloc float_SDSdataCLMAVG");
318 if ((float_SDSdataCLMSTD =
319 (float32 *) malloc(
sizeof (float32) * array_size)) ==
NULL)
320 pexit(
"malloc float_SDSdataCLMSTD");
322 if ((float_SDSdataQC =
323 (float32 *) malloc(
sizeof (float32) * array_size)) ==
NULL)
324 pexit(
"malloc float_SDSdataQC");
329 if ((
rdsds(infileRT, vgname, sdsname, inShape,
330 float_SDSdataRT)) != 0)
pexit(
"rdsds RT");
332 if ((inShape[0] != shape_nrt[0]) || (inShape[1] != shape_nrt[1])) {
333 printf(
"inShape[0] %d shape_nrt[0] %d inShape[1] %d shape_nrt[1] %d\n",
334 inShape[0], shape_nrt[0], inShape[1], shape_nrt[1]);
335 pexit(
"real-time dimensions not matching expected");
346 if (
rd_size(infileCLM, (
int*) &shape_clm[1], (
int*) &shape_clm[0]) != 0) {
347 pexit(
"metqc: Unable to read the climatology data field size");
353 (float32 *) malloc(
sizeof (float32) * shape_clm[0] * shape_clm[1]))
355 pexit(
"malloc float_SDSdataCLMAVG");
358 (float32 *) malloc(
sizeof (float32) * shape_clm[0] * shape_clm[1]))
360 pexit(
"malloc float_SDSdataCLMAVG");
364 strcpy(sdsname, vgroupname);
365 strcat(sdsname,
"_mean");
367 if ((
rdsds(infileCLM, vgname, sdsname, inShape,
368 c_avg)) != 0)
pexit(
"rdsds CLMAVG");
370 if ((inShape[0] != shape_clm[0]) || (inShape[1] != shape_clm[1])) {
371 printf(
"inShape[0] %d shape_clm[0] %d inShape[1] %d shape_clm[1] %d\n",
372 inShape[0], shape_clm[0], inShape[1], shape_clm[1]);
373 pexit(
"clim avg dimensions not matching expected");
376 strcpy(sdsname, vgroupname);
377 strcat(sdsname,
"_std_dev");
379 if ((
rdsds(infileCLM, vgname, sdsname, inShape,
380 c_std)) != 0)
pexit(
"rdsds CLMSTD");
382 if ((inShape[0] != shape_clm[0]) || (inShape[1] != shape_clm[1])) {
383 printf(
"inShape[0] %d shape_clm[0] %d inShape[1] %d shape_clm[1] %d\n",
384 inShape[0], shape_clm[0], inShape[1], shape_clm[1]);
385 pexit(
"clim std dimensions not matching expected");
402 printf(
"NRT field size is - pixels: %d, lines: %d\n",
403 shape_nrt[1], shape_nrt[0]);
404 if (shape_clm[0] != shape_nrt[0] || shape_clm[1] != shape_nrt[1]) {
405 printf(
"NOTE: Climatology field size is different\n");
406 printf(
" Pixels: %d, lines: %d\n", shape_clm[1], shape_clm[0]);
407 printf(
" Resizing climatology to match NRT\n");
410 if ((shape_clm[0] <= 0) || (shape_clm[1] <= 0) ||
411 (shape_nrt[0] <= 0) || (shape_nrt[1] <= 0)) {
412 pexit(
"metqc: one of shape_clm, shape_nrt dimensions <= 0!");
414 if (
resize_2d(c_avg, shape_clm[1], shape_clm[0],
415 shape_nrt[1], shape_nrt[0], float_SDSdataCLMAVG) != 0) {
416 pexit(
"metqc: resize_2d problem for CLMAVG, exiting");
419 if (
resize_2d(c_std, shape_clm[1], shape_clm[0],
420 shape_nrt[1], shape_nrt[0], float_SDSdataCLMSTD) != 0) {
421 pexit(
"metqc: resize_2d problem for CLMSTD, exiting");
444 for (
i = 0;
i < shape_nrt[0];
i++) {
445 for (
j = 0;
j < shape_nrt[1];
j++) {
446 float_SDSdataQC[
p] = 0.0;
453 if (float_SDSdataRT[
p] <=
VALIDMIN) rtmiss++;
454 if (float_SDSdataRT[
p] < minval) minval = float_SDSdataRT[
p];
455 if (float_SDSdataRT[
p] < loval) locnt++;
457 if ((!strcmp(gridflag,
"d")) || (!strcmp(gridflag,
"D"))) {
459 (float32) float_SDSdataRT[
p] - float_SDSdataCLMAVG[
p];
461 if (float_SDSdataCLMSTD[
p] > 0)
462 float_SDSdataQC[
p] = (float32)
463 (float_SDSdataRT[
p] - float_SDSdataCLMAVG[
p]) /
464 float_SDSdataCLMSTD[
p];
468 if (float_SDSdataRT[
p] < minval) minval = float_SDSdataRT[
p];
469 if (float_SDSdataRT[
p] > maxval) maxval = float_SDSdataRT[
p];
472 if (float_SDSdataRT[
p] < loval) locnt++;
473 if (float_SDSdataRT[
p] > hival) hicnt++;
475 if (float_SDSdataQC[
p] < -1.0) neg1cnt++;
476 if (float_SDSdataQC[
p] < -2.0) neg2cnt++;
477 if (float_SDSdataQC[
p] < -3.0) neg3cnt++;
479 if (float_SDSdataQC[
p] > 1.0) pos1cnt++;
480 if (float_SDSdataQC[
p] > 2.0) pos2cnt++;
481 if (float_SDSdataQC[
p] > 3.0) pos3cnt++;
497 totalpts = array_size - missing;
499 free(float_SDSdataRT);
500 free(float_SDSdataCLMAVG);
501 free(float_SDSdataCLMSTD);
505 for (
i = 0;
i < shape_nrt[0];
i++) {
506 for (
j = 0;
j < shape_nrt[1];
j++) {
507 float_SDSdataQC[
p] = 0.0;
522 if ((
startHDF(outfile, &sdfid, &fid, DFACC_CREATE)) != 0)
523 pexit(
"Fatal error starting output HDF file");
539 pexit(
"Fatal error writing geometry");
545 strcpy(datalabel,
"QC array");
546 datatype = DFNT_FLOAT;
548 dataunit =
"std dev diff";
550 if ((sdsid =
wrtsds(sdfid,
rank, shape_nrt, datatype, datalabel,
551 float_SDSdataQC)) < 0)
pexit(
"main wrtsds");
553 free(float_SDSdataQC);
559 if ((
result =
addAttr(sdsid, dataattr, DFNT_CHAR, dataunit)) != 0)
586 if (!strcmp(gridflag,
"s")) {
588 printf(
"-------------------------------------------------------------\n");
589 printf(
"Results of comparison of real-time and climatological files:\n");
590 printf(
"%s\n%s\n", infileRT, infileCLM);
591 printf(
"Month: %s\t\tParameter: %s\n", monthstr, argv[5]);
592 printf(
"Thresholds: %8.3f %8.3f %8.3f Max Missings: %d",
593 stdlimit1, stdlimit2, stdlimit3, maxmissing);
595 printf(
"Minimum value: %8.3f Maximum: %8.3f\n", minval, maxval);
597 printf(
"Total # non-missing values: %6d\n", totalpts);
598 printf(
"Total # values <%8.3f: %d >%8.3f: %d allowed: %d",
599 loval, locnt, hival, hicnt, outmax);
606 if (locnt >= outmax || hicnt >= outmax) {
614 printf(
"Total # points/percentage < -1 STD: %6d (%5.2f percent)\n",
615 neg1cnt, (
float) neg1cnt / totalpts * 100.0);
616 printf(
"Total # points/percentage +/- 1 STD: %6d (%5.2f percent)",
617 totalpts - (neg1cnt + pos1cnt),
618 (
float) (totalpts - (pos1cnt + neg1cnt)) / totalpts * 100.0);
620 if (((
float) (totalpts - (pos1cnt + neg1cnt)) / totalpts * 100.0) <
628 printf(
"Total # points/percentage > 1 STD: %6d (%5.2f percent)\n",
629 pos1cnt, (
float) pos1cnt / totalpts * 100.0);
632 printf(
"Total # points/percentage < -2 STD: %6d (%5.2f percent)\n",
633 neg2cnt, (
float) neg2cnt / totalpts * 100.0);
634 printf(
"Total # points/percentage +/- 2 STD: %6d (%5.2f percent)",
635 totalpts - (neg2cnt + pos2cnt),
636 (
float) (totalpts - (pos2cnt + neg2cnt)) / totalpts * 100.0);
638 if (((
float) (totalpts - (pos2cnt + neg2cnt)) / totalpts * 100.0) <
646 printf(
"Total # points/percentage > 2 STD: %6d (%5.2f percent)\n",
647 pos2cnt, (
float) pos2cnt / totalpts * 100.0);
650 printf(
"Total # points/percentage < -3 STD: %6d (%5.2f percent)\n",
651 neg3cnt, (
float) neg3cnt / totalpts * 100.0);
652 printf(
"Total # points/percentage +/- 3 STD: %6d (%5.2f percent)",
653 totalpts - (neg3cnt + pos3cnt),
654 (
float) (totalpts - (pos3cnt + neg3cnt)) / totalpts * 100.0);
656 if (((
float) (totalpts - (pos3cnt + neg3cnt)) / totalpts * 100.0) <
664 printf(
"Total # points/percentage > 3 STD: %6d (%5.2f percent)\n",
665 pos3cnt, (
float) pos3cnt / totalpts * 100.0);
668 printf(
"Total # missing values: %6d\n", missing);
669 printf(
"Total # missing real-time: %6d", rtmiss);
670 if (rtmiss > maxmissing) {
676 printf(
"-------------------------------------------------------------\n");
684 printf(
"Threshold Status: FAILED\n");
687 printf(
"Threshold Status: SUCCESS\n");
700 printf(
"\n\nUsage:\n");
701 printf(
"\t%s <file><file><outfile><monthstr><param><t1><t2><t2><maxmiss><lo><hi><outmax>[diff][type] \n", argv[0]);
702 printf(
"\nWhere:\n");
703 printf(
"\tfile: Real-time file to process\n");
704 printf(
"\tfile: Climatology file\n");
705 printf(
"\toutfile: Output QC file\n");
706 printf(
"\tmonthstr: Month in string form (i.e., January)\n");
707 printf(
"\tparam: Parameter to run: z_wind m_wind press rel_hum p_water\n");
708 printf(
"\tt1: Threshold percent of pts w/in 1 STD DEV of clim\n");
709 printf(
"\tt2: Threshold percent of pts w/in 2 STD DEV of clim\n");
710 printf(
"\tt3: Threshold percent of pts w/in 3 STD DEV of clim\n");
711 printf(
"\tmaxmiss: Max num of missing NRT points permitted\n");
712 printf(
"\tlo: Lowest allowed value in NRT file\n");
713 printf(
"\thi: Highest allowed value in NRT file\n");
714 printf(
"\tlo: Maximum # of NRT points allowed outside of lo/hi before returning an error flag\n");
715 printf(
"\t [optional] parameters follow. Preceeding args \n");
716 printf(
"\t are required for subsequent ones used\n");
717 printf(
"\tdiff: [optional] Enter 'd' to see a simple difference output\n");
718 printf(
"\t grid rather than the default STD variance 's'\n");
719 printf(
"\t Difference calculations do not product reports\n");
720 printf(
"\ttype [optional] Enter 1 if NRT file is actually a CLM file\n");
721 printf(
"\t This will allow the program to read the CLM\n");
722 printf(
"\t as of test of CLM to itself\n");