Skip to content

Commit

Permalink
Merge pull request #2719 from dtcenter/bugfix_2687_hpbl
Browse files Browse the repository at this point in the history 
Bugfix 2687 hpbl
  • Loading branch information
@hsoh-u
hsoh-u authored Oct 23, 2023
2 parents 9a71af5 + b916a30 commit 9fed892
Show file tree
Hide file tree
Showing 10 changed files with 230 additions and 221 deletions.
49 changes: 25 additions & 24 deletions scripts/python/met/dataplane.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
1 change: 0 additions & 1 deletion src/basic/vx_util/ascii_table.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1567,7 +1567,6 @@ switch ( just ) {
default:
mlog << Error << "\njustified_item() -> bad justification value\n\n";
exit ( 1 );
break;

} // switch

Expand Down
14 changes: 7 additions & 7 deletions src/libcode/vx_nc_util/nc_utils.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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
Expand All @@ -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];
Expand All @@ -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 ("
Expand Down
76 changes: 40 additions & 36 deletions src/libcode/vx_statistics/compute_stats.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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];
Expand All @@ -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];
Expand Down Expand Up @@ -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;
}

////////////////////////////////////////////////////////////////////////
5 changes: 3 additions & 2 deletions src/libcode/vx_statistics/compute_stats.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

////////////////////////////////////////////////////////////////////////
//
Expand Down
Loading

0 comments on commit 9fed892

Please sign in to comment.