From 1dcae35cbf3b1d813fe71e7824652cdf7ebd4890 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Wed, 6 Apr 2022 16:57:18 -0600 Subject: [PATCH] Per #2118, porting changes over from bugfix_2118_main_v10.1_grib1_rotll to fix this issue in the develop branch. --- met/src/libcode/vx_data2d_grib/data2d_grib.cc | 4 + .../libcode/vx_data2d_grib/grib_classes.cc | 21 ++++ met/src/libcode/vx_data2d_grib/grib_classes.h | 30 ++++- met/src/libcode/vx_data2d_grib/grib_utils.cc | 113 +++++++++++++++--- .../libcode/vx_data2d_grib2/data2d_grib2.cc | 2 +- .../libcode/vx_data2d_nc_met/get_met_grid.cc | 2 +- 6 files changed, 154 insertions(+), 18 deletions(-) diff --git a/met/src/libcode/vx_data2d_grib/data2d_grib.cc b/met/src/libcode/vx_data2d_grib/data2d_grib.cc index ef5e622b9c..c92e64569a 100644 --- a/met/src/libcode/vx_data2d_grib/data2d_grib.cc +++ b/met/src/libcode/vx_data2d_grib/data2d_grib.cc @@ -762,6 +762,10 @@ bool is_grid_relative(const GribRecord &r) { else if(r.gds->type == 5) { res_flag = r.gds->grid_type.stereographic.res_flag; } + // Rotated LatLon + else if(r.gds->type == 10) { + res_flag = r.gds->grid_type.rot_latlon_grid.res_flag; + } else { mlog << Error << "\nis_grid_relative() -> " << "Unsupported grid type value: " << r.gds->type diff --git a/met/src/libcode/vx_data2d_grib/grib_classes.cc b/met/src/libcode/vx_data2d_grib/grib_classes.cc index da1a5f0f76..d02445c6ef 100644 --- a/met/src/libcode/vx_data2d_grib/grib_classes.cc +++ b/met/src/libcode/vx_data2d_grib/grib_classes.cc @@ -1874,6 +1874,27 @@ else if ((h.type == 5)) file << " scan_flag: " << (int) h.grid_type.stereographic.scan_flag << "\n\n"; } +else if ((h.type == 10)) +{ + + file << " lat1: " << char3_to_int(h.grid_type.rot_latlon_grid.lat1) << "\n"; + file << " lon1: " << char3_to_int(h.grid_type.rot_latlon_grid.lon1) << "\n"; + + file << " res_flag: " << (int) h.grid_type.rot_latlon_grid.res_flag << "\n"; + + file << " lat2: " << char3_to_int(h.grid_type.rot_latlon_grid.lat2) << "\n"; + file << " lon2: " << char3_to_int(h.grid_type.rot_latlon_grid.lon2) << "\n"; + + file << " di: " << char2_to_int(h.grid_type.rot_latlon_grid.di) << "\n"; + file << " dj: " << char2_to_int(h.grid_type.rot_latlon_grid.dj) << "\n"; + + file << " scan_flag: " << (int) h.grid_type.rot_latlon_grid.scan_flag << "\n"; + + file << " lat_sp: " << char3_to_int(h.grid_type.rot_latlon_grid.lat_sp) << "\n"; + file << " lon_sp: " << char3_to_int(h.grid_type.rot_latlon_grid.lon_sp) << "\n"; + + file << " rotation: " << char4_to_dbl(h.grid_type.rot_latlon_grid.rotation) << "\n\n"; +} return ( file ); diff --git a/met/src/libcode/vx_data2d_grib/grib_classes.h b/met/src/libcode/vx_data2d_grib/grib_classes.h index 55b7046d65..3b0295cf45 100644 --- a/met/src/libcode/vx_data2d_grib/grib_classes.h +++ b/met/src/libcode/vx_data2d_grib/grib_classes.h @@ -162,7 +162,7 @@ struct Section1_Header { // PDS //////////////////////////////////////////////////////////////////////// -struct LatLon { // Latitude/Longitude Grid +struct LatLon { // Latitude/Longitude Grid unsigned char lat1[3]; // 11 - 13 unsigned char lon1[3]; // 14 - 16 @@ -179,6 +179,32 @@ struct LatLon { // Latitude/Longitude Grid unsigned char unused[14]; // 29 - 42 +}; + + // Reference: https://apps.ecmwf.int/codes/grib/format/grib1/grids/10 + +struct RotLatLon { // Rotated Latitude/Longitude Grid + + unsigned char lat1[3]; // 11 - 13 + unsigned char lon1[3]; // 14 - 16 + + unsigned char res_flag; // 17 + + unsigned char lat2[3]; // 18 - 20 + unsigned char lon2[3]; // 21 - 23 + + unsigned char di[2]; // 24 - 25 + unsigned char dj[2]; // 26 - 27 + + unsigned char scan_flag; // 28 + + unsigned char unused[4]; // 29 - 32 + + unsigned char lat_sp[3]; // 33 - 35 + unsigned char lon_sp[3]; // 36 - 38 + + unsigned char rotation[4]; // 39 - 42 + }; struct Mercator { // Mercator Grid @@ -282,7 +308,7 @@ struct Gaussian { union GridType { struct LatLon latlon_grid; // Latitude/Longitude Grid - // struct RotLatLon rot_latlon_grid; // Rotated Latitude/Longitude Grid + struct RotLatLon rot_latlon_grid; // Rotated Latitude/Longitude Grid struct Mercator mercator; // Mercator Grid struct LambertConf lambert_conf; // Lambert Conformal Secant Grid struct Stereographic stereographic; // Stereographic Grid diff --git a/met/src/libcode/vx_data2d_grib/grib_utils.cc b/met/src/libcode/vx_data2d_grib/grib_utils.cc index 6dc5b7c563..c1ef70f775 100644 --- a/met/src/libcode/vx_data2d_grib/grib_utils.cc +++ b/met/src/libcode/vx_data2d_grib/grib_utils.cc @@ -29,11 +29,13 @@ using namespace std; // grid types from the GDS section // -static const int latlon_type = 0; -static const int mercator_type = 1; -static const int lambert_type = 3; -static const int stereographic_type = 5; -static const int gaussian_type = 4; +static const int latlon_type = 0; +static const int mercator_type = 1; +static const int lambert_type = 3; +static const int gaussian_type = 4; +static const int stereographic_type = 5; +static const int rotated_latlon_type = 10; + //////////////////////////////////////////////////////////////////////// @@ -62,7 +64,7 @@ void gds_to_grid(const Section2_Header & gds, Grid & gr) LambertData lc_data; StereographicData st_data; LatLonData ll_data; -//RotatedLatLonData rll_data; +RotatedLatLonData rll_data; MercatorData mc_data; GaussianData g_data; @@ -84,11 +86,11 @@ if ( gds.type == latlon_type ) { gr.set(ll_data); -// else if ( gds.type == rotated_latlon_type ) { -// -// gds_to_rotated_latlon(gds, rll_data); -// -// gr.set(rll_data); +} else if ( gds.type == rotated_latlon_type ) { + + gds_to_rotated_latlon(gds, rll_data); + + gr.set(rll_data); } else if ( gds.type == mercator_type ) { @@ -143,6 +145,7 @@ void gds_to_order(const Section2_Header & gds, int & xdir, int & ydir, int & ord // Check GDS for the grid type. // The following Projection types are supported: // - Lat/Lon + // - Rotated Lat/Lon // - Mercator // - Lambert Conformal // - Polar Stereographic @@ -154,6 +157,10 @@ if ( gds.type == latlon_type ) { scan_flag_to_order(gds.grid_type.latlon_grid.scan_flag, xdir, ydir, order); +} else if ( gds.type == rotated_latlon_type ) { + + scan_flag_to_order(gds.grid_type.rot_latlon_grid.scan_flag, xdir, ydir, order); + } else if (gds.type == mercator_type ) { scan_flag_to_order(gds.grid_type.mercator.scan_flag, xdir, ydir, order); @@ -278,10 +285,88 @@ void gds_to_rotated_latlon(const Section2_Header & gds, RotatedLatLonData & data { -mlog << Error << "\ngds_to_rotated_latlon() -> " - << "Rotated Lat/Lon grids are not supported in GRIB version 1.\n\n"; +int xdir, ydir, order; +double d; + +scan_flag_to_order(gds.grid_type.rot_latlon_grid.scan_flag, xdir, ydir, order); + + // Store the grid name +data.name = rotated_latlon_proj_type; + + // + // Multiply longitude values by -1 since the NCAR code considers + // degrees west to be positive. + // + + // Latitude of the bottom left corner +data.rot_lat_ll = min(decode_lat_lon(gds.grid_type.rot_latlon_grid.lat1, 3), + decode_lat_lon(gds.grid_type.rot_latlon_grid.lat2, 3)); + + // Longitude of the bottom left corner +if ( xdir == 0 ) data.rot_lon_ll = -1.0*rescale_lon(decode_lat_lon(gds.grid_type.rot_latlon_grid.lon1, 3)); +else data.rot_lon_ll = -1.0*rescale_lon(decode_lat_lon(gds.grid_type.rot_latlon_grid.lon2, 3)); + + // Number of points in the Latitudinal (y) direction +data.Nlat = char2_to_int(gds.ny); + + // Number of points in the Longitudinal (x) direction +data.Nlon = char2_to_int(gds.nx); + + // Check for thinned lat/lon grids +if ( data.Nlon == 65535 ) { -exit ( 1 ); + mlog << Error << "\ngds_to_rotated_latlon() -> " + << "Thinned Lat/Lon grids are not supported for GRIB version 1.\n\n"; + + exit ( 1 ); + +} + + // Compute latitudinal increment from lat1 and lat2 +d = fabs(decode_lat_lon(gds.grid_type.rot_latlon_grid.lat1, 3) + - decode_lat_lon(gds.grid_type.rot_latlon_grid.lat2, 3)); + +data.delta_rot_lat = d/(data.Nlat); + + // If explicitly specified and non-zero, use it +if ( ! all_bits_set(gds.grid_type.rot_latlon_grid.dj, 2) ) { + + d = decode_lat_lon(gds.grid_type.rot_latlon_grid.dj, 2); + + if ( d > 0 ) data.delta_rot_lat = d; + +} + + // Compute longitudinal increment from lon1 and lon2 +d = fabs(decode_lat_lon(gds.grid_type.rot_latlon_grid.lon1, 3) + - decode_lat_lon(gds.grid_type.rot_latlon_grid.lon2, 3)); + +data.delta_rot_lon = d/(data.Nlon); + + // If explicitly specified and non-zero, use it +if ( ! all_bits_set(gds.grid_type.rot_latlon_grid.di, 2) ) { + + d = decode_lat_lon(gds.grid_type.rot_latlon_grid.di, 2); + + if ( d > 0 ) data.delta_rot_lon = d; + +} + + // Location of (rotated) south pole +data.true_lat_south_pole = + decode_lat_lon(gds.grid_type.rot_latlon_grid.lat_sp, 3); + +data.true_lon_south_pole = + -1.0*rescale_lon(decode_lat_lon(gds.grid_type.rot_latlon_grid.lon_sp, 3)); + + // Auxiliary rotation +data.aux_rotation = char4_to_dbl(gds.grid_type.rot_latlon_grid.rotation); + +data.dump(); + + // + // done + // return; diff --git a/met/src/libcode/vx_data2d_grib2/data2d_grib2.cc b/met/src/libcode/vx_data2d_grib2/data2d_grib2.cc index 0061d2fd00..65df4c7ef7 100644 --- a/met/src/libcode/vx_data2d_grib2/data2d_grib2.cc +++ b/met/src/libcode/vx_data2d_grib2/data2d_grib2.cc @@ -1100,7 +1100,7 @@ void MetGrib2DataFile::read_grib2_grid( gribfield *gfld) { data.true_lon_south_pole = s_lon; // - // auxilliary rotation around the rotated polar axis + // auxiliary rotation around the rotated polar axis // data.aux_rotation = (t[21])*angle_factor; diff --git a/met/src/libcode/vx_data2d_nc_met/get_met_grid.cc b/met/src/libcode/vx_data2d_nc_met/get_met_grid.cc index 8e9e68cf4d..3512c6aaed 100644 --- a/met/src/libcode/vx_data2d_nc_met/get_met_grid.cc +++ b/met/src/libcode/vx_data2d_nc_met/get_met_grid.cc @@ -309,7 +309,7 @@ data.rot_lon_ll *= -1.0; get_global_att(ncfile, string("true_lon_south_pole"), data.true_lon_south_pole); if ( !west_longitude_positive ) data.true_lon_south_pole *= -1.0; - // auxilliary rotation + // auxiliary rotation get_global_att(ncfile, string("aux_rotation"), data.aux_rotation);