Skip to content

Commit

Permalink
#2218 Removed a warning for ellipsoidal earth
Browse files Browse the repository at this point in the history
  • Loading branch information
Howard Soh committed Jan 31, 2023
1 parent 6ee856b commit e37c837
Showing 1 changed file with 58 additions and 22 deletions.
80 changes: 58 additions & 22 deletions src/libcode/vx_data2d_nccf/nccf_file.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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();

Expand All @@ -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";

Expand All @@ -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;
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down

0 comments on commit e37c837

Please sign in to comment.