Skip to content. | Skip to navigation

 Back to Classic
Personal tools

airs_rdr_test.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "airs_rad_typ.h"
#include "airs_rad_struct.h"
#include "airs_chan_props.h"

/*
* We're including 2 headers with arrays the size of the number of
* AIRS channels, and they define different constants for that purpose.
* Make sure they agree.
*/
#if AIRS_RAD_CHANNEL != NUM_AIRS_CHANS
#error "Values for # of AIRS channels do not agree in airs_rad.h and airs_chan_props.h"
#endif

/*
* Utility to check if a number is indistinguishable from
* the value AIRS uses to flag invalid data: -9999.0
*/
static int is_airs_bad_float(double arg) {
return (fabs(arg / -9999.0 - 1.0) < 1.0e-6);
}

/*
* Note: this struct is about 121 MBytes. On my system this works
* when allocated this way but not when allocated as an automatic
* local variable in main(). It could also be dynamically allocated
* with malloc().
*/
static airs_rad_gran_t airs_rad_gran;

int main(int argc, char * argv[]) {
int chan; /* 0-based channel index */
int nbad; /* count how many chans bad each way (reused many times) */
int track; /* 0-based index along track */
int xtrack;/* 0-based index across-track */
int chanisbad[AIRS_RAD_CHANNEL]; /* keep track of channels eliminated */
double mean = 0.0;
double dev = 0.0;
double rdiff_min = 0.0;
double rdiff_max = 0.0;
int ngood;
int nprocess=0, nspecial=0, nmissing=0, nother=0;
char * file_name = 0;
char cp_fname[256]; /* channel properties file name */
char *cp_dirname; /* channel properties directory name */
int have_cp_file = 0;
airs_chan_props_t chan_props[NUM_AIRS_CHANS];

if (argc != 2 && argc != 3) {
fprintf(stderr, "%s evaluates the quality of data in an AIRS Level-1B\n",
argv[0]);
fprintf(stderr, "file. It requires exactly one argument: the name of\n");
fprintf(stderr, "the file. An optional second argument gives the\n");
fprintf(stderr, "from which to read channel properties files.\n\n");
exit(EXIT_FAILURE);
}
file_name = argv[1];

memset(chan_props, 0, sizeof(chan_props));

printf("Scanning quality of data in AIRS Level-1B infrared radiance file:\n");
printf(" \"%s\"\n\n", file_name);

airs_rad_rdr(file_name, &airs_rad_gran);

printf("TAI Time = %f ... %f seconds past start of 1993.\n",
airs_rad_gran.Time[0][0],
airs_rad_gran.Time[AIRS_RAD_GEOTRACK-1][AIRS_RAD_GEOXTRACK-1]);
printf("This granule spans a range of %.2f seconds or %.1f minutes.\n",
airs_rad_gran.Time[AIRS_RAD_GEOTRACK-1][AIRS_RAD_GEOXTRACK-1]
- airs_rad_gran.Time[0][0],
(airs_rad_gran.Time[AIRS_RAD_GEOTRACK-1][AIRS_RAD_GEOXTRACK-1]
- airs_rad_gran.Time[0][0])/60.0);
printf("Nominal AIRS granule duration is 6.0 minutes.\n\n");

if (argc > 2) {
cp_dirname = argv[2];

/* find the correct channel properties file for this data time */
have_cp_file = !select_chan_props(cp_fname, cp_dirname,
airs_rad_gran.start_year,
airs_rad_gran.start_month,
airs_rad_gran.start_day);

/* if that worked, read the chan props file */
if (have_cp_file) {
char cp_full_path_fname[512];

sprintf(cp_full_path_fname, "%s/%s", cp_dirname, cp_fname);
have_cp_file = !read_chan_props(cp_full_path_fname, chan_props);
if (have_cp_file) {
printf("Read channel properties file:\n %s\n %s\n\n",
cp_dirname, cp_fname);
}
}
} else {
printf("Because no channel properties file directory was specified, that\n");
printf("information cannot be used in evaluating the data.\n\n");
}


printf("Checking channel attributes that are constant for granule...\n");
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
chanisbad[chan] = 0;

printf(" These channels should not be used because ExcludedChans indicates\n");
printf(" known radiometric problems. ExcludedChans echos information in\n");
printf(" field AB_State of the Channel Properties File:\n ");
nbad = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
{
if (airs_rad_gran.ExcludedChans[chan] > 2) {
chanisbad[chan] = 1;
printf("%4d ", chan + 1); /* 1-based index on output */
nbad++;
if ((nbad % 10) == 0)
printf("\n ");
}
}
if (nbad == 0) printf("(None!)\n");
printf("\n\n");

if (have_cp_file) {
nbad = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
if (chan_props[chan].AB_State > 2 && !chanisbad[chan])
nbad++;
printf(" Since we have access to the channel properties file, we'll also\n");
printf(" directly check AB_State from there. This should be identical to\n");
printf(" ExcludedChans, but it sometimes happens that we have a better\n");
printf(" channel properties file available locally than was used when\n");
printf(" the level-1B data was produced. In this case %d additional\n", nbad);
printf(" channels are flagged:\n ");
nbad = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
if (chan_props[chan].AB_State > 2 && !chanisbad[chan]) {
chanisbad[chan] = 1;
printf("%4d ", chan + 1); /* 1-based index on output */
nbad++;
if ((nbad % 10) == 0)
printf("\n ");
}
printf("\n\n");
}

printf(" These channels should not be used because CalChanSummary indicates\n");
printf(" a transient problem with high noise:\n ");
nbad = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
{
if (airs_rad_gran.CalChanSummary[chan] & 8) {
if (!chanisbad[chan]) { /* don't report a channel more than once */
printf("%4d ", chan + 1); /* 1-based index on output */
nbad++;
if ((nbad % 10) == 0)
printf("\n ");
chanisbad[chan] = 1;
}
}
}
if (nbad == 0) printf("(None!)\n");
printf("\n");
printf(" But note that the high noise flag only means that the granule's\n");
printf(" internal measurement of noise is much higher than the static\n");
printf(" value given in the Channel Properties File. Other channels\n");
printf(" may have higher noise levels, but are not flagged because they\n");
printf(" are conforming to their expected behavior.\n\n");

if (have_cp_file) {
nbad = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
if (chan_props[chan].NeDT > 0.7 && !chanisbad[chan])
nbad++;
printf(" Since we have access to the channel properties file, we'll also\n");
printf(" directly check the noise level from there. Let us say that for\n");
printf(" our application, we wish to exclude any channel with an expected\n");
printf(" noise level greater than 0.7 Kelvins. These %d otherwise good\n",
nbad);
printf(" channels exceed this level:\n ");
nbad = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
if (chan_props[chan].NeDT > 0.7 && !chanisbad[chan]) {
chanisbad[chan] = 1;
printf("%4d ", chan + 1); /* 1-based index on output */
nbad++;
if ((nbad % 10) == 0)
printf("\n ");
}
printf("\n\n");
}

printf(" These channels should not be used because CalChanSummary indicates\n");
printf(" a problem was encountered in calculating calibration coefficients:\n ");
nbad = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
{
/* 32 (bit 5) indicates gain anamoly, 64 (bit 6) for offset anamoly */
if (airs_rad_gran.CalChanSummary[chan] & (32 + 64)) {
if (!chanisbad[chan]) { /* don't report a channel more than once */
printf("%4d ", chan + 1); /* 1-based index on output */
nbad++;
if ((nbad % 10) == 0)
printf("\n ");
chanisbad[chan] = 1;
}
}
}
if (nbad == 0) printf("(None!)\n");
printf("\n\n");

printf("Each AIRS field-of-view also has a \"state\", indicating\n");
printf("whether the instrument was in science mode when the data\n");
printf("was taken and whether the data was successfully transmitted.\n");
printf("We'll only use fields-of-view with \"state\" = 0/\"PROCESS\"\n\n");

nbad = nprocess = nspecial = nmissing = nother = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
switch (airs_rad_gran.state[track][xtrack]) {
case 0: nprocess++; break;
case 1: nspecial++; break;
case 2: nbad++; break;
case 3: nmissing++; break;
default: nother++; break;
}
printf("Of %d fields-of-view in this granule:\n",
AIRS_RAD_GEOTRACK * AIRS_RAD_GEOXTRACK);
printf(" %5d fields of view are in state 0=Process\n", nprocess);
printf(" %5d fields of view are in state 1=Special Calibration sequence\n",
nspecial);
printf(" %5d fields of view are in state 2=Erroneous\n", nbad);
printf(" %5d fields of view are in state 3=Missing data\n", nmissing);
printf(" %5d fields of view have some other value for state "
"(should never happen)\n", nother);
printf("\n\n");

nbad = nother = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
if (airs_rad_gran.CalFlag[track][chan] & 16)
if (chanisbad[chan])
nbad++;
else
nother++;
printf("We also wish to exclude from our analysis any channels that have\n");
printf("experienced a sudden discontinuity in signal level, or \"pop\".\n");
printf("But for these events we only need to exclude a given channel\n");
printf("for the specific scan during which the event took place. In\n");
printf("this granule there were %d pops recorded. %d of these were in\n",
nbad + nother, nbad);
printf("channels we are already not using because of static or granule-\n");
printf("wide problems discussed above, leaving %d important \"pops\".\n",
nother);
if (nbad + nother) {
printf("The locations are:\n");
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
if (airs_rad_gran.CalFlag[track][chan] & 16)
printf(" Scan %3d channel %4d%s\n",
track+1, chan+1, /* +1 gives human-readable 1-based indices */
chanisbad[chan] ? " (known bad channel)" : "");
}
printf("\n\n");

printf("If you are sensitive to the shape of the spectral response\n");
printf("function because you are using a radiative transfer algorithm,\n");
printf("such as the AIRS RTA, that requires a well-charactersized shape,\n");
printf("then you should avoid using any channels that have the \"SRF Shape\"\n");
printf("comment in the Channel Properties File.\n\n");

if (have_cp_file) {
nbad = ngood = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
if (chan_props[chan].comment == SRF_Shape)
if (chanisbad[chan])
nbad++;
else
ngood++;
printf("There are %d channels with \"SRF Shape\" in the channel properties\n",
nbad+ngood);
printf("file but %d of them are channels that we have already eliminated from\n",
nbad);
printf("for one reason or another, leaving %d uniquely marked by "
"\"SRF Shape\"\n\n",
ngood);
}

printf("If you are sensitive to channel coregistration (Cij) then you can\n");
printf("filter for that one of two ways:\n");
printf(" - Permanently avoid the most misaligned channels, based on the \n");
printf(" boresight centroids in the channel properties file.\n");
printf(" - The field \"Cij\" in the channel properties file gives a measure\n");
printf(" of the spatial coregistration of a channel with reference channel.\n");
if (have_cp_file) {
nbad = ngood = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
/* NOTE Cij = 9.9999 when no good value available. This tests catches these. */
if (chan_props[chan].Cij < 0.96 || chan_props[chan].Cij > 1.1)
if (chanisbad[chan])
nbad++;
else
ngood++;
printf(" Requiring spatial coregistration to be 96%% or greater,\n");
printf(" eliminates %d of %d channels (%.1f%%). But %d\n",
ngood+nbad, AIRS_RAD_CHANNEL,
(ngood+(double)nbad)*100.0/AIRS_RAD_CHANNEL, nbad);
printf(" of these channels are already flagged for another reason,\n");
printf(" leaving %d uniquely flagged for Cij.\n", ngood);
}
printf(" - The fields \"Centroid X\" and \"Centroid Y\" in the channel\n");
printf(" properties file give the x- and y-offset in millidegrees from\n");
printf(" the boresight.\n");
if (have_cp_file) {
nbad = ngood = 0;
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++)
/*
* If we're equally sensitive to errors in X and Y directions and equally
* sensitive to positive and negative excursions then get overall
* distance as RSS of x & y coordinates.
*
* NOTE that centroids can be 9999.9 where not determined. These cases
* will all be rejected by this test without any additional logic.
*/
if (sqrt( chan_props[chan].centroid_x * chan_props[chan].centroid_x
+ chan_props[chan].centroid_y * chan_props[chan].centroid_y)
> 75.0)
if (chanisbad[chan])
nbad++;
else
ngood++;
printf(" Requiring the overall deviation to be less than 75 millidegrees\n");
printf(" eliminates %d of %d channels (%.1f%%). But %d\n",
ngood+nbad, AIRS_RAD_CHANNEL,
(ngood+(double)nbad)*100.0/AIRS_RAD_CHANNEL, nbad);
printf(" of these channels are already flagged for another reason,\n");
printf(" leaving %d uniquely flagged for centroid deviation.\n", ngood);
}
printf(" - Dynamically exclude scenes where the radiance difference\n");
printf(" between two channels with nearly the same frequency but\n");
printf(" boresight centroids indicates that misalignment is an issue.\n");
printf(" - SceneInhomogeneous gives one indication of such scenes.\n");
nbad = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
if ((airs_rad_gran.SceneInhomogeneous[track][xtrack] & (128 + 64)) == (128 + 64))
nbad++;
printf(" %5d of %d scenes in this granule trigger both bits\n",
nbad, AIRS_RAD_GEOTRACK * AIRS_RAD_GEOXTRACK);

nbad = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
if ((airs_rad_gran.SceneInhomogeneous[track][xtrack] & 64) == 64)
nbad++;
printf(" %5d of %d scenes in this granule trigger just the longwave "
"window bit\n", nbad, AIRS_RAD_GEOTRACK * AIRS_RAD_GEOXTRACK);

nbad = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
if ((airs_rad_gran.SceneInhomogeneous[track][xtrack] & 128) == 128)
nbad++;
printf(" %5d of %d scenes in this granule trigger just the shortwave "
"window bit\n", nbad, AIRS_RAD_GEOTRACK * AIRS_RAD_GEOXTRACK);

printf(" - Rdiff_lwindow and Rdiff_swindow give the actual radiance\n");
printf(" differences for two such pairs -- one in the longwave \n");
printf(" window spectral region and one in the shortwave window\n");
printf(" region.\n");
mean = 0.0;
dev = 0.0;
ngood = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
if ( (airs_rad_gran.state[track][xtrack] == 0)
&& !is_airs_bad_float(airs_rad_gran.Rdiff_lwindow[track][xtrack])) {
float rdiff = airs_rad_gran.Rdiff_lwindow[track][xtrack];
ngood++;
mean += rdiff;
dev += rdiff * rdiff;
if (ngood == 1) {
rdiff_max = rdiff;
rdiff_min = rdiff;
} else {
if (rdiff > rdiff_max) rdiff_max = rdiff;
if (rdiff < rdiff_min) rdiff_min = rdiff;
}
}
}
if (ngood > 0) {
if (ngood == 1)
dev = 0.0;
else {
double sqrt_arg = (dev - mean * mean / ngood) / (ngood - 1);
if (sqrt_arg > 0.0)
dev = sqrt(sqrt_arg);
else /* may just be barely negative from rounding errors */
dev = 0.0;
}
mean /= ngood;
printf(" In this granule Rdiff_lwindow ranges [%.3f, %.3f] "
"with a mean %.3f +/- %.3f\n", rdiff_min, rdiff_max, mean, dev);
} else
printf(" In this granule Rdiff_lwindow is never good.\n");
mean = 0.0;
dev = 0.0;
ngood = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
if ( (airs_rad_gran.state[track][xtrack] == 0)
&& !is_airs_bad_float(airs_rad_gran.Rdiff_swindow[track][xtrack])) {
float rdiff = airs_rad_gran.Rdiff_swindow[track][xtrack];
ngood++;
mean += rdiff;
dev += rdiff * rdiff;
if (ngood == 1) {
rdiff_max = rdiff;
rdiff_min = rdiff;
} else {
if (rdiff > rdiff_max) rdiff_max = rdiff;
if (rdiff < rdiff_min) rdiff_min = rdiff;
}
}
}
if (ngood > 0) {
if (ngood == 1)
dev = 0.0;
else {
double sqrt_arg = (dev - mean * mean / ngood) / (ngood - 1);
if (sqrt_arg > 0.0)
dev = sqrt(sqrt_arg);
else /* may just be barely negative from rounding errors */
dev = 0.0;
}
mean /= ngood;
printf(" In this granule Rdiff_swindow ranges [%.3f, %.3f] "
"with a mean %.3f +/- %.3f\n", rdiff_min, rdiff_max, mean, dev);
} else
printf(" In this granule Rdiff_swindow is never good.\n");

printf(" - Or you can craft your own check using the actual radiances of\n");
printf(" these channel pairs.\n");
printf(" Rdiff_lwindow uses channels:\n");
printf(" %4d @ %7.2f cm**-1\n",
airs_rad_gran.Rdiff_lwindow_M8_chan,
/* note "-1" below to convert from human-friendly 1-based index */
airs_rad_gran.nominal_freq[airs_rad_gran.Rdiff_lwindow_M8_chan-1]);
printf(" %4d @ %7.2f cm**-1\n",
airs_rad_gran.Rdiff_lwindow_M9_chan,
/* note "-1" below to convert from human-friendly 1-based index */
airs_rad_gran.nominal_freq[airs_rad_gran.Rdiff_lwindow_M9_chan-1]);
printf(" For the first field-of-view in the file these have\n");
printf(" radiances of %.2f and %.2f milliWatts/m**2/cm**-1/steradian\n",
airs_rad_gran.radiances[0][0][airs_rad_gran.Rdiff_lwindow_M8_chan-1],
airs_rad_gran.radiances[0][0][airs_rad_gran.Rdiff_lwindow_M9_chan-1]);
printf(" respectively, with noise levels of %.2f and %.2f respectively\n",
airs_rad_gran.NeN[airs_rad_gran.Rdiff_lwindow_M8_chan-1],
airs_rad_gran.NeN[airs_rad_gran.Rdiff_lwindow_M9_chan-1]);
printf("\n");
printf(" Rdiff_swindow uses channels:\n");
printf(" %4d @ %7.2f cm**-1\n",
airs_rad_gran.Rdiff_swindow_M1a_chan,
/* note "-1" below to convert from human-friendly 1-based index */
airs_rad_gran.nominal_freq[airs_rad_gran.Rdiff_swindow_M1a_chan-1]);
printf(" %4d @ %7.2f cm**-1\n",
airs_rad_gran.Rdiff_swindow_M2a_chan,
/* note "-1" below to convert from human-friendly 1-based index */
airs_rad_gran.nominal_freq[airs_rad_gran.Rdiff_swindow_M2a_chan-1]);
printf(" For the first field-of-view in the file these have\n");
printf(" radiances of %.4f and %.4f milliWatts/m**2/cm**-1/steradian\n",
airs_rad_gran.radiances[0][0][airs_rad_gran.Rdiff_swindow_M1a_chan-1],
airs_rad_gran.radiances[0][0][airs_rad_gran.Rdiff_swindow_M2a_chan-1]);
printf(" respectively, with noise levels of %.4f and %.4f respectively\n",
airs_rad_gran.NeN[airs_rad_gran.Rdiff_swindow_M1a_chan-1],
airs_rad_gran.NeN[airs_rad_gran.Rdiff_swindow_M2a_chan-1]);
printf("\n\n");
nbad = ngood = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++) {
if (is_airs_bad_float(airs_rad_gran.radiances[track][xtrack][chan])) {
if (chanisbad[chan]) {
/*
* Channel is bad for the whole granule. chanisbad[] wraps
* up the results of tests for high noise, bad calibration,
* and detectors marked unusable in the channel properties file.
*/
nbad++;
} else if (airs_rad_gran.state[track][xtrack] != 0) {
/* bad state for theis FOV */
nbad++;
} else if (airs_rad_gran.CalFlag[track][chan] & 16) {
/* pop on this scan */
nbad++;
} else
ngood++;
}
}
printf("When using radiances you also have to be aware that individual\n");
printf("readings will be flagged with a value of -9999.0 if no valid\n");
printf("radiance could be calculated. In this granule %d radiances "
"(%.2f%%) are\n",
nbad+ngood,
((double)nbad+ngood)
/((double)AIRS_RAD_GEOTRACK*AIRS_RAD_GEOXTRACK*AIRS_RAD_CHANNEL));
printf("-9999.0, but %d of these are in channels or fields-of-view that\n",
nbad);
printf("are eliminated for reasons discussed above, leaving %d marked only\n",
ngood);
printf("by having radiance = -9999.0\n\n");

nbad = ngood = 0;
for (track = 0; track < AIRS_RAD_GEOTRACK; track++)
for (xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
for (chan = 0; chan < AIRS_RAD_CHANNEL; chan++) {
/* look for cases where radiance is zero or negative but not -9999.0 */
if ( !is_airs_bad_float(airs_rad_gran.radiances[track][xtrack][chan])
&& (airs_rad_gran.radiances[track][xtrack][chan] <= 0.0)) {
if (chanisbad[chan]) {
/*
* Channel is bad for the whole granule. chanisbad[] wraps
* up the results of tests for high noise, bad calibration,
* and detectors marked unusable in the channel properties file.
*/
nbad++;
} else if (airs_rad_gran.state[track][xtrack] != 0) {
/* bad state for theis FOV */
nbad++;
} else if (airs_rad_gran.CalFlag[track][chan] & 16) {
/* pop on this scan */
nbad++;
} else
ngood++;
}
}
printf("In addition, remember that AIRS radiances can be zero or negative.\n");
printf("This typically happens when the signal level is less than or\n");
printf("on the order of the noise level. So it is more likely for\n");
printf("cold scenes and noisy channels. In this granule %d radiances "
"(%.2f%%) are\n",
nbad+ngood,
((double)nbad+ngood)
/((double)AIRS_RAD_GEOTRACK*AIRS_RAD_GEOXTRACK*AIRS_RAD_CHANNEL));
printf("nonpositive but not -9999.0, but %d of these are in channels or "
"fields-of-view that\n",
nbad);
printf("are eliminated for reasons discussed above, leaving %d \"real\"\n",
ngood);
printf("cases of radiance <= 0.0 but != -9999.0\n");

return 0;
}

Document Actions
NASA Logo - nasa.gov
NASA Privacy Policy and Important Notices
Last updated: Sep 09, 2009 02:25 PM ET
Top