Bugfix 2687 hpbl

scripts/python/met/dataplane.py

@@ -180,32 +180,33 @@ def validate_met_data(met_data, fill_value=None):
       from_ndarray = False
       if met_data is None:
          logger.quit(f"{method_name} The met_data is None")
-
-      nx, ny = met_data.shape
-      met_fill_value = dataplane.MET_FILL_VALUE
-      if dataplane.is_xarray_dataarray(met_data):
-         from_xarray = True
-         attrs = met_data.attrs
-         met_data = met_data.data
-         modified_met_data = True
-      if isinstance(met_data, np.ndarray):
-         from_ndarray = True
-         met_data = np.ma.array(met_data)
-
-      if isinstance(met_data, np.ma.MaskedArray):
-         is_int_data = dataplane.is_integer(met_data[0,0]) or dataplane.is_integer(met_data[int(nx/2),int(ny/2)])
-         met_data = np.ma.masked_equal(met_data, float('nan'))
-         met_data = np.ma.masked_equal(met_data, float('inf'))
-         if fill_value is not None:
-            met_data = np.ma.masked_equal(met_data, fill_value)
-         met_data = met_data.filled(int(met_fill_value) if is_int_data else met_fill_value)
       else:
-         logger.log_msg(f"{method_name} unknown datatype {type(met_data)}")
+         nx, ny = met_data.shape
+
+         met_fill_value = dataplane.MET_FILL_VALUE
+         if dataplane.is_xarray_dataarray(met_data):
+            from_xarray = True
+            attrs = met_data.attrs
+            met_data = met_data.data
+            modified_met_data = True
+         if isinstance(met_data, np.ndarray):
+            from_ndarray = True
+            met_data = np.ma.array(met_data)
+
+         if isinstance(met_data, np.ma.MaskedArray):
+            is_int_data = dataplane.is_integer(met_data[0,0]) or dataplane.is_integer(met_data[int(nx/2),int(ny/2)])
+            met_data = np.ma.masked_equal(met_data, float('nan'))
+            met_data = np.ma.masked_equal(met_data, float('inf'))
+            if fill_value is not None:
+               met_data = np.ma.masked_equal(met_data, fill_value)
+            met_data = met_data.filled(int(met_fill_value) if is_int_data else met_fill_value)
+         else:
+            logger.log_msg(f"{method_name} unknown datatype {type(met_data)}")
 
-      if dataplane.KEEP_XARRAY:
-         return xr.DataArray(met_data,attrs=attrs) if from_xarray else met_data
-      else:
-         return met_data
+         if dataplane.KEEP_XARRAY:
+            return xr.DataArray(met_data,attrs=attrs) if from_xarray else met_data
+         else:
+            return met_data
 
 
 def main(argv):

src/basic/vx_util/ascii_table.cc

@@ -1567,7 +1567,6 @@ switch ( just )  {
    default:
       mlog << Error << "\njustified_item() -> bad justification value\n\n";
       exit ( 1 );
-      break;
 
 }   //  switch
 

src/libcode/vx_nc_util/nc_utils.cc

@@ -8,21 +8,21 @@
 
 ////////////////////////////////////////////////////////////////////////
 
-using namespace std;
-
 #include <string.h>
 #include <cstring>
 #include <sys/stat.h>
 
 #include <netcdf>
-using namespace netCDF;
-using namespace netCDF::exceptions;
 
 #include "vx_log.h"
 #include "nc_utils.h"
 #include "util_constants.h"
 #include "vx_cal.h"
 
+using namespace std;
+using namespace netCDF;
+using namespace netCDF::exceptions;
+
 ////////////////////////////////////////////////////////////////////////
 
 void patch_nc_name(string *var_name) {
@@ -172,7 +172,7 @@ bool get_att_value_chars(const NcAtt *att, ConcatString &value) {
             att->getValues(att_value);
             value = att_value;
          }
-         catch (exceptions::NcChar ex) {
+         catch (exceptions::NcChar &ex) {
             value = "";
             // Handle netCDF::exceptions::NcChar:  NetCDF: Attempt to convert between text & numbers
             mlog << Warning << "\n" << method_name
@@ -188,7 +188,7 @@ bool get_att_value_chars(const NcAtt *att, ConcatString &value) {
             att->getValues(att_value);
             value = att_value;
          }
-         catch (exceptions::NcChar ex) {
+         catch (exceptions::NcChar &ex) {
             int num_elements_sub = 8096;
             int num_elements = att->getAttLength();;
             char *att_value[num_elements];
@@ -199,7 +199,7 @@ bool get_att_value_chars(const NcAtt *att, ConcatString &value) {
                att->getValues(att_value);
                value = att_value[0];
             }
-            catch (exceptions::NcException ex) {
+            catch (exceptions::NcException &ex) {
                mlog << Warning << "\n" << method_name
                     << "Exception: " << ex.what() << "\n"
                     << "Fail to read " << GET_NC_NAME_P(att) << " attribute ("

src/libcode/vx_statistics/compute_stats.cc

@@ -1449,49 +1449,55 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps) {
       seeps_mprs.push_back(seeps_mpr);
    }
    if (count > 0) {
-      double *density_vector;
+      vector<double> density_vector;
       double pvf[SEEPS_MATRIX_SIZE];
       double weighted_score, weight_sum, weight[count];
 
       seeps->n_obs = count;
       seeps->mean_fcst = fcst_sum / count;
       seeps->mean_obs = obs_sum / count;
       seeps->score = score_sum / count;
-      density_vector = compute_seeps_density_vector(pd, seeps);
 
       weighted_score = 0.;
       for (int i=0; i<SEEPS_MATRIX_SIZE; i++) pvf[i] = 0.;
-      if (density_vector != nullptr) {
-         //IDL: w = 1/d
-         weight_sum = 0.;
-         for (int i=0; i<count; i++) {
-            if (is_eq(density_vector[i], 0)) weight[i] = 0;
-            else {
-               weight[i] = 1 / density_vector[i];
-               weight_sum += weight[i];
-            }
+
+      compute_seeps_density_vector(pd, seeps, density_vector);
+      int density_cnt = density_vector.size();
+
+      //IDL: w = 1/d
+      weight_sum = 0.;
+      for (int i=0; i<count; i++) weight[i] = 0;
+      for (int i=0; i<density_cnt; i++) {
+         if (!is_eq(density_vector[i], 0)) {
+            weight[i] = 1 / density_vector[i];
+            weight_sum += weight[i];
          }
-         if (!is_eq(weight_sum, 0)) {
-            //IDL: w = w/sum(w)
-            for (int i=0; i<count; i++) weight[i] /= weight_sum;
-
-            //IDL: for i=0l, n-1l do begin
-            for (int i=0; i<seeps_mprs.size(); i++) {
-               seeps_mpr = seeps_mprs[i];
-               //IDL: s = s + c(4+cat(i) * w{i)
-               if (i < count) {
-                  weighted_score += seeps_mpr->score * weight[i];
-                  //IDL: svf(cat{i)) = svf(cat{i)) + c(4+cat(i) * w{i)
-                  //IDL: pvf(cat{i)) = pvf(cat{i)) + w{i)
-                  pvf[seeps_mpr->s_idx] += weight[i];
-               }
-               else mlog << Debug(1) << method_name
-                         << "the length of density vector (" << count << ") is less than MPR.\n";
+      }
+      if (!is_eq(weight_sum, 0)) {
+         //IDL: w = w/sum(w)
+         for (int i=0; i<count; i++) weight[i] /= weight_sum;
+
+         //IDL: for i=0l, n-1l do begin
+         for (int i=0; i<seeps_mprs.size(); i++) {
+            seeps_mpr = seeps_mprs[i];
+            //IDL: s = s + c(4+cat(i) * w{i)
+            if (i < density_cnt) {
+               weighted_score += seeps_mpr->score * weight[i];
+               //IDL: svf(cat{i)) = svf(cat{i)) + c(4+cat(i) * w{i)
+               //IDL: pvf(cat{i)) = pvf(cat{i)) + w{i)
+               pvf[seeps_mpr->s_idx] += weight[i];
+            }
+            else {
+               mlog << Debug(1) << method_name
+                    << "the length of density vector (" << density_cnt
+                    << ") is less than SEEPS MPR (" << seeps_mprs.size() << ").\n";
+               break;
             }
          }
-
-         if (density_vector != nullptr) delete [] density_vector;
       }
+
+      density_vector.clear();
+
       seeps_mprs.clear();
 
       // The weight for s12 to s32 should come from climo file, but not available yet
@@ -1707,7 +1713,7 @@ void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &ob
 // ; PV-WAVE prints: 2.00000      4.00000
 ////////////////////////////////////////////////////////////////////////
 
-double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps) {
+void compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps, vector<double> &density_vector) {
    int seeps_idx;
    SeepsScore *seeps_mpr;
    int seeps_cnt = seeps->n_obs;
@@ -1725,12 +1731,11 @@ double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *see
    if (seeps_cnt == 0) {
       mlog << Debug(1) << method_name
            << "no SEEPS_MPR available.\n";
-      return nullptr;
+      return;
    }
 
    // Get lat/lon & convert them to radian and get sin/cos values
    seeps_idx = 0;
-   double *density_vector = new double[seeps_cnt];
    for(int i=0; i<pd->n_obs; i++) {
       if (i >= pd->seeps_mpr.size()) break;
       seeps_mpr = pd->seeps_mpr[i];
@@ -1752,10 +1757,9 @@ double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *see
    // Initialize
    v_count = 0;
    if (seeps_idx < seeps_cnt) seeps_cnt = seeps_idx;
-   for(int i=0; i<seeps_cnt; i++) {
-      density_vector[i] = 0.;
-   }
+   density_vector.reserve(seeps_cnt);
    for(int j=0; j<seeps_cnt; j++) {
+      density_vector.push_back(0.);
       for(int i=0; i<seeps_cnt; i++) {
          // Make n by n matrix: clat_m = clat#transpose(clat), slat_m = slat#transpose(slat) by IDL
          clat_m[i][j] = clat[i] * clat[j];
@@ -1793,7 +1797,7 @@ double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *see
            << " density_vector[" << (seeps_cnt-1) << "]=" << density_vector[seeps_cnt-1] << "\n";
    }
 
-   return density_vector;
+   return;
 }
 
 ////////////////////////////////////////////////////////////////////////

src/libcode/vx_statistics/compute_stats.h

@@ -60,12 +60,13 @@ extern void   compute_i_mean_stdev(const NumArray &,
                                    bool, double, int,
                                    CIInfo &, CIInfo &);
 
-extern void   compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps);
-extern double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps);
+extern void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps);
 extern void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &obs_dp,
                                           DataPlane &seeps_dp, DataPlane &seeps_dp_fcat,
                                           DataPlane &seeps_dp_ocat,SeepsAggScore *seeps,
                                           int month, int hour, const SingleThresh &seeps_p1_thresh);
+extern void compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps,
+                                         std::vector<double> &density_vector);
 
 ////////////////////////////////////////////////////////////////////////
 //