Skip to content

Commit

Permalink
#2867 Read the ADP QC flag values and meanings attributes from DQF va…
Browse files Browse the repository at this point in the history 
…riable and set the QC high, meduium, low values to support Enterprise algorithm. Adjusted the ADP QC values by using AOD qc values
  • Loading branch information
Howard Soh committed Apr 30, 2024
1 parent 96fde58 commit e7e25a8
Showing 1 changed file with 106 additions and 27 deletions.
133 changes: 106 additions & 27 deletions src/tools/other/point2grid/point2grid.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,9 @@ static NcFile *nc_out = (NcFile *) nullptr;
static NcDim lat_dim ;
static NcDim lon_dim ;

static bool is_baseline_algorithm = true; /* false for Enterprise algorithm */
static int adp_qc_high; /* 3 as baseline algorithm, 0 for enterpirse algorithm */
static int adp_qc_medium; /* 1 as baseline algorithm, 1 for enterpirse algorithm */
static int adp_qc_low; /* 0 as baseline algorithm, 2 for enterpirse algorithm */

////////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -173,6 +175,7 @@ static void regrid_nc_variable(NcFile *nc_in, Met2dDataFile *fr_mtddf,
static bool keep_message_type(const int mt_index);

static bool has_lat_lon_vars(NcFile *nc_in);
static void set_adp_gc_values(NcVar var_adp_qc);

////////////////////////////////////////////////////////////////////////
// for GOES 16
Expand Down Expand Up @@ -1869,17 +1872,12 @@ void check_lat_lon(int data_size, float *latitudes, float *longitudes) {

int compute_adp_qc_flag(int adp_qc, int shift_bits) {
int particle_qc = ((adp_qc >> shift_bits) & 0x03);
int qc_for_flag = particle_qc;

if (is_baseline_algorithm) {
switch (particle_qc) {
case 3: qc_for_flag = 0; break; /* high */
case 1: qc_for_flag = 1; break; /* medium */
case 0: qc_for_flag = 2; break; /* low */
default: qc_for_flag = bad_data_int; break;
}
}
else if (3 == particle_qc) qc_for_flag = bad_data_int; /* Enterprise algorithm */
int qc_for_flag;

if (particle_qc == adp_qc_high) qc_for_flag = 0;
else if (particle_qc == adp_qc_medium) qc_for_flag = 1;
else if (particle_qc == adp_qc_low) qc_for_flag = 2;
else qc_for_flag = bad_data_int;

return qc_for_flag;
}
Expand Down Expand Up @@ -2364,12 +2362,17 @@ int get_lon_count(NcFile *_nc) {

static NcVar get_goes_nc_var(NcFile *nc, const ConcatString var_name,
bool exit_if_error) {
NcVar var_data = get_nc_var(nc, var_name.c_str(), false);
NcVar var_data;
static const char *method_name = "get_goes_nc_var() -> ";
if (has_var(nc, var_name.c_str())) var_data = get_nc_var(nc, var_name.c_str(), false);
if (IS_INVALID_NC(var_data)) {
mlog << Debug(4) << method_name
<< "The variable \"" << var_name << "\" does not exist. Find \""
<< var_name.split("_")[0] << "\" variable\n";
var_data = get_nc_var(nc, var_name.split("_")[0].c_str());
}
if (IS_INVALID_NC(var_data)) {
mlog << Error << "\nget_goes_nc_var() -> "
mlog << Error << "\n" << method_name
<< "The variable \"" << var_name << "\" does not exist\n\n";
if (exit_if_error) exit(1);
}
Expand Down Expand Up @@ -2444,6 +2447,7 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
DataPlane &fr_dp, DataPlane &to_dp,
Grid fr_grid, Grid to_grid, IntArray *cellMapping, NcFile *nc_adp) {

const int log_debug_level = 4;
bool has_qc_var = false;
bool has_adp_qc_var = false;
clock_t start_clock = clock();
Expand All @@ -2465,6 +2469,10 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
// -99 is arbitrary number as invalid QC value
memset(qc_data, -99, from_data_size*sizeof(uchar));

adp_qc_high = 3; /* 3 as baseline algorithm, 0 for enterpirse algorithm */
adp_qc_medium = 1; /* 1 as baseline algorithm, 1 for enterpirse algorithm */
adp_qc_low = 0; /* 0 as baseline algorithm, 2 for enterpirse algorithm */

NcVar var_qc;
NcVar var_adp;
NcVar var_adp_qc;
Expand All @@ -2487,7 +2495,7 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
else if (is_smoke_only) var_adp = get_goes_nc_var(nc_adp, vname_smoke);

if (IS_VALID_NC(var_adp)) {
get_nc_data(&var_adp, adp_data);
get_nc_data(&var_adp, adp_data, true);

//Smoke:ancillary_variables = "DQF" ; ubyte DQF(y, x) ;
if (get_att_value_string(&var_adp, (string)"ancillary_variables", qc_var_name)) {
Expand All @@ -2497,6 +2505,7 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
has_adp_qc_var = true;
mlog << Debug(5) << method_name << "found QC var: " << qc_var_name
<< " for " << GET_NC_NAME(var_adp) << ".\n";
set_adp_gc_values(var_adp_qc);
}
}
else {
Expand Down Expand Up @@ -2557,6 +2566,15 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
int shift_bits = 2;
if (is_dust_only) shift_bits += 2;

int cnt_aod_qc_low = 0;
int cnt_aod_qc_high = 0;
int cnt_aod_qc_medium = 0;
int cnt_aod_qc_nr = 0; // no_retrieval_qf
int cnt_adp_qc_low = 0;
int cnt_adp_qc_high = 0;
int cnt_adp_qc_medium = 0;
int cnt_adp_qc_nr = 0; // no_retrieval_qf
int cnt_adp_qc_adjused = 0;
for (int xIdx=0; xIdx<to_lon_count; xIdx++) {
for (int yIdx=0; yIdx<to_lat_count; yIdx++) {
int offset = to_dp.two_to_one(xIdx,yIdx);
Expand Down Expand Up @@ -2597,14 +2615,39 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
// Filter by QC flag
if (has_qc_var || has_adp_qc_var) {
qc_value = qc_data[from_index];
if (mlog.verbosity_level() >= log_debug_level) {
if (qc_min_value > qc_value) qc_min_value = qc_value;
if (qc_max_value < qc_value) qc_max_value = qc_value;
switch (qc_value) {
case 0: cnt_aod_qc_high++; break;
case 1: cnt_aod_qc_medium++; break;
case 2: cnt_aod_qc_low++; break;
default: cnt_aod_qc_nr++; break;
}
}
if (has_adp_qc_var) {
int qc_for_flag = compute_adp_qc_flag(adp_qc_data[from_index], shift_bits);
bool filter_out = is_eq(qc_for_flag, bad_data_int);

if (mlog.verbosity_level() >= log_debug_level) {
switch (qc_for_flag) {
case 0: cnt_adp_qc_high++; break;
case 1: cnt_adp_qc_medium++; break;
case 2: cnt_adp_qc_low++; break;
default: cnt_adp_qc_nr++; break;
}
}

if (!filter_out) {
/* Adjust the quality by AOD data QC */
if (1 == qc_value && 0 == qc_for_flag) qc_for_flag = 1; /* high to medium quality */
else if (2 == qc_value) qc_for_flag = 2; /* high/medium to low quality */
if (1 == qc_value && 0 == qc_for_flag) {
qc_for_flag = 1; /* high to medium quality */
cnt_adp_qc_adjused++;
}
else if (2 == qc_value) {
qc_for_flag = 2; /* high/medium to low quality */
cnt_adp_qc_adjused++;
}

if (has_qc_flags && !qc_flags.has(qc_for_flag)) filter_out = true;
}
Expand All @@ -2617,13 +2660,8 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
qc_filtered_count++;
continue;
}

dataArray.add(data_value);
if (mlog.verbosity_level() >= 4) {
if (qc_min_value > qc_value) qc_min_value = qc_value;
if (qc_max_value < qc_value) qc_max_value = qc_value;
}
}
dataArray.add(data_value);
valid_count++;
}

Expand Down Expand Up @@ -2658,14 +2696,20 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo,
delete [] from_data;
delete [] adp_qc_data;

mlog << Debug(4) << method_name << "Count: actual: " << to_cell_count
mlog << Debug(log_debug_level) << method_name << "Count: actual: " << to_cell_count
<< ", missing: " << missing_count << ", non_missing: " << non_missing_count
<< "\n\tFiltered: by QC: " << qc_filtered_count
<< "\n Filtered: by QC: " << qc_filtered_count
<< ", by adp QC: " << adp_qc_filtered_count
<< ", by absent: " << absent_count
<< ", total: " << (qc_filtered_count + adp_qc_filtered_count + absent_count)
<< "\n\tRange: data: [" << from_min_value << " - " << from_max_value
<< "] QC: [" << qc_min_value << " - " << qc_max_value << "]\n";
<< "\n Range: data: [" << from_min_value << " - " << from_max_value
<< "] QC: [" << qc_min_value << " - " << qc_max_value << "]"
<< "\n AOD QC: high=" << cnt_aod_qc_high
<< " medium=" << cnt_aod_qc_medium << ", low=" << cnt_aod_qc_low
<< ", no_retrieval=" << cnt_aod_qc_nr
<< "\n ADP QC: high=" << cnt_adp_qc_high << " medium=" << cnt_adp_qc_medium
<< ", low=" << cnt_adp_qc_low << ", no_retrieval=" << cnt_adp_qc_nr
<< ", adjusted=" << cnt_adp_qc_adjused << "\n";

if (to_cell_count == 0) {
mlog << Warning << "\n" << method_name
Expand Down Expand Up @@ -2866,6 +2910,41 @@ void usage() {
exit(1);
}

////////////////////////////////////////////////////////////////////////
const static ConcatString att_name_values = "flag_values";
const static ConcatString att_name_meanings = "flag_meanings";

void set_adp_gc_values(NcVar var_adp_qc) {
//NcVarAtt *nc_att_values = get_nc_att(&var_adp_qc, "flag_values");
//NcVarAtt *nc_att_meanings = get_nc_att(&var_adp_qc, "flag_meanings");
//unsigned short *flag_values = new unsigned short [256];
ConcatString att_flag_meanings;

if (get_nc_att_value(&var_adp_qc, (ConcatString)"flag_meanings", att_flag_meanings)) {
StringArray flag_meanings = to_lower(att_flag_meanings).split(" ");
unsigned short flag_values[256];
for (int i=0; i<flag_meanings.n(); i++) flag_values[i] = (unsigned short)-1;

if (get_nc_att_values(&var_adp_qc, att_name_values, flag_values)) {
int idx;
StringArray flag_meanings = to_lower(att_flag_meanings).split(" ");
if (flag_meanings.has("low_confidence_smoke_detection_qf", idx)) {
adp_qc_low = (flag_values[idx] >> 2) & 0x03;
}
if (flag_meanings.has("medium_confidence_smoke_detection_qf", idx)) {
adp_qc_medium = (flag_values[idx] >> 2) & 0x03;
}
if (flag_meanings.has("high_confidence_smoke_detection_qf", idx)) {
adp_qc_high = (flag_values[idx] >> 2) & 0x03;
}
}
mlog << Debug(4) << "set_adp_gc_values() "
<< " high_confidence = " << adp_qc_high
<< ", medium_confidence = " << adp_qc_medium
<< ", low_confidence = " << adp_qc_low << "\n";
}
}

////////////////////////////////////////////////////////////////////////

void set_field(const StringArray &a) {
Expand Down

0 comments on commit e7e25a8

Please sign in to comment.