diff --git a/.github/ISSUE_TEMPLATE/enhancement_request.md b/.github/ISSUE_TEMPLATE/enhancement_request.md index e1134781cd..78cd208f15 100644 --- a/.github/ISSUE_TEMPLATE/enhancement_request.md +++ b/.github/ISSUE_TEMPLATE/enhancement_request.md @@ -39,8 +39,8 @@ Consider breaking the enhancement down into sub-issues. - [ ] Select **requestor(s)** ### Milestone and Projects ### -- [ ] Select **Milestone** as the next official version or **Backlog of Development Ideas** -- [ ] For the next official version, select the **MET-X.Y.Z Development** project +- [ ] Select **Milestone** as a **MET-X.Y.Z** version, **Consider for Next Release**, or **Backlog of Development Ideas** +- [ ] For a **MET-X.Y.Z** version, select the **MET-X.Y.Z Development** project ## Define Related Issue(s) ## Consider the impact to the other METplus components. diff --git a/.github/ISSUE_TEMPLATE/new_feature_request.md b/.github/ISSUE_TEMPLATE/new_feature_request.md index e4adb302c7..ea95d70c13 100644 --- a/.github/ISSUE_TEMPLATE/new_feature_request.md +++ b/.github/ISSUE_TEMPLATE/new_feature_request.md @@ -43,8 +43,8 @@ Consider breaking the new feature down into sub-issues. - [ ] Select **requestor(s)** ### Milestone and Projects ### -- [ ] Select **Milestone** as the next official version or **Backlog of Development Ideas** -- [ ] For the next official version, select the **MET-X.Y.Z Development** project +- [ ] Select **Milestone** as a **MET-X.Y.Z** version, **Consider for Next Release**, or **Backlog of Development Ideas** +- [ ] For a **MET-X.Y.Z** version, select the **MET-X.Y.Z Development** project ## Define Related Issue(s) ## Consider the impact to the other METplus components. diff --git a/.github/ISSUE_TEMPLATE/sub-issue.md b/.github/ISSUE_TEMPLATE/sub-issue.md index e02109e963..45ee00432d 100644 --- a/.github/ISSUE_TEMPLATE/sub-issue.md +++ b/.github/ISSUE_TEMPLATE/sub-issue.md @@ -29,5 +29,5 @@ This is a sub-issue of #*List the parent issue number here*. - [ ] Select **requestor(s)** ### Milestone and Projects ### -- [ ] Select **Milestone** as the next official version or **Backlog of Development Ideas** -- [ ] For the next official version, select the **MET-X.Y.Z Development** project +- [ ] Select **Milestone** as a **MET-X.Y.Z** version, **Consider for Next Release**, or **Backlog of Development Ideas** +- [ ] For a **MET-X.Y.Z** version, select the **MET-X.Y.Z Development** project diff --git a/.github/ISSUE_TEMPLATE/task.md b/.github/ISSUE_TEMPLATE/task.md index 05ad99baf8..33685adf25 100644 --- a/.github/ISSUE_TEMPLATE/task.md +++ b/.github/ISSUE_TEMPLATE/task.md @@ -39,8 +39,8 @@ Consider breaking the task down into sub-issues. - [ ] Select **requestor(s)** ### Milestone and Projects ### -- [ ] Select **Milestone** as the next official version or **Backlog of Development Ideas** -- [ ] For the next official version, select the **MET-X.Y.Z Development** project +- [ ] Select **Milestone** as a **MET-X.Y.Z** version, **Consider for Next Release**, or **Backlog of Development Ideas** +- [ ] For a **MET-X.Y.Z** version, select the **MET-X.Y.Z Development** project ## Define Related Issue(s) ## Consider the impact to the other METplus components. diff --git a/.github/jobs/set_job_controls.sh b/.github/jobs/set_job_controls.sh index 76a123d18b..b8a66840a6 100755 --- a/.github/jobs/set_job_controls.sh +++ b/.github/jobs/set_job_controls.sh @@ -57,6 +57,12 @@ elif [ "${GITHUB_EVENT_NAME}" == "push" ]; then input_data_version=${branch_name:6} fi + # check for main_vX.Y in the branch name + elif [[ "${branch_name}" =~ .*(main_v)([0-9]+\.[0-9]+).* ]]; then + + truth_data_version=${BASH_REMATCH[1]}${BASH_REMATCH[2]} + input_data_version=${BASH_REMATCH[2]} + fi # check commit messages for skip or force keywords @@ -87,14 +93,27 @@ elif [ "${GITHUB_EVENT_NAME}" == "push" ]; then elif [ "${GITHUB_EVENT_NAME}" == "workflow_dispatch" ]; then + branch_name=`cut -d "/" -f3 <<< "${GITHUB_REF}"` + + # check for main_vX.Y in the branch name + if [[ "${branch_name}" =~ .*(main_v)([0-9]+\.[0-9]+).* ]]; then + + truth_data_version=${BASH_REMATCH[1]}${BASH_REMATCH[2]} + input_data_version=${BASH_REMATCH[2]} + + fi + if [ "${force_tests}" == "true" ]; then + run_diff=true + fi fi # if updating truth or running diff, run unit tests -if [ "$run_update_truth" == "true" ] || [ "$run_diff" == "true" ]; then +if [ "$run_update_truth" == "true" ] || + [ "$run_diff" == "true" ]; then run_unit_tests=true diff --git a/docs/Users_Guide/release-notes.rst b/docs/Users_Guide/release-notes.rst index 30a115abc1..01aee88695 100644 --- a/docs/Users_Guide/release-notes.rst +++ b/docs/Users_Guide/release-notes.rst @@ -9,6 +9,29 @@ When applicable, release notes are followed by the GitHub issue number which des enhancement, or new feature (`MET GitHub issues `_). Important issues are listed **in bold** for emphasis. +MET Version 11.1.1 release notes (20240621) +------------------------------------------- + + .. dropdown:: Bugfixes + + * Bugfix: Refine support for coordinate dimensions in CF-compliant NetCDF files (`#2638 `_). + * Bugfix: Fix support for NSIDC v4 Climate Data Record data on Polar Stereographic grids in CF-compliant NetCDF files (`#2652 `_). + * Bugfix: Fix the Point-Stat CNT header line typo causing duplicate SI_BCL column names (`#2730 `_). + * Bugfix: Fix MET to compile without the optional --enable-python configuration option (`#2760 `_). + * Bugfix: Fix the parsing of level values for GRIB2 template 4.48 data (`#2782 `_). + * **Bugfix: Fix the TC-Diag and TC-RMW tools to correctly handle the range and azimuth settings in range/azimuth grids** (`#2833 `_). + * **Bugfix: Fix TC-RMW to correct the tangential and radial wind computations** (`#2841 `_). + * Bugfix: Fix Ensemble-Stat's handling of climo data when verifying ensemble-derived probabilities (`#2856 `_). + * Bugfix: Update the logic of -qc option for ADP at point2grid (`#2867 `_). + * Bugfix: Some AOD data was filtered out without -qc option with point2grid (`#2884 `_). + * Bugfix: Fix inconsistent handling of point observation valid times processed through Python embedding (`#2897 `_). + * Bugfix: Fix PBL derivation bug in main_v11.1 by porting over fix from develop (`#2902 `_). + + .. dropdown:: Enhancements + + * **Eliminate the use of temporary files in the vx_config library** (`#2691 `_). + * Enhance Point2Grid to support modified quality control settings for smoke/dust AOD data in GOES-16/17 as of April 16, 2024 (`#2853 `_). + MET Version 11.1.0 release notes (20230731) ------------------------------------------- diff --git a/docs/Users_Guide/tc-diag.rst b/docs/Users_Guide/tc-diag.rst index 53c1ee1b0d..1ef12703d2 100644 --- a/docs/Users_Guide/tc-diag.rst +++ b/docs/Users_Guide/tc-diag.rst @@ -24,7 +24,7 @@ Originally developed for the Statistical Hurricane Intensity Prediction Scheme ( TC-Diag is run once for each initialization time to produce diagnostics for each user-specified combination of TC tracks and model fields. The user provides track data (such as one or more ATCF a-deck track files), along with track filtering criteria as needed, to select one or more tracks to be processed. The user also provides gridded model data from which diagnostics should be computed. Gridded data can be provided for multiple concurrent storms, multiple models, and/or multiple domains (i.e. parent and nest) in a single run. -TC-Diag first determines the list of valid times that appear in any one of the tracks. For each valid time, it processes all track points for that time. For each track point, it reads the gridded model fields requested in the configuration file and transforms the gridded data to a range-azimuth cylindrical coordinates grid. For each domain, it writes the range-azimuth data to a temporary NetCDF file. +TC-Diag first determines the list of valid times that appear in any one of the tracks. For each valid time, it processes all track points for that time. For each track point, it reads the gridded model fields requested in the configuration file and transforms the gridded data to a range-azimuth cylindrical coordinates grid, as described for the TC-RMW tool in :numref:`tc-rmw`. For each domain, it writes the range-azimuth data to a temporary NetCDF file. .. note:: The current version of the tool does not yet include the capabilities described in the next three paragraphs. These additional capabilities are planned to be added in the MET v12.0.0 release later in 2023. diff --git a/docs/Users_Guide/tc-rmw.rst b/docs/Users_Guide/tc-rmw.rst index d37ea66e42..9c67d56312 100644 --- a/docs/Users_Guide/tc-rmw.rst +++ b/docs/Users_Guide/tc-rmw.rst @@ -7,7 +7,7 @@ TC-RMW Tool Introduction ============ -The TC-RMW tool regrids tropical cyclone model data onto a moving range-azimuth grid centered on points along the storm track provided in ATCF format, most likely the adeck generated from the file. The radial grid spacing may be set as a factor of the radius of maximum winds (RMW). If wind fields are specified in the configuration file the radial and tangential wind components will be computed. Any regridding method available in MET can be used to interpolate data on the model output grid to the specified range-azimuth grid. The regridding will be done separately on each vertical level. The model data files must coincide with track points in a user provided ATCF formatted track file. +The TC-RMW tool regrids tropical cyclone model data onto a moving range-azimuth grid centered on points along the storm track provided in ATCF format, most likely the adeck generated from the file. The radial grid spacing can be defined in kilometers or as a factor of the radius of maximum winds (RMW). The azimuthal grid spacing is defined in degrees clockwise from due east. If wind vector fields are specified in the configuration file the radial and tangential wind components will be computed. Any regridding method available in MET can be used to interpolate data on the model output grid to the specified range-azimuth grid. The regridding will be done separately on each vertical level. The model data files must coincide with track points in a user provided ATCF formatted track file. Practical information ===================== diff --git a/docs/conf.py b/docs/conf.py index dc186985bc..e703398f92 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -20,11 +20,11 @@ project = 'MET' author = 'UCAR/NCAR, NOAA, CSU/CIRA, and CU/CIRES' author_list = 'Jensen, T., J. Prestopnik, H. Soh, L. Goodrich, B. Brown, R. Bullock, J. Halley Gotway, K. Newman, J. Opatz' -version = '11.1.0' +version = '11.1.1' verinfo = version release = f'{version}' -release_year = '2023' -release_date = f'{release_year}-07-31' +release_year = '2024' +release_date = f'{release_year}-06-25' copyright = f'{release_year}, {author}' # -- General configuration --------------------------------------------------- diff --git a/internal/scripts/environment/development.seneca b/internal/scripts/environment/development.seneca index 69f82ae047..c6f249befd 100644 --- a/internal/scripts/environment/development.seneca +++ b/internal/scripts/environment/development.seneca @@ -44,7 +44,7 @@ export MET_TEST_INPUT=${MET_PROJ_DIR}/MET_test_data/unit_test export MET_FONT_DIR=${MET_TEST_INPUT}/fonts # Define Rscript to use a version with the ncdf4 package 1.17 or later -export MET_TEST_RSCRIPT=/nrit/ral/R-4.3.1/bin/Rscript +export MET_TEST_RSCRIPT=/nrit/ral/R-4.4.0/bin/Rscript # Define runtime Python version export MET_TEST_MET_PYTHON_EXE=${MET_PYTHON_BIN_EXE} diff --git a/src/basic/vx_util/util_constants.h b/src/basic/vx_util/util_constants.h index 4014cd3a27..1d2ad91907 100644 --- a/src/basic/vx_util/util_constants.h +++ b/src/basic/vx_util/util_constants.h @@ -18,6 +18,7 @@ //////////////////////////////////////////////////////////////////////// // Released versions of MET +static const char met_version_11_1_1[] = "V11.1.1"; static const char met_version_11_1_0[] = "V11.1.0"; static const char met_version_11_0_0[] = "V11.0.0"; static const char met_version_10_1_0[] = "V10.1.0"; @@ -42,7 +43,7 @@ static const char met_version_1_1[] = "V1.1"; //////////////////////////////////////////////////////////////////////// -static const char * const met_version = met_version_11_1_0; +static const char * const met_version = met_version_11_1_1; static const char default_met_data_dir[] = "MET_BASE"; static const char txt_file_ext[] = ".txt"; static const char stat_file_ext[] = ".stat"; diff --git a/src/libcode/vx_grid/earth_rotation.cc b/src/libcode/vx_grid/earth_rotation.cc index 35049e0076..ec8b855379 100644 --- a/src/libcode/vx_grid/earth_rotation.cc +++ b/src/libcode/vx_grid/earth_rotation.cc @@ -179,7 +179,13 @@ M23 = clat*clon; M33 = slat; */ -set_np(lat_center, lon_center, lon_center - 180.0); + // + // MET #2841 define the rotation by subtracting 90 degrees + // instead of 180 to define TCRMW grids as pointing east + // instead of north. + // + +set_np(lat_center, lon_center, lon_center - 90.0); // // diff --git a/src/libcode/vx_grid/tcrmw_grid.cc b/src/libcode/vx_grid/tcrmw_grid.cc index 7ab34aa415..d0c32d4e51 100644 --- a/src/libcode/vx_grid/tcrmw_grid.cc +++ b/src/libcode/vx_grid/tcrmw_grid.cc @@ -136,7 +136,7 @@ Ir.set_xyz(1.0, 0.0, 0.0); Jr.set_xyz(0.0, 1.0, 0.0); Kr.set_xyz(0.0, 0.0, 1.0); -Range_n = 0; +Range_n = 0; Azimuth_n = 0; Range_max_km = 0.0; @@ -288,6 +288,7 @@ y = (lat_rot - RData.rot_lat_ll)/(RData.delta_rot_lat); x = lon_rot/(RData.delta_rot_lon); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise RotatedLatLonGrid::xy_to_latlon(x, y, lat, lon); @@ -310,6 +311,8 @@ const double range_max_deg = deg_per_km*Range_max_km; RotatedLatLonGrid::latlon_to_xy(lat, lon, x, y); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise + azi_deg = x*(RData.delta_rot_lon); range_deg = range_max_deg - y*(RData.delta_rot_lat); @@ -324,54 +327,23 @@ return; //////////////////////////////////////////////////////////////////////// -void TcrmwGrid::wind_ne_to_ra (const double lat, const double lon, - const double east_component, const double north_component, - double & radial_component, double & azimuthal_component) const +void TcrmwGrid::wind_ne_to_rt (const double azi_deg, + const double u_wind, const double v_wind, + double & radial_wind, double & tangential_wind) const { -Vector E, N, V; -Vector B_range, B_azi; -double azi_deg, range_deg, range_km; - - -latlon_to_range_azi(lat, lon, range_km, azi_deg); - -range_deg = deg_per_km*range_km; - -E = latlon_to_east (lat, lon); -N = latlon_to_north (lat, lon); - -V = east_component*E + north_component*N; - - -range_azi_to_basis(range_deg, azi_deg, B_range, B_azi); - - - - radial_component = dot(V, B_range); - -azimuthal_component = dot(V, B_azi); - - - - - -return; +double rcos = cosd(azi_deg); +double rsin = sind(azi_deg); +if (is_bad_data(u_wind) || is_bad_data(v_wind)) { + radial_wind = bad_data_double; + tangential_wind = bad_data_double; +} +else { + radial_wind = rcos*u_wind + rsin*v_wind; + tangential_wind = -1.0*rsin*u_wind + rcos*v_wind; } - - -//////////////////////////////////////////////////////////////////////// - - -void TcrmwGrid::wind_ne_to_ra_conventional (const double lat, const double lon, - const double east_component, const double north_component, - double & radial_component, double & azimuthal_component) const - -{ - -wind_ne_to_ra(lat, lon, east_component, north_component, radial_component, azimuthal_component); return; @@ -381,33 +353,17 @@ return; //////////////////////////////////////////////////////////////////////// -void TcrmwGrid::range_azi_to_basis(const double range_deg, const double azi_deg, Vector & B_range, Vector & B_azi) const +void TcrmwGrid::wind_ne_to_rt (const double lat, const double lon, + const double u_wind, const double v_wind, + double & radial_wind, double & tangential_wind) const { -double u, v, w; - - -u = cosd(range_deg)*sind(azi_deg); - -v = cosd(range_deg)*cosd(azi_deg); - -w = -sind(range_deg); - - - -B_range = u*Ir + v*Jr + w*Kr; - +double range_km, azi_deg; -u = cosd(azi_deg); - -v = -sind(azi_deg); - -w = 0.0; - - -B_azi = u*Ir + v*Jr + w*Kr; +latlon_to_range_azi(lat, lon, range_km, azi_deg); +wind_ne_to_rt(azi_deg, u_wind, v_wind, radial_wind, tangential_wind); return; @@ -425,8 +381,9 @@ RotatedLatLonGrid::latlon_to_xy(true_lat, true_lon, x, y); x -= Nx*floor(x/Nx); -x -= Nx*floor(x/Nx); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise +y -= Ny*floor(y/Ny); return; @@ -442,7 +399,9 @@ void TcrmwGrid::xy_to_latlon(double x, double y, double & true_lat, double & tru x -= Nx*floor(x/Nx); -x -= Nx*floor(x/Nx); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise + +y -= Ny*floor(y/Ny); RotatedLatLonGrid::xy_to_latlon(x, y, true_lat, true_lon); @@ -500,7 +459,3 @@ return; //////////////////////////////////////////////////////////////////////// - - - - diff --git a/src/libcode/vx_grid/tcrmw_grid.h b/src/libcode/vx_grid/tcrmw_grid.h index c102332408..3757111e9e 100644 --- a/src/libcode/vx_grid/tcrmw_grid.h +++ b/src/libcode/vx_grid/tcrmw_grid.h @@ -35,11 +35,8 @@ class TcrmwGrid : public RotatedLatLonGrid { void calc_ijk(); // calculate rotated basis vectors - void range_azi_to_basis(const double range_deg, const double azi_deg, Vector & B_range, Vector & B_azi) const; - TcrmwData TData; - Vector Ir, Jr, Kr; int Range_n, Azimuth_n; // # of points in the radial and azimuthal directions @@ -89,23 +86,15 @@ class TcrmwGrid : public RotatedLatLonGrid { void xy_to_latlon(double x, double y, double & true_lat, double & true_lon) const; + void wind_ne_to_rt(const double azi_deg, + const double u_wind, const double v_wind, + double & radial_wind, double & tangential_wind) const; + void wind_ne_to_rt(const double lat, const double lon, + const double u_wind, const double v_wind, + double & radial_wind, double & tangential_wind) const; - void wind_ne_to_ra(const double lat, const double lon, - const double east_component, const double north_component, - double & radial_component, double & azimuthal_component) const; - - - // - // possibly toggles the signs of the radial and/or azimuthal components - // - // to align with the conventions used in the TC community - // - void wind_ne_to_ra_conventional (const double lat, const double lon, - const double east_component, const double north_component, - double & radial_component, double & azimuthal_component) const; - }; diff --git a/src/libcode/vx_pointdata_python/python_pointdata.cc b/src/libcode/vx_pointdata_python/python_pointdata.cc index cf6b2b2418..cb41a05bf8 100644 --- a/src/libcode/vx_pointdata_python/python_pointdata.cc +++ b/src/libcode/vx_pointdata_python/python_pointdata.cc @@ -351,8 +351,10 @@ bool process_point_data_list(PyObject *python_point_data, MetPointDataPython &me // get valid time index vld_time = obs.getValidTime(); if ( !header_data->vld_num_array.has(vld_time, vld_idx) ) { + // MET #2897 keep vld_array and vld_num_array in sync + header_data->vld_array.add(obs.getValidTimeString()); header_data->vld_num_array.add(vld_time); - header_data->vld_num_array.has(vld_time, vld_idx); + vld_idx = header_data->vld_num_array.n() - 1; } if (!is_eq(prev_lat, lat) || !is_eq(prev_lon, lon) || !is_eq(prev_elv, elv) @@ -364,7 +366,6 @@ bool process_point_data_list(PyObject *python_point_data, MetPointDataPython &me header_data->sid_idx_array.add(sid); header_data->typ_idx_array.add(typ_idx); header_data->vld_idx_array.add(vld_idx); - header_data->vld_array.add(obs.getValidTimeString()); prev_lat = lat; prev_lon = lon; @@ -765,7 +766,7 @@ void print_met_data(MetPointObsData *obs_data, MetPointHeader *header_data, << header_data->vld_idx_array.n() << ", lat=" << header_data->lat_array.n() << ", lon=" << header_data->lon_array.n() << ", elv=" - << header_data->elv_array.n() << ", message_type=" + << header_data->elv_array.n() << ", message_type=" << header_data->typ_array.n() << ", station_id=" << header_data->sid_array.n() << ", valid_time=" << header_data->vld_array.n() << ", prpt=" @@ -776,7 +777,7 @@ void print_met_data(MetPointObsData *obs_data, MetPointHeader *header_data, log_count = (header_data->hdr_count > min_count) ? min_count : header_data->hdr_count; mlog << Debug(debug_level) << method_name - << "header_data: message_type,station_id,time_time,lat,lon.elv\n"; + << "header_data: message_type,station_id,time_time,lat,lon,elv\n"; for (int idx=0; idx " << "Could not find data for " << data_info->magic_time_str() - << " in file list:\n:" << write_css(search_files) << "\n\n"; + << " in file list:\n" << write_css(search_files) << "\n\n"; exit(1); } diff --git a/src/libcode/vx_stat_out/stat_columns.cc b/src/libcode/vx_stat_out/stat_columns.cc index 1fd22b47e6..87a30118c5 100644 --- a/src/libcode/vx_stat_out/stat_columns.cc +++ b/src/libcode/vx_stat_out/stat_columns.cc @@ -3335,11 +3335,13 @@ void write_pjc_cols(const PCTInfo &pct_info, col++; at.set_entry(r, col, // OY_TP - pct_info.pct.event_count_by_row(i)/(double) n); + (n == 0 ? bad_data_double : + pct_info.pct.event_count_by_row(i)/(double) n)); col++; at.set_entry(r, col, // ON_TP - pct_info.pct.nonevent_count_by_row(i)/(double) n); + (n == 0 ? bad_data_double : + pct_info.pct.nonevent_count_by_row(i)/(double) n)); col++; at.set_entry(r, col, // CALIBRATION @@ -3430,7 +3432,7 @@ void write_eclv_cols(const TTContingencyTable &ct, // // Economic Cost/Loss Value // Dump out the ECLV line: - // TOTAL, BASER, BASER_VALUE, + // TOTAL, BASER, VALUE_BASER, // N_PNT, [CL_], [VALUE_] (for each point) // at.set_entry(r, c+0, // Total Number of pairs diff --git a/src/libcode/vx_statistics/contable_nx2.cc b/src/libcode/vx_statistics/contable_nx2.cc index 2c2be22d62..55b270511e 100644 --- a/src/libcode/vx_statistics/contable_nx2.cc +++ b/src/libcode/vx_statistics/contable_nx2.cc @@ -468,7 +468,12 @@ return ( x ); double Nx2ContingencyTable::baser() const { - return ( (double) event_col_total()/n() ); + double v; + + if( n() == 0 ) v = bad_data_double; + else v = (double) event_col_total()/n(); + + return ( v ); } @@ -476,7 +481,7 @@ double Nx2ContingencyTable::baser() const { double Nx2ContingencyTable::baser_ci(double alpha, - double &cl, double &cu) const { + double &cl, double &cu) const { double v; v = baser(); diff --git a/src/libcode/vx_statistics/contable_stats.cc b/src/libcode/vx_statistics/contable_stats.cc index 9df41040a0..e9ed4feaf6 100644 --- a/src/libcode/vx_statistics/contable_stats.cc +++ b/src/libcode/vx_statistics/contable_stats.cc @@ -737,6 +737,8 @@ double TTContingencyTable::sedi_ci(double alpha, double TTContingencyTable::cost_loss(double r) const { double num, den, h, m, f, b, v; + if(n() == 0) return(bad_data_double); + // Total proportion of hits, misses, false alarms, and observations h = (double) fy_oy() / n(); m = (double) fn_oy() / n(); diff --git a/src/libcode/vx_tc_util/vx_tc_nc_util.cc b/src/libcode/vx_tc_util/vx_tc_nc_util.cc index 4c5297d8a1..9278fd5118 100644 --- a/src/libcode/vx_tc_util/vx_tc_nc_util.cc +++ b/src/libcode/vx_tc_util/vx_tc_nc_util.cc @@ -158,6 +158,21 @@ void write_tc_rmw(NcFile* nc_out, //////////////////////////////////////////////////////////////////////// +bool has_pressure_level(vector levels) { + + bool status = false; + + for (int j = 0; j < levels.size(); j++) { + if (levels[j].substr(0, 1) == "P") { + status = true; + break; + } + } + + return status; +} + +//////////////////////////////////////////////////////////////////////// set get_pressure_level_strings( map > variable_levels) { @@ -323,7 +338,7 @@ void def_tc_range_azimuth(NcFile* nc_out, add_att(&range_var, "_FillValue", bad_data_double); add_att(&azimuth_var, "long_name", "azimuth"); - add_att(&azimuth_var, "units", "degrees_clockwise_from_north"); + add_att(&azimuth_var, "units", "degrees_clockwise_from_east"); add_att(&azimuth_var, "standard_name", "azimuth"); add_att(&azimuth_var, "_FillValue", bad_data_double); @@ -538,7 +553,7 @@ void def_tc_variables(NcFile* nc_out, string long_name = variable_long_names[i->first]; string units = variable_units[i->first]; - if (levels.size() > 1) { + if (has_pressure_level(levels)) { data_var = nc_out->addVar( var_name, ncDouble, dims_3d); add_att(&data_var, "long_name", long_name); diff --git a/src/libcode/vx_tc_util/vx_tc_nc_util.h b/src/libcode/vx_tc_util/vx_tc_nc_util.h index 981d55e08e..e0239b531f 100644 --- a/src/libcode/vx_tc_util/vx_tc_nc_util.h +++ b/src/libcode/vx_tc_util/vx_tc_nc_util.h @@ -35,6 +35,8 @@ extern void write_tc_track_point(netCDF::NcFile*, extern void write_tc_rmw(netCDF::NcFile*, const netCDF::NcDim&, const TrackInfo&); +extern bool has_pressure_level(std::vector); + extern std::set get_pressure_level_strings( std::map >); diff --git a/src/tools/core/ensemble_stat/ensemble_stat.cc b/src/tools/core/ensemble_stat/ensemble_stat.cc index 62885a5f5e..98c92a0fc5 100644 --- a/src/tools/core/ensemble_stat/ensemble_stat.cc +++ b/src/tools/core/ensemble_stat/ensemble_stat.cc @@ -70,6 +70,7 @@ // 038 09/06/22 Halley Gotway MET #1908 Remove ensemble processing logic. // 039 09/29/22 Halley Gotway MET #2286 Refine GRIB1 table lookup logic. // 040 10/03/22 Prestopnik MET #2227 Remove using namespace netCDF from header files +// 041 06/17/24 Halley Gotway MET #2856 Reinitialize climo_cdf pointer // //////////////////////////////////////////////////////////////////////// @@ -2221,6 +2222,7 @@ void do_pct_cat_thresh(const EnsembleStatVxOpt &vx_opt, // Re-initialize pd_pnt.erase(); + pd_pnt.set_climo_cdf_info_ptr(&vx_opt.cdf_info); // Process the observations for(i_obs=0; i_obs, - map, map &); +static int combine_tqz_and_uv(map, map, + vector &); static float compute_pbl(map pqtzuv_map_tq, map pqtzuv_map_uv); static void copy_pqtzuv(float *to_pqtzuv, float *from_pqtzuv, bool copy_all=true); @@ -2958,14 +2958,16 @@ void copy_pqtzuv(float *to_pqtzuv, float *from_pqtzuv, bool copy_all) { //////////////////////////////////////////////////////////////////////// int combine_tqz_and_uv(map pqtzuv_map_tq, - map pqtzuv_map_uv, map &pqtzuv_map_merged) { + map pqtzuv_map_uv, vector &pqtzuv_merged_array) { static const char *method_name = "combine_tqz_and_uv() "; int tq_count = pqtzuv_map_tq.size(); int uv_count = pqtzuv_map_uv.size(); + map pqtzuv_map_merged; + pqtzuv_merged_array.clear(); if (tq_count > 0 && uv_count > 0) { IntArray common_levels, tq_levels; float *pqtzuv_tq, *pqtzuv_uv; - float *pqtzuv_merged = (float *) 0; + float *pqtzuv_merged = (float *) nullptr; float *next_pqtzuv, *prev_pqtzuv; float tq_pres_max, tq_pres_min, uv_pres_max, uv_pres_min; map::iterator it, it_tq, it_uv; @@ -3054,9 +3056,10 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, } interpolate_pqtzuv(prev_pqtzuv, pqtzuv_merged, next_pqtzuv); } - float first_pres = (pqtzuv_merged[0] == 0 ? bad_data_float : pqtzuv_merged[0]); + float first_pres = (pqtzuv_merged[0] < 0 || is_eq(pqtzuv_merged[0], 0.) + ? bad_data_float : pqtzuv_merged[0]); pqtzuv_map_merged[first_pres] = pqtzuv_merged; - mlog << Debug(9) << method_name << "Added " << first_pres << " to merged records\n"; + mlog << Debug(9) << method_name << "Added " << first_pres << " to merged records (first record)\n"; if (pqtzuv_merged != 0) { //Merge UV into TQZ records @@ -3064,6 +3067,12 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, //Merge TQZ into UV records merge_records(pqtzuv_merged, pqtzuv_map_uv, pqtzuv_map_tq, pqtzuv_map_merged); } + for (map::iterator it=pqtzuv_map_merged.begin(); + it!=pqtzuv_map_merged.end(); ++it) { + float *new_pqtzuv = new float[mxr8vt]; + for (int i=0; isecond[i]; + pqtzuv_merged_array.push_back(new_pqtzuv); + } if(mlog.verbosity_level() >= PBL_DEBUG_LEVEL) { log_merged_tqz_uv(pqtzuv_map_tq, pqtzuv_map_uv, pqtzuv_map_merged, method_name); @@ -3071,7 +3080,7 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, delete [] pqtzuv_merged; } - return pqtzuv_map_merged.size(); + return pqtzuv_merged_array.size(); } //////////////////////////////////////////////////////////////////////// @@ -3092,12 +3101,13 @@ float compute_pbl(map pqtzuv_map_tq, mlog << Debug(7) << method_name << "is called: TQZ: " << tq_count << " UV: " << uv_count << "\n"; if (tq_count > 0 || uv_count > 0) { + float *pqtzuv = nullptr; int hgt_cnt, spfh_cnt; IntArray selected_levels; - map pqtzuv_map_merged; + vector pqtzuv_merged_array; pbl_level = combine_tqz_and_uv(pqtzuv_map_tq, pqtzuv_map_uv, - pqtzuv_map_merged); + pqtzuv_merged_array); mlog << Debug(7) << method_name << "pbl_level= " << pbl_level << " from TQ (" << tq_count << ") and UV (" << uv_count << ")\n"; @@ -3111,60 +3121,62 @@ float compute_pbl(map pqtzuv_map_tq, << "Skip CALPBL because of only one available record after combining TQZ and UV\n"; } else { - // Order all observations by pressure from bottom to top - index = pbl_level - 1; + // Reverse the order all observations by pressure from bottom to top + index = 0; hgt_cnt = spfh_cnt = 0; - for (it=pqtzuv_map_merged.begin(); it!=pqtzuv_map_merged.end(); ++it) { - if (index < 0) { - mlog << Error << "\n" << method_name << "negative index: " << index << "\n\n"; - break; - } - - if (index < MAX_PBL_LEVEL) { - float *pqtzuv = it->second; - pbl_data_pres[index] = pqtzuv[0]; - pbl_data_spfh[index] = pqtzuv[1]; - pbl_data_temp[index] = pqtzuv[2]; - pbl_data_hgt[index] = pqtzuv[3]; - pbl_data_ugrd[index] = pqtzuv[4]; - pbl_data_vgrd[index] = pqtzuv[5]; - if (is_valid_pb_data(pbl_data_spfh[index])) spfh_cnt++; - if (is_valid_pb_data(pbl_data_hgt[index])) hgt_cnt++; - selected_levels.add(nint(it->first)); - } - - index--; - } - if (index != -1) { - mlog << Error << "\n" << method_name << "Missing some levels (" << index << ")\n"; - } - - if (pbl_level > MAX_PBL_LEVEL) { - it = pqtzuv_map_tq.begin(); + int start_offset = (MAX_PBL_LEVEL >= pbl_level) ? 0 : (pbl_level-MAX_PBL_LEVEL); + for (int i=(pbl_level-1); i>=start_offset; i--,index++) { + pqtzuv = pqtzuv_merged_array[i]; + pbl_data_pres[index] = pqtzuv[0]; + pbl_data_pres[index] = pqtzuv[0]; + pbl_data_spfh[index] = pqtzuv[1]; + pbl_data_temp[index] = pqtzuv[2]; + pbl_data_hgt[index] = pqtzuv[3]; + pbl_data_ugrd[index] = pqtzuv[4]; + pbl_data_vgrd[index] = pqtzuv[5]; + if (is_valid_pb_data(pbl_data_spfh[index])) spfh_cnt++; + if (is_valid_pb_data(pbl_data_hgt[index])) hgt_cnt++; + selected_levels.add(nint(pqtzuv[0])); + } + if (start_offset > 0) { + // Replace the interpolated records with common records. + mlog << Error << "\n" << method_name << "Excluded " << start_offset << " records\n"; // Find vertical levels with both data float highest_pressure = bad_data_float; - for (; it!=pqtzuv_map_tq.end(); ++it) { + for (it = pqtzuv_map_tq.begin(); it!=pqtzuv_map_tq.end(); ++it) { if (pqtzuv_map_uv.count(it->first) > 0) { highest_pressure = it->first; break; } } if (!is_bad_data(highest_pressure)) { + bool found; + int vector_idx = start_offset - 1; index = MAX_PBL_LEVEL - 1; for (; it!=pqtzuv_map_tq.end(); ++it) { int pres_level = nint(it->first); + // Stop replacing if already exists at input list if (selected_levels.has(pres_level)) break; - float *pqtzuv = pqtzuv_map_merged[it->first]; - pbl_data_pres[index] = pqtzuv[0]; - pbl_data_spfh[index] = pqtzuv[1]; - pbl_data_temp[index] = pqtzuv[2]; - pbl_data_hgt[index] = pqtzuv[3]; - pbl_data_ugrd[index] = pqtzuv[4]; - pbl_data_vgrd[index] = pqtzuv[5]; - mlog << Debug(6) << method_name << "Force to add " - << pres_level << " into " << index << "\n"; - index--; + found = false; + for (; vector_idx>=0; vector_idx--) { + if (is_eq(pqtzuv_merged_array[vector_idx][0], it->first)) { + pqtzuv = pqtzuv_merged_array[vector_idx]; + pbl_data_pres[index] = pqtzuv[0]; + pbl_data_spfh[index] = pqtzuv[1]; + pbl_data_temp[index] = pqtzuv[2]; + pbl_data_hgt[index] = pqtzuv[3]; + pbl_data_ugrd[index] = pqtzuv[4]; + pbl_data_vgrd[index] = pqtzuv[5]; + mlog << Debug(6) << method_name << "Force to add " + << pres_level << " into " << index << "\n"; + vector_idx--; + found = true; + break; + } + } + if (vector_idx < 0) break; + if(found) index--; } } } @@ -3204,9 +3216,11 @@ float compute_pbl(map pqtzuv_map_tq, if (!is_valid_pb_data(hpbl)) mlog << Debug(5) << method_name << " fail to compute PBL. TQ records: " << tq_count << " UV records: " << uv_count << " merged records: " - << pqtzuv_map_merged.size() << "\n"; + << pqtzuv_merged_array.size() << "\n"; } } + for (int i; i= log_debug_level) { if (qc_min_value > qc_value) qc_min_value = qc_value; @@ -2626,8 +2628,6 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, } 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; @@ -2637,15 +2637,19 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, } } + bool filter_out = is_eq(qc_for_flag, bad_data_int); 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++; + if (2 == qc_value) { + if (2 > qc_for_flag) { + if (0 == qc_for_flag) cnt_adp_qc_high_to_low++; + else if (1 == qc_for_flag) cnt_adp_qc_medium_to_low++; + qc_for_flag = 2; /* high/medium to low quality */ + } } - else if (2 == qc_value) { - qc_for_flag = 2; /* high/medium to low quality */ - cnt_adp_qc_adjused++; + else if (1 == qc_value && 0 == qc_for_flag) { + qc_for_flag = 1; /* high to medium quality */ + cnt_adp_qc_high_to_medium++; } if (has_qc_flags && !qc_flags.has(qc_for_flag)) filter_out = true; } @@ -2654,7 +2658,7 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, continue; } } - else if (has_qc_var && has_qc_flags && !qc_flags.has(qc_value)) { + else if (has_qc_var && !qc_flags.has(qc_value)) { qc_filtered_count++; continue; } @@ -2694,6 +2698,10 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, delete [] from_data; delete [] adp_qc_data; + int cnt_adjused_low = cnt_adp_qc_low + cnt_adp_qc_high_to_low + cnt_adp_qc_medium_to_low; + int cnt_adjused_high = cnt_adp_qc_high - cnt_adp_qc_high_to_medium - cnt_adp_qc_high_to_low; + int cnt_adjused_medium = cnt_adp_qc_medium + cnt_adp_qc_high_to_medium - cnt_adp_qc_medium_to_low; + int cnt_adjused_total = cnt_adp_qc_high_to_medium + cnt_adp_qc_high_to_low + cnt_adp_qc_medium_to_low; mlog << Debug(log_debug_level) << method_name << "Count: actual: " << to_cell_count << ", missing: " << missing_count << ", non_missing: " << non_missing_count << "\n Filtered: by QC: " << qc_filtered_count @@ -2701,13 +2709,22 @@ void regrid_goes_variable(NcFile *nc_in, VarInfo *vinfo, << ", by absent: " << absent_count << ", total: " << (qc_filtered_count + adp_qc_filtered_count + absent_count) << "\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"; + << "] QC: [" << qc_min_value << " - " << qc_max_value << "]\n"; + if (has_qc_flags) + mlog << Debug(log_debug_level) + << " 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_adjused_high << " (" << cnt_adp_qc_high + << "), medium=" << cnt_adjused_medium << " (" << cnt_adp_qc_medium + << "), low=" << cnt_adjused_low << " (" << cnt_adp_qc_low + << "), no_retrieval=" << cnt_adp_qc_nr + << "\n adjusted: high to medium=" << cnt_adp_qc_high_to_medium + << ", high to low=" << cnt_adp_qc_high_to_low + << ", medium to low=" << cnt_adp_qc_medium_to_low + << ", total=" << cnt_adjused_total + << "\n"; + if (to_cell_count == 0) { mlog << Warning << "\n" << method_name @@ -2922,7 +2939,6 @@ void set_adp_gc_values(NcVar var_adp_qc) { 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; } diff --git a/src/tools/tc_utils/tc_diag/tc_diag.cc b/src/tools/tc_utils/tc_diag/tc_diag.cc index c808faa29a..11bb4c04e8 100644 --- a/src/tools/tc_utils/tc_diag/tc_diag.cc +++ b/src/tools/tc_utils/tc_diag/tc_diag.cc @@ -847,7 +847,7 @@ void compute_lat_lon(TcrmwGrid& grid, ia * grid.azimuth_delta_deg(), lat, lon); lat_arr[i] = lat; - lon_arr[i] = -lon; // degrees east to west + lon_arr[i] = -lon; // degrees west to east } } diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw.cc b/src/tools/tc_utils/tc_rmw/tc_rmw.cc index 5479136ed0..5c69775609 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw.cc @@ -67,22 +67,21 @@ static bool file_is_ok(const ConcatString &, const GrdFileType); static void process_rmw(); static void process_tracks(TrackInfoArray&); static void get_atcf_files(const StringArray&, - const StringArray&, StringArray&, StringArray&); + const StringArray&, StringArray&, StringArray&); static void process_track_files(const StringArray&, - const StringArray&, TrackInfoArray&); + const StringArray&, TrackInfoArray&); static bool is_keeper(const ATCFLineBase *); static void set_deck(const StringArray&); static void set_atcf_source(const StringArray&, - StringArray&, StringArray&); + StringArray&, StringArray&); static void set_data_files(const StringArray&); static void set_config(const StringArray&); static void set_out(const StringArray&); static void setup_grid(); static void setup_nc_file(); static void build_outfile_name(const ConcatString&, - const char*, ConcatString&); -static void compute_lat_lon(TcrmwGrid&, - double*, double*); + const char*, ConcatString&); +static void compute_lat_lon(TcrmwGrid&, double*, double*); static void process_fields(const TrackInfoArray&); //////////////////////////////////////////////////////////////////////// @@ -550,19 +549,19 @@ void set_out(const StringArray& a) { void setup_grid() { - grid_data.name = "TCRMW"; - grid_data.range_n = conf_info.n_range; - grid_data.azimuth_n = conf_info.n_azimuth; + tcrmw_data.name = "TCRMW"; + tcrmw_data.range_n = conf_info.n_range; + tcrmw_data.azimuth_n = conf_info.n_azimuth; // Define the maximum range in km based on the fixed increment if(is_bad_data(conf_info.rmw_scale)) { - grid_data.range_max_km = + tcrmw_data.range_max_km = conf_info.delta_range_km * (conf_info.n_range - 1); } - tcrmw_grid.set_from_data(grid_data); - grid.set(grid_data); + tcrmw_grid.set_from_data(tcrmw_data); + grid_out.set(tcrmw_data); } //////////////////////////////////////////////////////////////////////// @@ -612,11 +611,11 @@ void setup_nc_file() { // Get VarInfo data_info = conf_info.data_info[i_var]; mlog << Debug(4) << "Processing field: " << data_info->magic_str() << "\n"; - string fname = data_info->name_attr(); + string fname = data_info->name_attr(); variable_levels[fname].push_back(data_info->level_attr()); variable_long_names[fname] = data_info->long_name_attr(); variable_units[fname] = data_info->units_attr(); - wind_converter.update_input(fname, data_info->units_attr()); + wind_converter.update_input(fname, data_info->units_attr()); } // Define pressure levels @@ -649,8 +648,8 @@ void compute_lat_lon(TcrmwGrid& tcrmw_grid, ir * tcrmw_grid.range_delta_km(), ia * tcrmw_grid.azimuth_delta_deg(), lat, lon); - lat_arr[i] = lat; - lon_arr[i] = - lon; + lat_arr[i] = lat; + lon_arr[i] = -lon; // degrees west to east } } } @@ -691,12 +690,12 @@ void process_fields(const TrackInfoArray& tracks) { << point.lon() << ").\n"; // Set grid center - grid_data.lat_center = point.lat(); - grid_data.lon_center = -1.0*point.lon(); // internal sign change + tcrmw_data.lat_center = point.lat(); + tcrmw_data.lon_center = -1.0*point.lon(); // internal sign change // Define the maximum range in km relative to the radius of maximum winds if(!is_bad_data(conf_info.rmw_scale)) { - grid_data.range_max_km = + tcrmw_data.range_max_km = conf_info.rmw_scale * point.mrd() * tc_km_per_nautical_miles * (conf_info.n_range - 1); @@ -704,9 +703,9 @@ void process_fields(const TrackInfoArray& tracks) { // Re-define the range/azimuth grid tcrmw_grid.clear(); - tcrmw_grid.set_from_data(grid_data); - grid.clear(); - grid.set(grid_data); + tcrmw_grid.set_from_data(tcrmw_data); + grid_out.clear(); + grid_out.set(tcrmw_data); // Compute lat and lon coordinate arrays compute_lat_lon(tcrmw_grid, lat_arr, lon_arr); @@ -734,35 +733,33 @@ void process_fields(const TrackInfoArray& tracks) { data_info->set_valid(valid_time); // Find data for this track point - get_series_entry(i_point, data_info, data_files, ftype, data_dp, latlon_arr); - - // Check data range - double data_min, data_max; - data_dp.data_range(data_min, data_max); - mlog << Debug(4) << "data_min:" << data_min << "\n"; - mlog << Debug(4) << "data_max:" << data_max << "\n"; - - // Regrid data - data_dp = met_regrid(data_dp, latlon_arr, grid, data_info->regrid()); - data_dp.data_range(data_min, data_max); - mlog << Debug(4) << "data_min:" << data_min << "\n"; - mlog << Debug(4) << "data_max:" << data_max << "\n"; - - // if this is "U", setup everything for matching "V" and compute the radial/tangential - if(wind_converter.compute_winds_if_input_is_u(i_point, sname, slevel, valid_time, data_files, ftype, - latlon_arr, lat_arr, lon_arr, grid, data_dp, tcrmw_grid)) { + get_series_entry(i_point, data_info, data_files, ftype, data_dp, grid_in); + + // Regrid data and log the range of values before and after + double dmin, dmax, dmin_rgd, dmax_rgd; + data_dp.data_range(dmin, dmax); + data_dp = met_regrid(data_dp, grid_in, grid_out, data_info->regrid()); + data_dp.data_range(dmin_rgd, dmax_rgd); + + mlog << Debug(4) << data_info->magic_str() + << " input range (" << dmin << ", " << dmax + << "), regrid range (" << dmin_rgd << ", " << dmax_rgd << ")\n"; + + // if this is "U", setup everything for matching "V" and compute the radial/tangential + if(wind_converter.compute_winds_if_input_is_u(i_point, sname, slevel, valid_time, data_files, ftype, + grid_in, grid_out, data_dp, tcrmw_grid)) { write_tc_pressure_level_data(nc_out, tcrmw_grid, pressure_level_indices, data_info->level_attr(), i_point, - data_3d_vars[conf_info.radial_velocity_field_name.string()], - wind_converter.get_wind_r_arr()); - write_tc_pressure_level_data(nc_out, tcrmw_grid, - pressure_level_indices, data_info->level_attr(), i_point, - data_3d_vars[conf_info.tangential_velocity_field_name.string()], - wind_converter.get_wind_t_arr()); + data_3d_vars[conf_info.radial_velocity_field_name.string()], + wind_converter.get_wind_r_arr()); + write_tc_pressure_level_data(nc_out, tcrmw_grid, + pressure_level_indices, data_info->level_attr(), i_point, + data_3d_vars[conf_info.tangential_velocity_field_name.string()], + wind_converter.get_wind_t_arr()); } - + // Write data - if(variable_levels[data_info->name_attr()].size() > 1) { + if(has_pressure_level(variable_levels[data_info->name_attr()])) { write_tc_pressure_level_data(nc_out, tcrmw_grid, pressure_level_indices, data_info->level_attr(), i_point, data_3d_vars[data_info->name_attr()], data_dp.data()); diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw.h b/src/tools/tc_utils/tc_rmw/tc_rmw.h index 87e37b74be..85ffef1c64 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw.h +++ b/src/tools/tc_utils/tc_rmw/tc_rmw.h @@ -84,7 +84,6 @@ static TCRMWConfInfo conf_info; static GrdFileType ftype; static TCRMW_WindConverter wind_converter; - // Optional arguments static ConcatString out_dir; static ConcatString out_prefix; @@ -136,20 +135,15 @@ static std::map pressure_level_indices; // //////////////////////////////////////////////////////////////////////// -static DataPlane dp; -static Grid latlon_arr; -static TcrmwData grid_data; +static Grid grid_in; +static TcrmwData tcrmw_data; static TcrmwGrid tcrmw_grid; -static Grid grid; +static Grid grid_out; // Grid coordinate arrays static double* lat_arr; static double* lon_arr; -// Wind arrays -/* static double* wind_r_arr; */ -/* static double* wind_t_arr; */ - //////////////////////////////////////////////////////////////////////// #endif // __TC_RMW_H__ diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc index f94994b314..76b54037c7 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc @@ -17,6 +17,7 @@ // ---- ---- ---- ----------- // 000 05/11/22 Albo Pulled the wind conversion into a class // 001 09/28/22 Prestopnik MET #2227 Remove namespace std from header files +// 002 05/03/24 Halley Gotway MET #2841 Fix radial and tangential winds // //////////////////////////////////////////////////////////////////////// @@ -28,14 +29,13 @@ using namespace std; //////////////////////////////////////////////////////////////////////// -static void wind_ne_to_ra(const TcrmwGrid&, - const DataPlane&, const DataPlane&, - const double*, const double*, double*, double*); +static void wind_ne_to_rt(const TcrmwGrid&, + const DataPlane&, const DataPlane&, + double*, double*); //////////////////////////////////////////////////////////////////////// -void TCRMW_WindConverter::_free_winds_arrays(void) -{ +void TCRMW_WindConverter::_free_winds_arrays(void) { if (_windR != nullptr) { delete [] _windR; _windR = nullptr; @@ -94,6 +94,7 @@ void TCRMW_WindConverter::init(const TCRMWConfInfo *conf) { _vIndexMap[varlevel] = i_var; } } + // test for consistency if (_uIndexMap.size() != _vIndexMap.size()) { mlog << Warning << "Uneven number of u/v wind inputs, no wind conversion will be done:\n" @@ -101,12 +102,16 @@ void TCRMW_WindConverter::init(const TCRMWConfInfo *conf) { << _conf->v_wind_field_name.string() << " has " << _vIndexMap.size() << " inputs\n"; _computeWinds = false; } - map::const_iterator iu, iv; - for (iu=_uIndexMap.begin(), iv=_vIndexMap.begin(); iu!=_uIndexMap.end(); ++iu, ++iv) { - if (iu->first != iv->first) { - mlog << Warning << "Ordering of u/v wind input levels not the same, not implemented, no wind conversions will be done:\n" - << " " << iu->first << " " << iv->first << "\n"; - _computeWinds = false; + if (_computeWinds) { + map::const_iterator iu, iv; + for (iu=_uIndexMap.begin(), iv=_vIndexMap.begin(); iu!=_uIndexMap.end(); ++iu, ++iv) { + if (iu->first != iv->first) { + mlog << Warning << "Ordering of u/v wind input levels not the same, " + << "not implemented, no wind conversions will be done:\n" + << " " << iu->first << " " << iv->first << "\n"; + _computeWinds = false; + break; + } } } } @@ -128,11 +133,9 @@ void TCRMW_WindConverter::update_input(const string &variableName, const string //////////////////////////////////////////////////////////////////////// void TCRMW_WindConverter::append_nc_output_vars(map > &variable_levels, - map &variable_long_names, - map &variable_units) { - if (!_computeWinds) { - return; - } + map &variable_long_names, + map &variable_units) { + if (!_computeWinds) return; if (_foundUInInput && _foundVInInput) { variable_levels[_conf->tangential_velocity_field_name] = variable_levels[_conf->u_wind_field_name.string()]; @@ -144,14 +147,15 @@ void TCRMW_WindConverter::append_nc_output_vars(map > &va } else { if (!_foundUInInput) { - mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " - << "field not found in input \"" << _conf->u_wind_field_name << "\"\n\n"; + mlog << Warning << "\nTCWRMW_WindConverter::append_nc_output_vars() -> " + << "field not found in input \"" << _conf->u_wind_field_name << "\"\n\n"; } if (!_foundVInInput) { - mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " - << "field not found in input \"" << _conf->v_wind_field_name << "\"\n\n"; + mlog << Warning << "\nTCWRMW_WindConverter::append_nc_output_vars() -> " + << "field not found in input \"" << _conf->v_wind_field_name << "\"\n\n"; } - mlog << Warning << "\nNot computing radial and tangential winds\n\n"; + mlog << Warning << "\nTCWRMW_WindConverter::append_nc_output_vars() -> " + << "Not computing radial and tangential winds\n\n"; _computeWinds = false; } } @@ -159,84 +163,90 @@ void TCRMW_WindConverter::append_nc_output_vars(map > &va //////////////////////////////////////////////////////////////////////// bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point, - const string &varName, - const string &varLevel, - unixtime valid_time, - const StringArray &data_files, - const GrdFileType &ftype, - const Grid &latlon_arr, - const double *lat_arr, - const double *lon_arr, - const Grid &grid, - const DataPlane &data_dp, - const TcrmwGrid &tcrmw_grid) { + const string &varName, + const string &varLevel, + unixtime valid_time, + const StringArray &data_files, + const GrdFileType &ftype, + const Grid &grid_in, + const Grid &grid_out, + const DataPlane &u_wind_dp, + const TcrmwGrid &tcrmw_grid) { if (!_computeWinds) { return false; } int uIndex = -1; int vIndex = -1; - VarInfo *data_infoV = (VarInfo *) 0; + VarInfo *v_wind_info = (VarInfo *) 0; if (varName == _conf->u_wind_field_name.string()) { uIndex = _uIndexMap[varLevel]; vIndex = _vIndexMap[varLevel]; - data_infoV = _conf->data_info[vIndex]; - data_infoV->set_valid(valid_time); + v_wind_info = _conf->data_info[vIndex]; + v_wind_info->set_valid(valid_time); } else { // not the U input return false; } - DataPlane data_dpV; - Grid latlon_arrV; - get_series_entry(i_point, data_infoV, data_files, ftype, data_dpV, - latlon_arrV); - double data_min, data_max; - data_dpV.data_range(data_min, data_max); - mlog << Debug(4) << "V data_min:" << data_min << "\n"; - mlog << Debug(4) << "V data_max:" << data_max << "\n"; - data_dpV = met_regrid(data_dpV, latlon_arr, grid, data_infoV->regrid()); - data_dpV.data_range(data_min, data_max); - mlog << Debug(4) << "V data_min:" << data_min << "\n"; - mlog << Debug(4) << "V data_max:" << data_max << "\n"; - - // here's the conversion, at last - wind_ne_to_ra(tcrmw_grid, data_dp, data_dpV, lat_arr, lon_arr, - _windR, _windT); - // _windR and _windT now set + DataPlane v_wind_dp; + Grid v_wind_grid; + get_series_entry(i_point, v_wind_info, data_files, ftype, + v_wind_dp, v_wind_grid); + double dmin, dmax, dmin_rgd, dmax_rgd; + v_wind_dp.data_range(dmin, dmax); + v_wind_dp = met_regrid(v_wind_dp, v_wind_grid, grid_out, v_wind_info->regrid()); + v_wind_dp.data_range(dmin_rgd, dmax_rgd); + + mlog << Debug(4) << v_wind_info->magic_str() + << " input range (" << dmin << ", " << dmax + << "), regrid range (" << dmin_rgd << ", " << dmax_rgd << ")\n"; + + // Compute the radial and tangential winds and store in _windR and _windT + wind_ne_to_rt(tcrmw_grid, u_wind_dp, v_wind_dp, _windR, _windT); + return true; } - //////////////////////////////////////////////////////////////////////// -void wind_ne_to_ra(const TcrmwGrid& tcrmw_grid, - const DataPlane& u_dp, const DataPlane& v_dp, - const double* lat_arr, const double* lon_arr, - double* wind_r_arr, double* wind_t_arr) { - // Transform (u, v) to (radial, azimuthal) - for(int ir = 0; ir < tcrmw_grid.range_n(); ir++) { - for(int ia = 0; ia < tcrmw_grid.azimuth_n(); ia++) { - int i = ir * tcrmw_grid.azimuth_n() + ia; - double lat = lat_arr[i]; - double lon = - lon_arr[i]; - double u = u_dp.data()[i]; - double v = v_dp.data()[i]; - double wind_r; - double wind_t; - if(is_bad_data(u) || is_bad_data(v)) { - mlog << Debug(4) << "wind_ne_to_ra: latlon:" << lat << "," << lon << " winds are missing\n"; - wind_r = bad_data_double; - wind_t = bad_data_double; - } else { - tcrmw_grid.wind_ne_to_ra(lat, lon, u, v, wind_r, wind_t); - mlog << Debug(4) << "wind_ne_to_ra: latlon:" << lat << "," << lon << " uv:" << u << "," - << v << ", rt:" << wind_r << "," << wind_t <<"\n"; - } - wind_r_arr[i] = wind_r; - wind_t_arr[i] = wind_t; - } - } +void wind_ne_to_rt(const TcrmwGrid& tcrmw_grid, + const DataPlane& u_dp, const DataPlane& v_dp, + double* wind_r_arr, double* wind_t_arr) { + + int n_rng = tcrmw_grid.range_n(); + int n_azi = tcrmw_grid.azimuth_n(); + + // Transform (u, v) to (radial, tangential) winds + for(int ir = 0; ir < n_rng; ir++) { + for(int ia = 0; ia < n_azi; ia++) { + + // Store data in reverse order + int i_rev = (n_rng - ir - 1) * n_azi + ia; + + double azi_deg = ia * tcrmw_grid.azimuth_delta_deg(); + double range_km = ir * tcrmw_grid.range_delta_km(); + + double lat, lon; + tcrmw_grid.range_azi_to_latlon(range_km, azi_deg, lat, lon); + + tcrmw_grid.wind_ne_to_rt(azi_deg, u_dp.data()[i_rev], v_dp.data()[i_rev], + wind_r_arr[i_rev], wind_t_arr[i_rev]); + + mlog << Debug(4) << "wind_ne_to_rt() -> " + << "center lat/lon (" << tcrmw_grid.lat_center_deg() + << ", " << tcrmw_grid.lon_center_deg() + << "), range (km): " << range_km + << ", azimuth (deg): " << azi_deg + << ", point lat/lon (" << lat << ", " << lon + << "), uv (" << u_dp.data()[i_rev] << ", " << v_dp.data()[i_rev] + << "), radial wind: " << wind_r_arr[i_rev] + << ", tangential wind: " << wind_t_arr[i_rev] << "\n"; + } // end for ia + } // end for ir + + return; } +//////////////////////////////////////////////////////////////////////// diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h index 8ad0642a17..d80b8ed15d 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h @@ -15,7 +15,7 @@ // // Mod# Date Name Description // ---- ---- ---- ----------- -// 000 05/11/22 DaveAlbo New +// 000 05/11/22 Albo New // //////////////////////////////////////////////////////////////////////// @@ -32,7 +32,6 @@ using std::map; using std::string; - //////////////////////////////////////////////////////////////////////// // // Constants @@ -98,8 +97,8 @@ class TCRMW_WindConverter { // if configured to compute winds, but didn't find U or V, turn off // the wind computations and report an error void append_nc_output_vars(std::map > &variable_levels, - std::map &variable_long_names, - std::map &variable_units); + std::map &variable_long_names, + std::map &variable_units); // Check input varName against U, and if it's a match, lookup V using the // map members, and then compute tangential and radial winds if it is so @@ -109,17 +108,15 @@ class TCRMW_WindConverter { // If true if returned, the winds can be accessed by calls to // get_wind_t_arr() and get_wind_r_arr() bool compute_winds_if_input_is_u(int i_point, - const string &varName, - const string &varLevel, - unixtime valid_time, - const StringArray &data_files, - const GrdFileType &ftype, - const Grid &latlon_arr, - const double *lat_arr, - const double *lon_arr, - const Grid &grid, - const DataPlane &data_dp, - const TcrmwGrid &tcrmw_grid); + const string &varName, + const string &varLevel, + unixtime valid_time, + const StringArray &data_files, + const GrdFileType &ftype, + const Grid &grid_in, + const Grid &grid_out, + const DataPlane &u_wind_dp, + const TcrmwGrid &tcrmw_grid); };