From eb35a2ea45f5d1fa9034e96e6af804b9fd03f000 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Fri, 3 May 2024 14:13:43 -0600 Subject: [PATCH] Bugfix 2867 point2grid qc flag main v11.1 (#2874) * #2867 Added point2grid_GOES_16_ADP_Enterprise_high and changed qc flags for ADP * #2867 Added get_nc_att_values_ * #2867 Added get_nc_att_values * #2867 Get the ADP QC flag values from the varibale attributes (support GOES16 Enterprise allgorithm) and apply them for QC Flags. Adjusted the ADP QC flags based on the variable QC values --------- Co-authored-by: Howard Soh --- internal/test_unit/xml/unit_point2grid.xml | 29 +++- src/libcode/vx_nc_util/nc_utils.cc | 28 ++- src/libcode/vx_nc_util/nc_utils.h | 5 +- src/libcode/vx_nc_util/nc_utils.hpp | 23 +++ src/tools/other/point2grid/point2grid.cc | 189 ++++++++++++++++----- 5 files changed, 218 insertions(+), 56 deletions(-) diff --git a/internal/test_unit/xml/unit_point2grid.xml b/internal/test_unit/xml/unit_point2grid.xml index dd1792ccf0..5221ea4460 100644 --- a/internal/test_unit/xml/unit_point2grid.xml +++ b/internal/test_unit/xml/unit_point2grid.xml @@ -137,7 +137,7 @@ &OUTPUT_DIR;/point2grid/point2grid_GOES_16_AOD_TO_G212_compute.nc - + &MET_BIN;/point2grid @@ -155,7 +155,7 @@ &OUTPUT_DIR;/point2grid/point2grid_GOES_16_AOD_TO_G212_gaussian.nc - + &MET_BIN;/point2grid @@ -174,7 +174,7 @@ &OUTPUT_DIR;/point2grid/point2grid_GOES_16_ADP.nc - + &MET_BIN;/point2grid @@ -185,7 +185,7 @@ G212 \ &OUTPUT_DIR;/point2grid/point2grid_GOES_16_AOD_TO_G212_grid_map.nc \ -field 'name="AOD"; level="(*,*)";' \ - -qc 1,2,3 -method MAX \ + -qc 0,1,2 -method MAX \ -v 1 @@ -205,7 +205,7 @@ G212 \ &OUTPUT_DIR;/point2grid/point2grid_GOES_16_AOD_TO_G212.nc \ -field 'name="AOD"; level="(*,*)";' \ - -qc 1,2,3 -method MAX \ + -qc 0,1,2 -method MAX \ -v 1 @@ -213,6 +213,25 @@ + + &MET_BIN;/point2grid + + MET_TMP_DIR &OUTPUT_DIR;/point2grid + + \ + &DATA_DIR_MODEL;/goes_16/OR_ABI-L2-AODC-M6_G16_s20241100001171_e20241100003544_c20241100006242.nc \ + G212 \ + &OUTPUT_DIR;/point2grid/point2grid_GOES_16_ADP_Enterprise_high.nc \ + -field 'name="AOD_Smoke"; level="(*,*)";' \ + -adp &DATA_DIR_MODEL;/goes_16/OR_ABI-L2-ADPC-M6_G16_s20241100001171_e20241100003544_c20241100006361.nc \ + -qc 0,1 -method MAX \ + -v 1 + + + &OUTPUT_DIR;/point2grid/point2grid_GOES_16_ADP_Enterprise_high.nc + + + &MET_BIN;/point2grid diff --git a/src/libcode/vx_nc_util/nc_utils.cc b/src/libcode/vx_nc_util/nc_utils.cc index 6350b6576a..83a6f47950 100644 --- a/src/libcode/vx_nc_util/nc_utils.cc +++ b/src/libcode/vx_nc_util/nc_utils.cc @@ -537,6 +537,17 @@ bool get_nc_att_value(const NcVar *var, const ConcatString &att_name, //////////////////////////////////////////////////////////////////////// +bool get_nc_att_values(const NcVar *var, const ConcatString &att_name, + unsigned short *att_val, bool exit_on_error) { + static const char *method_name = "get_nc_att_value(NcVar,float) -> "; + bool status = get_nc_att_values_(var, att_name, att_val, exit_on_error, + (unsigned short)bad_data_int, method_name); + return status; +} + +//////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////// + bool get_nc_att_value(const NcVarAtt *att, ConcatString &att_val) { bool status = false; @@ -1729,18 +1740,21 @@ bool get_nc_data(NcVar *var, char **data) { //////////////////////////////////////////////////////////////////////// -bool get_nc_data(NcVar *var, uchar *data) { +bool get_nc_data(NcVar *var, uchar *data, bool allow_conversion) { bool return_status = false; int data_type = GET_NC_TYPE_ID_P(var); static const char *method_name = "get_nc_data(NcVar *, uchar *) -> "; if (NC_UBYTE == data_type) return_status = get_nc_data_t(var, data); - else if (NC_BYTE == data_type && has_unsigned_attribute(var)) { - int cell_count = get_data_size(var); - ncbyte *signed_data = new ncbyte[cell_count]; - return_status = get_nc_data_t(var, signed_data); - for (int idx=0; idx +bool get_nc_att_values_(const netCDF::NcVar *var, const ConcatString &att_name, + T *att_vals, bool exit_on_error, + T bad_data, const char *caller_name) { + // caller should initialize att_vals + + // + // Retrieve the NetCDF variable attribute. + // + netCDF::NcVarAtt *att = get_nc_att(var, att_name); + bool status = IS_VALID_NC_P(att); + if (status) att->getValues(att_vals); + else { + mlog << Error << "\n" << caller_name + << get_log_msg_for_att(att, GET_SAFE_NC_NAME_P(var), att_name); + } + if (att) delete att; + if (!status && exit_on_error) exit(1); + return status; +} + +//////////////////////////////////////////////////////////////////////// + template bool get_nc_att_value_(const netCDF::NcVarAtt *att, T &att_val, bool exit_on_error, T bad_data, const char *caller_name) { diff --git a/src/tools/other/point2grid/point2grid.cc b/src/tools/other/point2grid/point2grid.cc index 1861a1ca66..edf18accaa 100644 --- a/src/tools/other/point2grid/point2grid.cc +++ b/src/tools/other/point2grid/point2grid.cc @@ -124,6 +124,9 @@ static NcFile *nc_out = (NcFile *) 0; static NcDim lat_dim ; static NcDim lon_dim ; +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 */ //////////////////////////////////////////////////////////////////////// @@ -170,6 +173,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 @@ -940,8 +944,6 @@ void process_point_met_data(MetPointData *met_point_obs, MetConfig &config, VarI int valid_count = 0; int absent_count = 0; int censored_count = 0; - int qc_filtered_count = 0; - int adp_qc_filtered_count = 0; float data_value; float from_min_value = 10e10; float from_max_value = -10e10; @@ -1185,7 +1187,6 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, void process_point_python(string python_command, MetConfig &config, VarInfo *vinfo, const Grid to_grid, bool use_xarray) { - int idx, hdr_idx; ConcatString vname, vname_cnt, vname_mask; DataPlane fr_dp, to_dp; DataPlane cnt_dp, mask_dp; @@ -1193,8 +1194,6 @@ void process_point_python(string python_command, MetConfig &config, VarInfo *vin NcVar var_obs_gc, var_obs_var; clock_t start_clock = clock(); - bool has_prob_thresh = !prob_cat_thresh.check(bad_data_double); - unixtime requested_valid_time, valid_time; static const char *method_name = "process_point_python() -> "; static const char *method_name_s = "process_point_python()"; @@ -1863,6 +1862,24 @@ void check_lat_lon(int data_size, float *latitudes, float *longitudes) { << method_name << "LONG: " << min_lon << " to " << max_lon << "\n"; } +//////////////////////////////////////////////////////////////////////// +// QC flags: 0=high, 1=medium, 2=low +// Enterpise algorithm: 0=high, 1=medium, 2=low +// Baseline algorithm: 3=high, 1=medium, 0=low (high=12/48, medium=4/16) +// returns bad_data_int if it does not belong to high, mediuam, or low. + +int compute_adp_qc_flag(int adp_qc, int shift_bits) { + int particle_qc = ((adp_qc >> shift_bits) & 0x03); + 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; +} + //////////////////////////////////////////////////////////////////////// static unixtime compute_unixtime(NcVar *time_var, unixtime var_value) { @@ -2342,12 +2359,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)) { - var_data = get_nc_var(nc, var_name.split("_")[0].c_str()); + 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); } @@ -2424,6 +2446,7 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, bool has_qc_var = false; bool has_adp_qc_var = false; + const int log_debug_level = 4; clock_t start_clock = clock(); int to_lat_count = to_grid.ny(); int to_lon_count = to_grid.nx(); @@ -2443,6 +2466,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; @@ -2465,13 +2492,14 @@ 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)) { var_adp_qc = get_nc_var(nc_adp, qc_var_name.c_str()); if (IS_VALID_NC(var_adp_qc)) { get_nc_data(&var_adp_qc, adp_qc_data); + set_adp_gc_values(var_adp_qc); has_adp_qc_var = true; mlog << Debug(5) << method_name << "found QC var: " << qc_var_name << " for " << GET_NC_NAME(var_adp) << ".\n"; @@ -2521,7 +2549,6 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, int non_missing_count = 0; int qc_filtered_count = 0; int adp_qc_filtered_count = 0; - float data_value; float from_min_value = 10e10; float from_max_value = -10e10; float qc_min_value = 10e10; @@ -2534,6 +2561,18 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, missing_count = non_missing_count = 0; to_dp.set_constant(bad_data_double); + 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; xIdxcensor_thresh().n(); i++) { + // Break out after the first match. + if(vinfo->censor_thresh()[i].check(data_value)) { + data_value = vinfo->censor_val()[i]; + censored_count++; + break; + } + } + + // Check the data existance (always 1 if ADP variable does not exist) + if (0 == adp_data[from_index]) { + absent_count++; + continue; + } + // Filter by QC flag - qc_value = qc_data[from_index]; - if (!has_qc_var || !has_qc_flags || qc_flags.has(qc_value)) { - for(int i=0; icensor_thresh().n(); i++) { - // Break out after the first match. - if(vinfo->censor_thresh()[i].check(data_value)) { - data_value = vinfo->censor_val()[i]; - censored_count++; - break; + 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 (0 == adp_data[from_index]) { - absent_count++; - continue; - } + 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 (has_adp_qc_var && has_qc_flags) { - int shift_bits = 2; - if (is_dust_only) shift_bits += 2; - particle_qc = ((adp_qc_data[from_index] >> shift_bits) & 0x03); - int qc_for_flag = 3 - particle_qc; // high = 3, qc_flag for high = 0 - if (!qc_flags.has(qc_for_flag)) { + 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 */ + 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; + } + if (filter_out) { adp_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; + else if (has_qc_var && has_qc_flags && !qc_flags.has(qc_value)) { + qc_filtered_count++; + continue; } } - else { - qc_filtered_count++; - } + dataArray.add(data_value); valid_count++; } @@ -2618,7 +2686,6 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, << data_count << " data values.\n"; } } - else {} } } @@ -2627,14 +2694,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 @@ -2835,6 +2908,38 @@ 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) { + 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> 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) {