From 4f81f59e37ac2ed42694ded0fb97ba8f93699fa5 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Fri, 20 Jan 2023 16:08:37 -0700 Subject: [PATCH 01/11] #2228 Added parameters for ellipsoidal earth --- src/libcode/vx_grid/st_grid_defs.h | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/libcode/vx_grid/st_grid_defs.h b/src/libcode/vx_grid/st_grid_defs.h index 66596de81a..ab44bdcb10 100644 --- a/src/libcode/vx_grid/st_grid_defs.h +++ b/src/libcode/vx_grid/st_grid_defs.h @@ -25,9 +25,9 @@ struct StereographicData { const char * name; - char hemisphere; // 'N' or 'S' + char hemisphere; // 'N' or 'S' from latitude_of_projection_origin - double scale_lat; + double scale_lat; //latitude_of_projection_origin double lat_pin; double lon_pin; @@ -35,11 +35,17 @@ struct StereographicData { double x_pin; double y_pin; - double lon_orient; + double lon_orient; // -straight_vertical_longitude_from_pole double d_km; - double r_km; + double r_km; // semi_major_axis + + // ellipsoidal earth + double eccentricity; // 0 for shperical earth + double false_east; + double false_north; + double scale_factor; int nx; int ny; From d06a6f543c5f36a94c724f1302a25636a8fea67c Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Fri, 20 Jan 2023 16:11:32 -0700 Subject: [PATCH 02/11] #2228 Added functions for ellipsoidal earth --- src/libcode/vx_grid/st_grid.cc | 156 +++++++++++++++++++++++++++++---- src/libcode/vx_grid/st_grid.h | 14 ++- 2 files changed, 152 insertions(+), 18 deletions(-) diff --git a/src/libcode/vx_grid/st_grid.cc b/src/libcode/vx_grid/st_grid.cc index b0b79bd12d..4e71f1dbdc 100644 --- a/src/libcode/vx_grid/st_grid.cc +++ b/src/libcode/vx_grid/st_grid.cc @@ -1,5 +1,3 @@ - - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2022 // ** University Corporation for Atmospheric Research (UCAR) @@ -108,6 +106,7 @@ Ny = data.ny; Name = data.name; Lon_orient = data.lon_orient; +Data = data; // @@ -116,6 +115,7 @@ Lon_orient = data.lon_orient; Alpha = (1.0 + H*sind(data.scale_lat))*((data.r_km)/(data.d_km)); + // // Calculate Bx, By // @@ -129,8 +129,6 @@ theta0 = H*(Lon_orient - data.lon_pin); Bx = data.x_pin - Alpha*r0*H*sind(theta0); By = data.y_pin + Alpha*r0*H*cosd(theta0); -Data = data; - // // Done // @@ -203,13 +201,23 @@ const double H = ( IsNorthHemisphere ? 1.0 : -1.0 ); reduce(lon); -r = st_func(lat, is_north()); +if(is_eq(Data.eccentricity, 0.0)) { -theta = H*(Lon_orient - lon); + r = st_func(lat, is_north()); -x = Bx + Alpha*r*H*sind(theta); + theta = H*(Lon_orient - lon); -y = By - Alpha*r*H*cosd(theta); + x = Bx + Alpha*r*H*sind(theta); + + y = By - Alpha*r*H*cosd(theta); +} +else { + st_latlon_to_xy_func(lat, lon, x, y, Data.scale_factor, (Data.lon_orient*-1.0), + (Data.r_km*m_per_km), Data.false_east, Data.false_north, + Data.eccentricity, IsNorthHemisphere); + x = x / (Data.d_km*m_per_km) - Data.x_pin; + y = y / (Data.d_km*m_per_km) - Data.y_pin; +} return; @@ -227,19 +235,27 @@ double r, theta; const double H = ( IsNorthHemisphere ? 1.0 : -1.0 ); -x = (x - Bx)/Alpha; -y = (y - By)/Alpha; +if(is_eq(Data.eccentricity, 0.0)) { + x = (x - Bx)/Alpha; + y = (y - By)/Alpha; + + r = sqrt( x*x + y*y ); -r = sqrt( x*x + y*y ); + lat = st_inv_func(r, is_north()); -lat = st_inv_func(r, is_north()); + if ( fabs(r) < 1.0e-5 ) theta = 0.0; + else theta = atan2d(H*x, -H*y); // NOT atan2d(y, x); -if ( fabs(r) < 1.0e-5 ) theta = 0.0; -else theta = atan2d(H*x, -H*y); // NOT atan2d(y, x); + if ( is_south() ) theta = -theta; -if ( is_south() ) theta = -theta; + lon = Lon_orient - theta; +} +else { + st_xy_to_latlon_func(x, y, lat, lon, Data.scale_factor, (Data.r_km*m_per_km), + (-1.0*Data.lon_orient), Data.false_east, Data.false_north, + Data.eccentricity, IsNorthHemisphere); +} -lon = Lon_orient - theta; reduce(lon); @@ -608,7 +624,7 @@ return ( p ); //////////////////////////////////////////////////////////////////////// -double st_func(double lat, bool is_north_hemisphere) +double st_func(double lat, bool is_north_hemisphere, double eccentricity) { @@ -617,6 +633,9 @@ double r; if ( is_north_hemisphere ) r = tand(45.0 - 0.5*lat); else r = tand(45.0 + 0.5*lat); +if (!is_eq(eccentricity, 0.)) { + r *= pow(((1 + eccentricity*sind(lat)) / (1 - eccentricity*sind(lat))),(eccentricity/2)); +} return ( r ); } @@ -660,6 +679,109 @@ return ( a ); //////////////////////////////////////////////////////////////////////// +bool st_latlon_to_xy_func(double lat, double lon, double &x_m, double &y_m, + double scale_factor, double scale_lat, double semi_major_axis, + double false_east, double false_north, + double e, bool is_north_hemisphere) +{ + + const double lat_rad = lat * rad_per_deg; + const double lon_rad = lon * rad_per_deg; + const double lat_sin = sin(lat_rad); + const double lonO_rad = scale_lat * rad_per_deg; + const double H = (is_north_hemisphere? 1.0 : -1.0 ); + double t = tan(M_PI/4 - H*lat_rad/2) * pow((1 + e*lat_sin)/(1 - e*lat_sin), e/2); + double rho = (2 * semi_major_axis * scale_factor * t) + / sqrt(pow(1+e,1+e) * pow(1-e,1-e)); + // meters in polar stereographics, not index + x_m = false_east + rho*sin(lon_rad - lonO_rad); + y_m = false_north - H*rho*cos(lon_rad - lonO_rad); + + return true; +} + +//////////////////////////////////////////////////////////////////////// + + +bool st_xy_to_latlon_func(double x_m, double y_m, double &lat, double &lon, + double scale_factor, double semi_major_axis, + double proj_vertical_lon, double false_east, double false_north, + double eccentricity, bool is_north_hemisphere) +{ +bool result = true; +double chi; +double x_diff = x_m - false_east; +double y_diff = y_m - false_north; +double lonO_rad = proj_vertical_lon * rad_per_deg; +double r_rho = sqrt(x_diff*x_diff + y_diff*y_diff); +double t = r_rho * sqrt(pow((1+eccentricity),(1+eccentricity)) * pow((1-eccentricity),(1-eccentricity))) + / (2*semi_major_axis*scale_factor); +if (is_north_hemisphere) chi = M_PI/2 - 2 * atan(t); +else chi = 2 * atan(t) - M_PI/2; + +lat = chi + (eccentricity*eccentricity/2 + 5*pow(eccentricity,4)/24 + + pow(eccentricity,6)/12 + 13*pow(eccentricity,8)/360*sin(2 * chi) + + (7*pow(eccentricity,4)/48 + 29*pow(eccentricity,6)/240 + 811*pow(eccentricity,8)/11520)*sin(4*chi) + + (7*pow(eccentricity,6)/120 + 81*pow(eccentricity,8)/1120)*sin(6*chi) + + (4279*pow(eccentricity,8)/161280)*sin(8*chi)); + +if (x_m == false_east) lon = lonO_rad; +else if (is_north_hemisphere) lon = lonO_rad + atan2(x_diff,-y_diff); +else lon = lonO_rad + atan2(x_diff,y_diff); + +lat /= rad_per_deg; +lon /= rad_per_deg; +return result; + +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + +double st_eccentricity_func(double semi_major_axis, double semi_minor_axis, + double inverse_flattening) +{ + +double c; +double eccentricity = bad_data_double; + +if (is_eq(semi_minor_axis, bad_data_double)) { + semi_minor_axis = semi_major_axis - semi_major_axis/inverse_flattening; +} +eccentricity = sqrt(semi_major_axis*semi_major_axis - semi_minor_axis*semi_minor_axis) / semi_major_axis; + +return ( eccentricity ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +double st_sf_func(double standard_parallel, double eccentricity, bool is_north_hemisphere) +{ + +double scale_factor; +double tF, mF, temp1, temp2; +double sp_rad = standard_parallel * rad_per_deg; +double lat_sin = sin(sp_rad); +double lat_cos = cos(sp_rad); + +temp1 = eccentricity * lat_sin; +temp2 = pow((1 + temp1)/(1 - temp1), eccentricity/2); +if (is_north_hemisphere) tF = tan(M_PI/4 - sp_rad/2) * temp2; +else tF = tan(M_PI/4 + sp_rad/2) / temp2; +mF = lat_cos / sqrt(1 - eccentricity*eccentricity * lat_sin*lat_sin); +scale_factor = mF * sqrt(pow((1+eccentricity),(1+eccentricity)) * pow((1-eccentricity),(1-eccentricity))) / (2 * tF); + +return ( scale_factor ); + +} + +//////////////////////////////////////////////////////////////////////// + + double stereographic_alpha(double scale_lat, double r_km, double d_km) { diff --git a/src/libcode/vx_grid/st_grid.h b/src/libcode/vx_grid/st_grid.h index 53d5150b9b..ff9a1e3997 100644 --- a/src/libcode/vx_grid/st_grid.h +++ b/src/libcode/vx_grid/st_grid.h @@ -121,11 +121,23 @@ inline double StereographicGrid::scale_km () const { return ( Data.d_km ); } extern double stereographic_alpha(double scale_lat, double r_km, double d_km); -extern double st_func (double lat, bool is_north_hemisphere); +extern double st_func (double lat, bool is_north_hemisphere, double eccentricity = 0.); extern double st_der_func (double lat, bool is_north_hemisphere); extern double st_inv_func (double r, bool is_north_hemisphere); +extern double st_eccentricity_func(double semi_major_axis, double semi_minor_axis, + double inverse_flattening); +extern double st_sf_func(double standard_parallel, double eccentricity, bool is_north_hemisphere); +extern bool st_latlon_to_xy_func(double lat, double lon, double &x_m, double &y_m, + double scale_factor, double scale_lat, double semi_major_axis, + double false_east, double false_north, + double e, bool is_north_hemisphere); +extern bool st_xy_to_latlon_func(double x_m, double y_m, double &lat, double &lon, + double scale_factor, double semi_major_axis, + double proj_vertical_lon, double false_east, double false_north, + double eccentricity, bool is_north_hemisphere); + //////////////////////////////////////////////////////////////////////// From 56f79f7b8e7637aec26a8cfc908b69e29142bbc4 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Fri, 20 Jan 2023 16:15:00 -0700 Subject: [PATCH 03/11] #2218 Support polar stereographics without scale_factor and ellipsoidal earth --- src/libcode/vx_data2d_nccf/nccf_file.cc | 282 +++++++++++++++++++----- 1 file changed, 232 insertions(+), 50 deletions(-) diff --git a/src/libcode/vx_data2d_nccf/nccf_file.cc b/src/libcode/vx_data2d_nccf/nccf_file.cc index 3e3cfcdbad..df3642d469 100644 --- a/src/libcode/vx_data2d_nccf/nccf_file.cc +++ b/src/libcode/vx_data2d_nccf/nccf_file.cc @@ -37,6 +37,7 @@ using namespace netCDF; static const int max_met_args = 30; static const double DELTA_TOLERANCE_PERCENT = 0.05; +static const double EARTH_MAJOR_AXIS_km = 6371.20; const double NcCfFile::DELTA_TOLERANCE = 15.0; @@ -53,7 +54,8 @@ static const char var_lon_nt[] = "tlon"; static ConcatString x_dim_var_name; static ConcatString y_dim_var_name; -static double get_nc_var_att_double(const NcVar *nc_var, const char *att_name); +static double get_nc_var_att_double(const NcVar *nc_var, const char *att_name, + bool is_required=true); #define USE_BUFFER 1 @@ -1648,7 +1650,9 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin << "assuming meters.\n\n"; } else { - if (x_coord_units_name == "m" ) x_coord_to_m_cf = 1.0; + if ( x_coord_units_name == "m" || + x_coord_units_name == "meter" || + x_coord_units_name == "meters") x_coord_to_m_cf = 1.0; else if (x_coord_units_name == "km") x_coord_to_m_cf = 1000.0; else { mlog << Error << "\n" << method_name << " -> " @@ -1671,7 +1675,9 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin << "assuming meters.\n\n"; } else { - if (y_coord_units_name == "m" ) y_coord_to_m_cf = 1.0; + if ( y_coord_units_name == "m" || + y_coord_units_name == "meter" || + y_coord_units_name == "meters" ) y_coord_to_m_cf = 1.0; else if (y_coord_units_name == "km") y_coord_to_m_cf = 1000.0; else { mlog << Error << "\n" << method_name << " -> " @@ -1705,8 +1711,10 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin double dx_m = (x_values[x_counts-1] - x_values[0]) / (x_counts - 1); double dy_m = (y_values[y_counts-1] - y_values[0]) / (y_counts - 1); + double dx_m_a = fabs(dx_m); + double dy_m_a = fabs(dy_m); - if (fabs(dx_m - dy_m) > DELTA_TOLERANCE) + if (fabs(dx_m_a - dy_m_a) > DELTA_TOLERANCE) { mlog << Error << "\n" << method_name << " -> " << "MET can only process Lambert Conformal files where the x-axis and y-axis deltas are the same\n\n"; @@ -1718,8 +1726,8 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin for (int i = 1; i < (int)x_counts; ++i) { - double curr_delta = x_values[i] - x_values[i-1]; - if (fabs(curr_delta - dx_m) > DELTA_TOLERANCE) + double curr_delta = fabs(x_values[i] - x_values[i-1]); + if (fabs(curr_delta - dx_m_a) > DELTA_TOLERANCE) { mlog << Error << "\n" << method_name << " -> " << "MET can only process Lambert Conformal files where the delta along the x-axis is constant\n\n"; @@ -1729,8 +1737,8 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin for (int i = 1; i < (int)y_counts; ++i) { - double curr_delta = y_values[i] - y_values[i-1]; - if (fabs(curr_delta - dy_m) > DELTA_TOLERANCE) + double curr_delta = fabs(y_values[i] - y_values[i-1]); + if (fabs(curr_delta - dy_m_a) > DELTA_TOLERANCE) { mlog << Error << "\n" << method_name << " -> " << "MET can only process Lambert Conformal files where the delta along the y-axis is constant\n\n"; @@ -1767,12 +1775,14 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin data.y_pin = y_pin; data.lon_orient = -double_datas[0]; data.d_km = dx_m / 1000.0; - data.r_km = 6371.20; + data.r_km = EARTH_MAJOR_AXIS_km; data.nx = _xDim->getSize(); data.ny = _yDim->getSize(); data.so2_angle = 0.0; grid.set(data); + if (dy_m < 0) grid.set_swap_to_north(true); + if(std_parallel_att) delete std_parallel_att; if(central_lon_att) delete central_lon_att; if(proj_origin_lat_att) delete proj_origin_lat_att; @@ -1937,16 +1947,19 @@ void NcCfFile::get_grid_mapping_orthographic(const NcVar *grid_mapping_var) //////////////////////////////////////////////////////////////////////// -double get_nc_var_att_double(const NcVar *nc_var, const char *att_name) +double get_nc_var_att_double(const NcVar *nc_var, const char *att_name, bool is_required) { NcVarAtt *nc_att = get_nc_att(nc_var, (string)att_name); if(IS_INVALID_NC_P(nc_att)) { - mlog << Error << "\nget_nc_var_att_double() -> " - << "Cannot get \"" << att_name << "\" from " - << GET_NC_NAME_P(nc_var) << " variable.\n\n"; - exit(1); + if (is_required) { + mlog << Error << "\nget_nc_var_att_double() -> " + << "Cannot get \"" << att_name << "\" from " + << GET_NC_NAME_P(nc_var) << " variable.\n\n"; + exit(1); + } + else return bad_data_double; } double att_val = get_att_value_double(nc_att); if (nc_att) delete nc_att; @@ -1960,30 +1973,48 @@ double get_nc_var_att_double(const NcVar *nc_var, const char *att_name) void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_var) { - static const string method_name = "NcCfFile::get_grid_mapping_polar_stereographic()"; double x_coord_to_m_cf = 1.0; double y_coord_to_m_cf = 1.0; + static const string method_name = "NcCfFile::get_grid_mapping_polar_stereographic() --> "; // Get projection attributes - - double proj_origin_lat = - get_nc_var_att_double(grid_mapping_var, - "latitude_of_projection_origin"); - double proj_origin_lon = - get_nc_var_att_double(grid_mapping_var, - "longitude_of_projection_origin"); - double proj_vertical_lon = - get_nc_var_att_double(grid_mapping_var, - "straight_vertical_longitude_from_pole"); - double proj_origin_scale_factor = - get_nc_var_att_double(grid_mapping_var, - "scale_factor_at_projection_origin"); + // proj_origin_lat: either 90.0 or -90.0, to decide the northern/southern hemisphere + bool is_spherical_earch = true; + double proj_origin_lat = get_nc_var_att_double( + grid_mapping_var, "latitude_of_projection_origin"); + double proj_vertical_lon = get_nc_var_att_double( + grid_mapping_var, "straight_vertical_longitude_from_pole"); + double proj_origin_lon = get_nc_var_att_double( + grid_mapping_var, "longitude_of_projection_origin", false); + double proj_standard_parallel = get_nc_var_att_double( + grid_mapping_var, "standard_parallel", false); + double proj_origin_scale_factor = get_nc_var_att_double( + grid_mapping_var, "scale_factor_at_projection_origin", false); + double semi_major_axis = get_nc_var_att_double( + grid_mapping_var, "semi_major_axis", false); + double semi_minor_axis = get_nc_var_att_double( + grid_mapping_var, "semi_minor_axis", false); + double inverse_flattening = get_nc_var_att_double( + grid_mapping_var, "inverse_flattening", false); + bool has_scale_factor = !is_eq(proj_origin_scale_factor, bad_data_double); + bool has_standard_parallel = !is_eq(proj_standard_parallel, bad_data_double); // Check that the scale factor at the origin is 1. - if(!is_eq(proj_origin_scale_factor, 1.0)) - { - mlog << Error << "\n" << method_name << " -> " + if (!is_eq(inverse_flattening, bad_data_double) || + (!is_eq(semi_minor_axis, bad_data_double) && !is_eq(semi_minor_axis,semi_major_axis))) { + is_spherical_earch = false; + mlog << Warning << "\n" << method_name + << "This is an ellipsoidal earth which is not fully supported.\n\n"; + } + else if(!has_scale_factor && !has_standard_parallel) { + mlog << Error << "\n" << method_name + << "The attribute \"scale_factor_at_projection_origin\" and \"standard_parallel\" of the " + << GET_NC_NAME_P(grid_mapping_var) << " variable do not exist.\n\n"; + exit(1); + } + else if(has_scale_factor && !is_eq(proj_origin_scale_factor, 1.0)) { + mlog << Error << "\n" << method_name << "Unexpected attribute value of " << proj_origin_scale_factor << " for the scale_factor_at_projection_origin attribute of the " << GET_NC_NAME_P(grid_mapping_var) << " variable.\n\n"; @@ -2042,21 +2073,21 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va if (_xDim == 0) { - mlog << Error << "\n" << method_name << " -> " + mlog << Error << "\n" << method_name << "Didn't find X dimension (projection_x_coordinate) in netCDF file.\n\n"; exit(1); } if (_yDim == 0) { - mlog << Error << "\n" << method_name << " -> " + mlog << Error << "\n" << method_name << "Didn't find Y dimension (projection_y_coordinate) in netCDF file.\n\n"; exit(1); } if (_xCoordVar == 0) { - mlog << Error << "\n" << method_name << " -> " + mlog << Error << "\n" << method_name << "Didn't find X coord variable (" << GET_NC_NAME_P(_xDim) << ") in netCDF file.\n\n"; exit(1); @@ -2064,7 +2095,7 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va if (_yCoordVar == 0) { - mlog << Error << "\n" << method_name << " -> " + mlog << Error << "\n" << method_name << "Didn't find Y coord variable (" << GET_NC_NAME_P(_yDim) << ") in netCDF file.\n\n"; exit(1); @@ -2073,7 +2104,7 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va if (get_data_size(_xCoordVar) != (int) GET_NC_SIZE_P(_xDim) || get_data_size(_yCoordVar) != (int) GET_NC_SIZE_P(_yDim)) { - mlog << Error << "\n" << method_name << " -> " + mlog << Error << "\n" << method_name << "Coordinate variables don't match dimension sizes in netCDF file.\n\n"; exit(1); } @@ -2085,21 +2116,24 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va ConcatString x_coord_units_name; const NcVarAtt *x_coord_units_att = get_nc_att(_xCoordVar, units_att_name); if (!get_var_units(_xCoordVar, x_coord_units_name)) { - mlog << Warning << "\n" << method_name << " -> " + mlog << Warning << "\n" << method_name << "Units not given for X coordinate variable -- assuming meters.\n\n"; } else { if (0 == x_coord_units_name.length()) { - mlog << Warning << "\n" << method_name << " -> " + mlog << Warning << "\n" << method_name << "Cannot extract X coordinate units from netCDF file -- " << "assuming meters.\n\n"; } else { - if ( x_coord_units_name == "m" ) x_coord_to_m_cf = 1.0; + if ( x_coord_units_name == "m" || + x_coord_units_name == "meter" || + x_coord_units_name == "meters") x_coord_to_m_cf = 1.0; else if ( x_coord_units_name == "km") x_coord_to_m_cf = 1000.0; else { - mlog << Error << "\n" << method_name << " -> " - << "The X coordinates must be in meters or kilometers for MET.\n\n"; + mlog << Error << "\n" << method_name + << "The X coordinates (" << x_coord_units_name + << ") must be in meters or kilometers for MET.\n\n"; exit(1); } } @@ -2108,20 +2142,22 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va ConcatString y_coord_units_name; const NcVarAtt *y_coord_units_att = get_nc_att(_yCoordVar, units_att_name); if (!get_var_units(_yCoordVar, y_coord_units_name)) { - mlog << Warning << "\n" << method_name << " -> " + mlog << Warning << "\n" << method_name << "Units not given for Y coordinate variable -- assuming meters.\n\n"; } else { if (0 == y_coord_units_name.length()) { - mlog << Warning << "\n" << method_name << " -> " + mlog << Warning << "\n" << method_name << "Cannot extract Y coordinate units from netCDF file -- " << "assuming meters.\n\n"; } else { - if ( y_coord_units_name == "m" ) y_coord_to_m_cf = 1.0; + if ( y_coord_units_name == "m" || + y_coord_units_name == "meter" || + y_coord_units_name == "meters" ) y_coord_to_m_cf = 1.0; else if ( y_coord_units_name == "km") y_coord_to_m_cf = 1000.0; else { - mlog << Error << "\n" << method_name << " -> " + mlog << Error << "\n" << method_name << "The X coordinates must be in meters or kilometers for MET.\n\n"; exit(1); } @@ -2150,14 +2186,20 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va double dx_m = (x_values[x_counts-1] - x_values[0]) / (x_counts - 1); double dy_m = (y_values[y_counts-1] - y_values[0]) / (y_counts - 1); + double dx_m_a = fabs(dx_m); + double dy_m_a = fabs(dy_m); - if (fabs(dx_m - dy_m) > DELTA_TOLERANCE) + if (fabs(dx_m_a - dy_m_a) > DELTA_TOLERANCE) { - mlog << Error << "\n" << method_name << " -> " + mlog << Error << "\n" << method_name << "MET can only process Polar Stereographic files where the x-axis and y-axis deltas are the same.\n\n"; exit(1); } + if (is_eq(semi_major_axis, bad_data_double)) + semi_major_axis = EARTH_MAJOR_AXIS_km * m_per_km; + else semi_major_axis *= x_coord_to_m_cf; // meters + // Calculate the pin indices. The pin will be located at the grid's reference // location since that's the only lat/lon location we know about. @@ -2169,17 +2211,157 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va StereographicData data; data.name = stereographic_proj_type; data.lat_pin = proj_origin_lat; - data.lon_pin = -1.0 * proj_origin_lon; - data.hemisphere = (data.lat_pin > 0 ? 'N' : 'S'); + data.lon_pin = -1.0 * (is_eq(proj_origin_lon, bad_data_double) ? 0 : proj_origin_lon); + data.hemisphere = (proj_origin_lat > 0 ? 'N' : 'S'); data.x_pin = x_pin; data.y_pin = y_pin; data.scale_lat = proj_origin_lat; data.lon_orient = -1.0 * proj_vertical_lon; data.d_km = dx_m / 1000.0; - data.r_km = 6371.20; + data.r_km = semi_major_axis / 1000.0; data.nx = _xDim->getSize(); data.ny = _yDim->getSize(); + + bool is_north_hemisphere = proj_origin_lat > 0; + double eccentricity, false_east,false_north, scale_factor; + scale_factor = proj_origin_scale_factor; + eccentricity = false_east = false_north = 0.; + if(!has_scale_factor && has_standard_parallel) { + double lat, lon; + double x, y, x2, y2; + + false_east = get_nc_var_att_double(grid_mapping_var, "false_east", false); + false_north = get_nc_var_att_double(grid_mapping_var, "false_north", false); + + if (is_eq(false_east, bad_data_double)) false_east = 0.; + if (is_eq(false_north, bad_data_double)) false_north = 0.; + if(!is_spherical_earch) eccentricity = st_eccentricity_func(semi_major_axis, semi_minor_axis, + inverse_flattening); + + x = (dx_m > 0) ? x_values[0] : x_values[x_counts-1]; + y = (dy_m > 0) ? y_values[0] : y_values[y_counts-1]; + scale_factor = st_sf_func(proj_standard_parallel, eccentricity, is_north_hemisphere); + st_xy_to_latlon_func(x, y, lat, lon, scale_factor, semi_major_axis, + proj_vertical_lon, false_east, false_north, + eccentricity, is_north_hemisphere); + if(is_eq(eccentricity, 0.0)) { + data.x_pin = 0.; + data.y_pin = 0.; + } + else { + data.x_pin = ((dx_m>0) ? x_values[0] : x_values[x_counts-1]) / dx_m_a; + data.y_pin = ((dy_m>0) ? y_values[0] : y_values[x_counts-1]) / dy_m_a; + } + data.lat_pin = lat; + data.lon_pin = -lon; + + st_latlon_to_xy_func(lat, lon, x2, y2, scale_factor, proj_vertical_lon, + semi_major_axis, false_east, false_north, + eccentricity, is_north_hemisphere); + mlog << Debug(6) << method_name + << "pin: (x,y)m -> (lat,lon) -> (x2,y2)m: (" << x << "," << y << ") -> (" + << lat << "," << lon << ") -> (" << x2 << "," << y2 + << ") Diff (m): x=" << (x-x2) <<", y=" << (y-y2) << "\n"; + + } + + // ellipsoidal earth + data.scale_factor = scale_factor; + data.eccentricity = eccentricity; + data.false_east = false_east; + data.false_north = false_north; + grid.set(data); + + if (dy_m < 0) grid.set_swap_to_north(true); + + if(mlog.verbosity_level() >= 10) { + double lat1, lon1; + double x1, y1, x2, y2; + + mlog << Debug(15) << method_name + << "dx_m=" << dx_m << ", dy_m=" << dy_m << "\n"; + + x1 = (dx_m > 0) ? x_values[0] : x_values[x_counts-1]; + y1 = (dy_m > 0) ? y_values[0] : y_values[y_counts-1]; + st_xy_to_latlon_func(x1, y1, lat1, lon1, scale_factor, semi_major_axis, + proj_vertical_lon, false_east, false_north, + eccentricity, is_north_hemisphere); + st_latlon_to_xy_func(lat1, lon1, x2, y2, scale_factor, proj_vertical_lon, + semi_major_axis, false_east, false_north, + eccentricity, is_north_hemisphere); + mlog << Debug(10) << method_name + << "left bottom: (x,y) -> (lat,lon) -> (x2,y2)m: (" + << x1 << "," << y1 << ") -> (" + << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 + << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) + << "\n"; + + x1 = x_values[int(x_counts/2)]; + y1 = y_values[int(y_counts/2)]; + st_xy_to_latlon_func(x1, y1, lat1, lon1, scale_factor, semi_major_axis, + proj_vertical_lon, false_east, false_north, + eccentricity, is_north_hemisphere); + st_latlon_to_xy_func(lat1, lon1, x2, y2, scale_factor, proj_vertical_lon, + semi_major_axis, false_east, false_north, + eccentricity, is_north_hemisphere); + mlog << Debug(10) << method_name + << " center: (x,y)m -> (lat,lon) -> (x2,y2)m: (" << x1 + << "," << y1 << ") -> (" + << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 + << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) + << "\n"; + + x1 = (dx_m > 0) ? x_values[x_counts-1] : x_values[0]; + y1 = (dy_m > 0) ? y_values[y_counts-1] : y_values[0]; + st_xy_to_latlon_func(x1, y1, lat1, lon1, scale_factor, semi_major_axis, + proj_vertical_lon, false_east, false_north, + eccentricity, is_north_hemisphere); + st_latlon_to_xy_func(lat1, lon1, x2, y2, scale_factor, proj_vertical_lon, + semi_major_axis, false_east, false_north, + eccentricity, is_north_hemisphere); + mlog << Debug(10) << method_name + << " top-right: (x,y)m -> (lat,lon) -> (x2,y2)m: (" << x1 + << "," << y1 << ") -> (" + << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 + << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) + << "\n"; + + if(mlog.verbosity_level() >= 15) { + mlog << Debug(15) << method_name + << "data.scale_factor=" << data.scale_factor + <<", -data.lon_orient=" << -1.0*data.lon_orient << ", data.r_km=" << data.r_km + << ", data.false_east=" << data.false_east << ", data.false_north=" << data.false_north + << ", data.eccentricity=" << data.eccentricity + << ", is_north_hemisphere=" << is_north_hemisphere << "\n"; + for (int ix=0; ix (lat,lon) -> (x2,y2)m: (" << x1 + << "," << y1 << ") -> (" + << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 + << ") Diff (m): x=" << x_diff <<", y=" << y_diff + << ((abs(x_diff/dx_m) > 0.5) ? " [x_delta > dx/2] " : " ") + << ((abs(y_diff/dy_m) > 0.5) ? " [y_delta > dy/2] " : " ") + << ((abs(x_diff) > dx_m_a || abs(y_diff) > dy_m_a) ? " *** check dx/dy ***" : " ") + << "\n"; + } + } + } + } + } From 29945fe49693c22506883e51271b8d2c409a06a0 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Fri, 20 Jan 2023 17:06:20 -0700 Subject: [PATCH 04/11] #2218 Addedmrtwo unit tests for polar stereographics --- .../test_unit/xml/unit_plot_data_plane.xml | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/internal/test_unit/xml/unit_plot_data_plane.xml b/internal/test_unit/xml/unit_plot_data_plane.xml index 3dc7c53f59..ba8203e15d 100644 --- a/internal/test_unit/xml/unit_plot_data_plane.xml +++ b/internal/test_unit/xml/unit_plot_data_plane.xml @@ -507,5 +507,33 @@ &OUTPUT_DIR;/plot_data_plane/nbm_2022021513_F119_GRIB2_FICEAC_A48_Perc-10.ps - + + + &MET_BIN;/plot_data_plane + \ + &DATA_DIR_MODEL;/nccf/radolan_sp.nc \ + &OUTPUT_DIR;/plot_data_plane/radolan_sp_PREC.ps \ + 'name="PREC"; level="(0,*,*)";' \ + -title "NCCF Polar Stereographic Precipitation" \ + -v 1 + + + &OUTPUT_DIR;/plot_data_plane/radolan_sp_PREC.ps + + + + + &MET_BIN;/plot_data_plane + \ + &DATA_DIR_MODEL;/nccf/NSIDC0081_SEAICE_PS_N25km_20220723_v2.0.nc \ + &OUTPUT_DIR;/plot_data_plane/NSIDC0081_SEAICE_PS_N25km.ps \ + 'name="F16_ICECON"; level="(0,*,*)";' \ + -title "Sea Ice Concentration" \ + -v 1 + + + &OUTPUT_DIR;/plot_data_plane/NSIDC0081_SEAICE_PS_N25km.ps + + + From dc7ec946515f02f0dc9c095ec2024cad27cb444f Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Mon, 23 Jan 2023 10:30:40 -0700 Subject: [PATCH 05/11] #2218 convert index to meter at xy_to_latlon --- src/libcode/vx_grid/st_grid.cc | 71 ++++++++++++++++++++-------------- 1 file changed, 43 insertions(+), 28 deletions(-) diff --git a/src/libcode/vx_grid/st_grid.cc b/src/libcode/vx_grid/st_grid.cc index 4e71f1dbdc..6b7db8ec7f 100644 --- a/src/libcode/vx_grid/st_grid.cc +++ b/src/libcode/vx_grid/st_grid.cc @@ -201,7 +201,7 @@ const double H = ( IsNorthHemisphere ? 1.0 : -1.0 ); reduce(lon); -if(is_eq(Data.eccentricity, 0.0)) { +if(is_eq(Data.eccentricity, 0.0) || Data.eccentricity >= 0.1) { r = st_func(lat, is_north()); @@ -235,7 +235,7 @@ double r, theta; const double H = ( IsNorthHemisphere ? 1.0 : -1.0 ); -if(is_eq(Data.eccentricity, 0.0)) { +if(is_eq(Data.eccentricity, 0.0) || Data.eccentricity >= 0.1) { x = (x - Bx)/Alpha; y = (y - By)/Alpha; @@ -251,7 +251,9 @@ if(is_eq(Data.eccentricity, 0.0)) { lon = Lon_orient - theta; } else { - st_xy_to_latlon_func(x, y, lat, lon, Data.scale_factor, (Data.r_km*m_per_km), + double x1 = Data.x_pin + x * (Data.d_km*m_per_km); // index to meters + double y1 = Data.y_pin + y * (Data.d_km*m_per_km); // index to meters + st_xy_to_latlon_func(x1, y1, lat, lon, Data.scale_factor, (Data.r_km*m_per_km), (-1.0*Data.lon_orient), Data.false_east, Data.false_north, Data.eccentricity, IsNorthHemisphere); } @@ -677,7 +679,7 @@ return ( a ); //////////////////////////////////////////////////////////////////////// - +// lat/lon to meters bool st_latlon_to_xy_func(double lat, double lon, double &x_m, double &y_m, double scale_factor, double scale_lat, double semi_major_axis, @@ -685,6 +687,12 @@ bool st_latlon_to_xy_func(double lat, double lon, double &x_m, double &y_m, double e, bool is_north_hemisphere) { +if (e >= 0.1) { + mlog << Error << "\nst_latlon_to_xy_func() -> " + << "eccentricity (" << e << ") should be less than 0.1 for earth\n\n"; + exit( 1 ); +} +else { const double lat_rad = lat * rad_per_deg; const double lon_rad = lon * rad_per_deg; const double lat_sin = sin(lat_rad); @@ -696,12 +704,12 @@ bool st_latlon_to_xy_func(double lat, double lon, double &x_m, double &y_m, // meters in polar stereographics, not index x_m = false_east + rho*sin(lon_rad - lonO_rad); y_m = false_north - H*rho*cos(lon_rad - lonO_rad); - +} return true; } //////////////////////////////////////////////////////////////////////// - +// meters to lat/lon bool st_xy_to_latlon_func(double x_m, double y_m, double &lat, double &lon, double scale_factor, double semi_major_axis, @@ -709,28 +717,35 @@ bool st_xy_to_latlon_func(double x_m, double y_m, double &lat, double &lon, double eccentricity, bool is_north_hemisphere) { bool result = true; -double chi; -double x_diff = x_m - false_east; -double y_diff = y_m - false_north; -double lonO_rad = proj_vertical_lon * rad_per_deg; -double r_rho = sqrt(x_diff*x_diff + y_diff*y_diff); -double t = r_rho * sqrt(pow((1+eccentricity),(1+eccentricity)) * pow((1-eccentricity),(1-eccentricity))) - / (2*semi_major_axis*scale_factor); -if (is_north_hemisphere) chi = M_PI/2 - 2 * atan(t); -else chi = 2 * atan(t) - M_PI/2; - -lat = chi + (eccentricity*eccentricity/2 + 5*pow(eccentricity,4)/24 - + pow(eccentricity,6)/12 + 13*pow(eccentricity,8)/360*sin(2 * chi) - + (7*pow(eccentricity,4)/48 + 29*pow(eccentricity,6)/240 + 811*pow(eccentricity,8)/11520)*sin(4*chi) - + (7*pow(eccentricity,6)/120 + 81*pow(eccentricity,8)/1120)*sin(6*chi) - + (4279*pow(eccentricity,8)/161280)*sin(8*chi)); - -if (x_m == false_east) lon = lonO_rad; -else if (is_north_hemisphere) lon = lonO_rad + atan2(x_diff,-y_diff); -else lon = lonO_rad + atan2(x_diff,y_diff); - -lat /= rad_per_deg; -lon /= rad_per_deg; +if (eccentricity >= 0.1) { + mlog << Error << "\nst_xy_to_latlon_func() -> " + << "eccentricity (" << eccentricity << ") should be less than 0.1 for earth\n\n"; + exit( 1 ); +} +else { + double chi; + double x_diff = x_m - false_east; + double y_diff = y_m - false_north; + double lonO_rad = proj_vertical_lon * rad_per_deg; + double r_rho = sqrt(x_diff*x_diff + y_diff*y_diff); + double t = r_rho * sqrt(pow((1+eccentricity),(1+eccentricity)) * pow((1-eccentricity),(1-eccentricity))) + / (2*semi_major_axis*scale_factor); + if (is_north_hemisphere) chi = M_PI/2 - 2 * atan(t); + else chi = 2 * atan(t) - M_PI/2; + + lat = chi + (eccentricity*eccentricity/2 + 5*pow(eccentricity,4)/24 + + pow(eccentricity,6)/12 + 13*pow(eccentricity,8)/360*sin(2 * chi) + + (7*pow(eccentricity,4)/48 + 29*pow(eccentricity,6)/240 + 811*pow(eccentricity,8)/11520)*sin(4*chi) + + (7*pow(eccentricity,6)/120 + 81*pow(eccentricity,8)/1120)*sin(6*chi) + + (4279*pow(eccentricity,8)/161280)*sin(8*chi)); + + if (x_m == false_east) lon = lonO_rad; + else if (is_north_hemisphere) lon = lonO_rad + atan2(x_diff,-y_diff); + else lon = lonO_rad + atan2(x_diff,y_diff); + + lat /= rad_per_deg; + lon /= rad_per_deg; +} return result; } From 2dcc27409c568cff02f22934c2b1ce88af43b9b6 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Mon, 23 Jan 2023 15:53:47 -0700 Subject: [PATCH 06/11] #2218 Check eccentricity --- src/libcode/vx_grid/st_grid.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/libcode/vx_grid/st_grid.cc b/src/libcode/vx_grid/st_grid.cc index 6b7db8ec7f..8e0262597f 100644 --- a/src/libcode/vx_grid/st_grid.cc +++ b/src/libcode/vx_grid/st_grid.cc @@ -201,7 +201,7 @@ const double H = ( IsNorthHemisphere ? 1.0 : -1.0 ); reduce(lon); -if(is_eq(Data.eccentricity, 0.0) || Data.eccentricity >= 0.1) { +if(is_eq(Data.eccentricity, 0.0)) { r = st_func(lat, is_north()); @@ -235,7 +235,7 @@ double r, theta; const double H = ( IsNorthHemisphere ? 1.0 : -1.0 ); -if(is_eq(Data.eccentricity, 0.0) || Data.eccentricity >= 0.1) { +if(is_eq(Data.eccentricity, 0.0)) { x = (x - Bx)/Alpha; y = (y - By)/Alpha; From a04888ef9f92d89e09b8b39a20ba221b7ddaa174 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Mon, 23 Jan 2023 15:53:47 -0700 Subject: [PATCH 07/11] #2218 Changed the order of members --- src/libcode/vx_grid/st_grid_defs.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/libcode/vx_grid/st_grid_defs.h b/src/libcode/vx_grid/st_grid_defs.h index ab44bdcb10..47c9b172e7 100644 --- a/src/libcode/vx_grid/st_grid_defs.h +++ b/src/libcode/vx_grid/st_grid_defs.h @@ -41,15 +41,15 @@ struct StereographicData { double r_km; // semi_major_axis + int nx; + int ny; + // ellipsoidal earth double eccentricity; // 0 for shperical earth double false_east; double false_north; double scale_factor; - int nx; - int ny; - void dump(); }; From e452aef42816eb498beecd75c1a323c4e72ca272 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Mon, 23 Jan 2023 15:53:47 -0700 Subject: [PATCH 08/11] #2218 Initialize newly added members --- src/libcode/vx_grid/find_grid_by_name.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/libcode/vx_grid/find_grid_by_name.cc b/src/libcode/vx_grid/find_grid_by_name.cc index 57e5a3de2e..0f018cfeeb 100644 --- a/src/libcode/vx_grid/find_grid_by_name.cc +++ b/src/libcode/vx_grid/find_grid_by_name.cc @@ -450,6 +450,11 @@ if ( !west_longitude_positive ) { } +sdata.eccentricity = 0.; +sdata.false_east = 0.; +sdata.false_north = 0.; +sdata.scale_factor = 1.0; + ToGrid = new Grid ( sdata ); g = *ToGrid; From fe1495fdaa7fa966a07159a7a8d3823cc97a909b Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Mon, 30 Jan 2023 23:12:40 -0700 Subject: [PATCH 09/11] #2218 Added dy_km --- src/libcode/vx_grid/find_grid_by_name.cc | 1 + src/libcode/vx_grid/grid_base.cc | 1 + src/libcode/vx_grid/st_grid_defs.h | 1 + 3 files changed, 3 insertions(+) diff --git a/src/libcode/vx_grid/find_grid_by_name.cc b/src/libcode/vx_grid/find_grid_by_name.cc index 0f018cfeeb..9d8fe9acf1 100644 --- a/src/libcode/vx_grid/find_grid_by_name.cc +++ b/src/libcode/vx_grid/find_grid_by_name.cc @@ -454,6 +454,7 @@ sdata.eccentricity = 0.; sdata.false_east = 0.; sdata.false_north = 0.; sdata.scale_factor = 1.0; +sdata.dy_km = sdata.d_km; ToGrid = new Grid ( sdata ); diff --git a/src/libcode/vx_grid/grid_base.cc b/src/libcode/vx_grid/grid_base.cc index 4dd2604dca..b909547805 100644 --- a/src/libcode/vx_grid/grid_base.cc +++ b/src/libcode/vx_grid/grid_base.cc @@ -1313,6 +1313,7 @@ if ( st1->nx == st2->nx && is_eq (rescale_lon(st1->lon_orient), rescale_lon(st2->lon_orient), loose_tol) && is_eq (st1->d_km, st2->d_km, loose_tol) && + is_eq (st1->dy_km, st2->dy_km, loose_tol) && is_eq (st1->r_km, st2->r_km, loose_tol) ) status = true; return ( status ); diff --git a/src/libcode/vx_grid/st_grid_defs.h b/src/libcode/vx_grid/st_grid_defs.h index 47c9b172e7..dcd5111eeb 100644 --- a/src/libcode/vx_grid/st_grid_defs.h +++ b/src/libcode/vx_grid/st_grid_defs.h @@ -49,6 +49,7 @@ struct StereographicData { double false_east; double false_north; double scale_factor; + double dy_km; void dump(); }; From 6ee856b3864a03cd539563f1354aa327f698a9b2 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Mon, 30 Jan 2023 23:15:25 -0700 Subject: [PATCH 10/11] #2218 Adjusted meters to index at latlon_to_xy & xy_to_latlon --- src/libcode/vx_grid/st_grid.cc | 59 +++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 26 deletions(-) diff --git a/src/libcode/vx_grid/st_grid.cc b/src/libcode/vx_grid/st_grid.cc index 8e0262597f..952568e6b4 100644 --- a/src/libcode/vx_grid/st_grid.cc +++ b/src/libcode/vx_grid/st_grid.cc @@ -212,14 +212,17 @@ if(is_eq(Data.eccentricity, 0.0)) { y = By - Alpha*r*H*cosd(theta); } else { - st_latlon_to_xy_func(lat, lon, x, y, Data.scale_factor, (Data.lon_orient*-1.0), + double delta_sign; + st_latlon_to_xy_func(lat, -lon, x, y, Data.scale_factor, (Data.lon_orient*-1.0), (Data.r_km*m_per_km), Data.false_east, Data.false_north, Data.eccentricity, IsNorthHemisphere); - x = x / (Data.d_km*m_per_km) - Data.x_pin; - y = y / (Data.d_km*m_per_km) - Data.y_pin; -} -return; + delta_sign = ((Data.d_km > 0) ? 1.0 : -1.0 ); + x = delta_sign * ((x / (fabs(Data.d_km)*m_per_km)) - Data.x_pin); // meters to index + delta_sign = ((Data.dy_km > 0) ? 1.0 : -1.0 ); + y = delta_sign * ((y / (fabs(Data.dy_km)*m_per_km)) - Data.y_pin); // meters to index + +} } @@ -251,14 +254,13 @@ if(is_eq(Data.eccentricity, 0.0)) { lon = Lon_orient - theta; } else { - double x1 = Data.x_pin + x * (Data.d_km*m_per_km); // index to meters - double y1 = Data.y_pin + y * (Data.d_km*m_per_km); // index to meters + double x1 = (Data.x_pin + x) * (Data.d_km*m_per_km); // index to meters + double y1 = (Data.y_pin + y) * (Data.dy_km*m_per_km); // index to meters st_xy_to_latlon_func(x1, y1, lat, lon, Data.scale_factor, (Data.r_km*m_per_km), (-1.0*Data.lon_orient), Data.false_east, Data.false_north, Data.eccentricity, IsNorthHemisphere); } - reduce(lon); // @@ -681,36 +683,39 @@ return ( a ); //////////////////////////////////////////////////////////////////////// // lat/lon to meters + bool st_latlon_to_xy_func(double lat, double lon, double &x_m, double &y_m, double scale_factor, double scale_lat, double semi_major_axis, double false_east, double false_north, double e, bool is_north_hemisphere) { -if (e >= 0.1) { - mlog << Error << "\nst_latlon_to_xy_func() -> " - << "eccentricity (" << e << ") should be less than 0.1 for earth\n\n"; - exit( 1 ); -} -else { - const double lat_rad = lat * rad_per_deg; - const double lon_rad = lon * rad_per_deg; - const double lat_sin = sin(lat_rad); - const double lonO_rad = scale_lat * rad_per_deg; - const double H = (is_north_hemisphere? 1.0 : -1.0 ); - double t = tan(M_PI/4 - H*lat_rad/2) * pow((1 + e*lat_sin)/(1 - e*lat_sin), e/2); - double rho = (2 * semi_major_axis * scale_factor * t) - / sqrt(pow(1+e,1+e) * pow(1-e,1-e)); - // meters in polar stereographics, not index - x_m = false_east + rho*sin(lon_rad - lonO_rad); - y_m = false_north - H*rho*cos(lon_rad - lonO_rad); -} + if (e >= 0.1) { + mlog << Error << "\nst_latlon_to_xy_func() -> " + << "eccentricity (" << e << ") should be less than 0.1 for earth\n\n"; + exit( 1 ); + } + else { + const double lat_rad = lat * rad_per_deg; + const double lon_rad = lon * rad_per_deg; + const double lat_sin = sin(lat_rad); + const double lonO_rad = scale_lat * rad_per_deg; + const double H = (is_north_hemisphere? 1.0 : -1.0 ); + double t = tan(M_PI/4 - H*lat_rad/2) * pow((1 + e*lat_sin)/(1 - e*lat_sin), e/2); + double rho = (2 * semi_major_axis * scale_factor * t) + / sqrt(pow(1+e,1+e) * pow(1-e,1-e)); + // meters in polar stereographics, not index + x_m = false_east + rho*sin(lon_rad - lonO_rad); + y_m = false_north - H*rho*cos(lon_rad - lonO_rad); + } + return true; } //////////////////////////////////////////////////////////////////////// // meters to lat/lon + bool st_xy_to_latlon_func(double x_m, double y_m, double &lat, double &lon, double scale_factor, double semi_major_axis, double proj_vertical_lon, double false_east, double false_north, @@ -745,6 +750,8 @@ else { lat /= rad_per_deg; lon /= rad_per_deg; + reduce(lon); + } return result; From e37c837190b3fdc0b1b3a14a28b809194ddd46a1 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Tue, 31 Jan 2023 10:02:19 -0700 Subject: [PATCH 11/11] #2218 Removed a warning for ellipsoidal earth --- src/libcode/vx_data2d_nccf/nccf_file.cc | 80 ++++++++++++++++++------- 1 file changed, 58 insertions(+), 22 deletions(-) diff --git a/src/libcode/vx_data2d_nccf/nccf_file.cc b/src/libcode/vx_data2d_nccf/nccf_file.cc index df3642d469..ad8098d868 100644 --- a/src/libcode/vx_data2d_nccf/nccf_file.cc +++ b/src/libcode/vx_data2d_nccf/nccf_file.cc @@ -2004,8 +2004,8 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va if (!is_eq(inverse_flattening, bad_data_double) || (!is_eq(semi_minor_axis, bad_data_double) && !is_eq(semi_minor_axis,semi_major_axis))) { is_spherical_earch = false; - mlog << Warning << "\n" << method_name - << "This is an ellipsoidal earth which is not fully supported.\n\n"; + mlog << Debug(2) << "\n" << method_name + << "This is an ellipsoidal earth.\n\n"; } else if(!has_scale_factor && !has_standard_parallel) { mlog << Error << "\n" << method_name @@ -2198,7 +2198,7 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va if (is_eq(semi_major_axis, bad_data_double)) semi_major_axis = EARTH_MAJOR_AXIS_km * m_per_km; - else semi_major_axis *= x_coord_to_m_cf; // meters + else if (semi_major_axis < 10000.0) semi_major_axis *= m_per_km; // meters // Calculate the pin indices. The pin will be located at the grid's reference // location since that's the only lat/lon location we know about. @@ -2217,8 +2217,9 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va data.y_pin = y_pin; data.scale_lat = proj_origin_lat; data.lon_orient = -1.0 * proj_vertical_lon; - data.d_km = dx_m / 1000.0; - data.r_km = semi_major_axis / 1000.0; + data.d_km = dx_m / m_per_km; + data.dy_km = dy_m / m_per_km; + data.r_km = semi_major_axis / m_per_km; data.nx = _xDim->getSize(); data.ny = _yDim->getSize(); @@ -2238,28 +2239,33 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va if(!is_spherical_earch) eccentricity = st_eccentricity_func(semi_major_axis, semi_minor_axis, inverse_flattening); - x = (dx_m > 0) ? x_values[0] : x_values[x_counts-1]; - y = (dy_m > 0) ? y_values[0] : y_values[y_counts-1]; + x = x_values[0]; + y = y_values[0]; scale_factor = st_sf_func(proj_standard_parallel, eccentricity, is_north_hemisphere); st_xy_to_latlon_func(x, y, lat, lon, scale_factor, semi_major_axis, proj_vertical_lon, false_east, false_north, eccentricity, is_north_hemisphere); - if(is_eq(eccentricity, 0.0)) { - data.x_pin = 0.; - data.y_pin = 0.; - } - else { - data.x_pin = ((dx_m>0) ? x_values[0] : x_values[x_counts-1]) / dx_m_a; - data.y_pin = ((dy_m>0) ? y_values[0] : y_values[x_counts-1]) / dy_m_a; - } data.lat_pin = lat; data.lon_pin = -lon; + if(!has_scale_factor && has_standard_parallel) { + if(is_eq(eccentricity, 0.0)) { + data.x_pin = data.y_pin = 0.; + } + else { + data.x_pin = x / dx_m_a; + data.y_pin = y / dy_m_a; + } + } + st_latlon_to_xy_func(lat, lon, x2, y2, scale_factor, proj_vertical_lon, semi_major_axis, false_east, false_north, eccentricity, is_north_hemisphere); mlog << Debug(6) << method_name - << "pin: (x,y)m -> (lat,lon) -> (x2,y2)m: (" << x << "," << y << ") -> (" + << " x: between " << x_values[0] <<" and " << x_values[x_counts-1] + << ", y: between " << y_values[0] <<" and " << y_values[y_counts-1] << "\n"; + mlog << Debug(6) << method_name + << "pin: (x,y) -> (lat,lon) -> (x2,y2): (" << x << "," << y << ") -> (" << lat << "," << lon << ") -> (" << x2 << "," << y2 << ") Diff (m): x=" << (x-x2) <<", y=" << (y-y2) << "\n"; @@ -2273,7 +2279,7 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va grid.set(data); - if (dy_m < 0) grid.set_swap_to_north(true); + //Note: do not set grid.set_swap_to_north() if(mlog.verbosity_level() >= 10) { double lat1, lon1; @@ -2291,7 +2297,22 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va semi_major_axis, false_east, false_north, eccentricity, is_north_hemisphere); mlog << Debug(10) << method_name - << "left bottom: (x,y) -> (lat,lon) -> (x2,y2)m: (" + << " bottom left: (x,y) -> (lat,lon) -> (x2,y2): (" + << x1 << "," << y1 << ") -> (" + << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 + << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) + << "\n"; + + x1 = (dx_m > 0) ? x_values[0] : x_values[x_counts-1]; + y1 = (dy_m > 0) ? y_values[y_counts-1] : y_values[0]; + st_xy_to_latlon_func(x1, y1, lat1, lon1, scale_factor, semi_major_axis, + proj_vertical_lon, false_east, false_north, + eccentricity, is_north_hemisphere); + st_latlon_to_xy_func(lat1, lon1, x2, y2, scale_factor, proj_vertical_lon, + semi_major_axis, false_east, false_north, + eccentricity, is_north_hemisphere); + mlog << Debug(10) << method_name + << " top left: (x,y) -> (lat,lon) -> (x2,y2): (" << x1 << "," << y1 << ") -> (" << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) @@ -2306,7 +2327,7 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va semi_major_axis, false_east, false_north, eccentricity, is_north_hemisphere); mlog << Debug(10) << method_name - << " center: (x,y)m -> (lat,lon) -> (x2,y2)m: (" << x1 + << " center: (x,y) -> (lat,lon) -> (x2,y2): (" << x1 << "," << y1 << ") -> (" << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) @@ -2321,7 +2342,22 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va semi_major_axis, false_east, false_north, eccentricity, is_north_hemisphere); mlog << Debug(10) << method_name - << " top-right: (x,y)m -> (lat,lon) -> (x2,y2)m: (" << x1 + << " top right: (x,y) -> (lat,lon) -> (x2,y2): (" << x1 + << "," << y1 << ") -> (" + << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 + << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) + << "\n"; + + x1 = (dx_m > 0) ? x_values[x_counts-1] : x_values[0]; + y1 = (dy_m > 0) ? y_values[0] : y_values[y_counts-1]; + st_xy_to_latlon_func(x1, y1, lat1, lon1, scale_factor, semi_major_axis, + proj_vertical_lon, false_east, false_north, + eccentricity, is_north_hemisphere); + st_latlon_to_xy_func(lat1, lon1, x2, y2, scale_factor, proj_vertical_lon, + semi_major_axis, false_east, false_north, + eccentricity, is_north_hemisphere); + mlog << Debug(10) << method_name + << "bottom right: (x,y) -> (lat,lon) -> (x2,y2): (" << x1 << "," << y1 << ") -> (" << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 << ") Diff (m): x=" << (x1-x2) <<", y=" << (y1-y2) @@ -2348,8 +2384,8 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va data.eccentricity, is_north_hemisphere); x_diff = (x1-x2); y_diff = (y1-y2); - mlog << Debug(15) << method_name - << "index=("<< ix << "," << iy << "): (x,y)m -> (lat,lon) -> (x2,y2)m: (" << x1 + mlog << Debug(15) + << " index=("<< ix << "," << iy << "): (x,y) -> (lat,lon) -> (x2,y2): (" << x1 << "," << y1 << ") -> (" << lat1 << "," << lon1 << ") -> (" << x2 << "," << y2 << ") Diff (m): x=" << x_diff <<", y=" << y_diff