Skip to content

Commit

Permalink
Fixes opposite pole coordinates (#202)
Browse files Browse the repository at this point in the history 
* Fixes opposite pole coordinates
* Rerun tests with larger tolerance
  • Loading branch information
@benjaminmenetrier
benjaminmenetrier authored Jun 17, 2024
1 parent ed9eb7e commit 3f0055c
Show file tree
Hide file tree
Showing 11 changed files with 3,431 additions and 1,772 deletions.
6 changes: 3 additions & 3 deletions doc/example-grids/classic_gaussian_2.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ projection :

check :
size : 1688
lonlat(first) : [3,44.8796]
lonlat(last) : [-172.453,-54.9736]
uid : 64d609c6ee4b036b209047aef97a10eb
lonlat(first) : [3,49.1204]
lonlat(last) : [179.6485,-38.8936]
uid : 02ef7ed866af557e04882342a90fddfe
bounding_box(n,w,s,e) : [90,0,-90,360]
2 changes: 1 addition & 1 deletion doc/example-grids/custom_structured_5.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,5 @@ check :
size : 144
lonlat(first) : [3,40]
lonlat(last) : [-177,-40]
uid : d55211f27c4a8c6ce94ea67f0c706e22
uid : 23f0a23085bbea10277b457f482927c0
bounding_box(n,w,s,e) : [90,0,-90,360]
6 changes: 3 additions & 3 deletions doc/example-grids/octahedral_gaussian_2.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ projection : { type: "rotated_schmidt", stretching_factor : 4.0, north_pole : [

check :
size : 1600
lonlat(first) : [3,45.9397]
lonlat(last) : [-165.776,-62.6128]
uid : 20308244515b7c78e22709ebaa842246
lonlat(first) : [3,48.0603]
lonlat(last) : [177.0166,-30.8]
uid : 1297f9ab088b831da3139eb5f44e20a9
bounding_box(n,w,s,e) : [90,0,-90,360]
6 changes: 3 additions & 3 deletions doc/example-grids/octahedral_gaussian_4.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ projection : { type: "rotated_schmidt", stretching_factor: 4.0, north_pole: [3.

check :
size : 1600
lonlat(first) : [3,45.9397]
lonlat(last) : [-165.776,-62.6128]
uid : 20308244515b7c78e22709ebaa842246
lonlat(first) : [3,48.0603]
lonlat(last) : [177.0166,-30.8]
uid : 1297f9ab088b831da3139eb5f44e20a9
bounding_box(n,w,s,e) : [90,0,-90,360]
2 changes: 1 addition & 1 deletion doc/example-grids/regional_rotated_mercator_1.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,5 +14,5 @@ check :
size : 2000
lonlat(first) : [-10.319,40.2109]
lonlat(last) : [24.4206,57.2263]
uid : 9d4b8c38db08c6f3fe3221c8770d8d57
uid : 88d856331d46d30703b38e3ed2e8b2de
bounding_box(n,w,s,e) : [58.734,-16.4206,40.2109,24.4206]
6 changes: 3 additions & 3 deletions doc/example-grids/regular_gaussian_2.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ projection :

check :
size : 128
lonlat(first) : [3,38.8589]
lonlat(last) : [-135.011,-72.4668]
uid : d179fd0b56811242a1a0a8e83dade6a1
lonlat(first) : [3,55.1411]
lonlat(last) : [170.8438,-16.8513]
uid : 2efc83bfae617ec2fe1a7cebe523cf41
bounding_box(n,w,s,e) : [90,0,-90,360]
18 changes: 9 additions & 9 deletions src/atlas/util/Rotation.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,9 +126,9 @@ void Rotation::precompute() {

Rotation::Rotation(const PointLonLat& south_pole, double rotation_angle) {
spole_ = south_pole;
npole_ = PointLonLat(spole_.lon() - 180., spole_.lat() + 180.);
if (npole_.lat() > 90) {
npole_.lon() += 180.;
npole_ = PointLonLat(spole_.lon() - 180., -spole_.lat());
if (npole_.lon() < 0.) {
npole_.lon() += 360.;
}
angle_ = wrap_angle(rotation_angle);

Expand All @@ -144,16 +144,16 @@ Rotation::Rotation(const eckit::Parametrisation& p) {
std::vector<double> pole(2);
if (p.get("north_pole", pole)) {
npole_ = PointLonLat(pole.data());
spole_ = PointLonLat(npole_.lon() + 180., npole_.lat() - 180.);
if (spole_.lat() < -90) {
spole_.lon() -= 180.;
spole_ = PointLonLat(npole_.lon() - 180., -npole_.lat());
if (spole_.lon() < 0.) {
spole_.lon() += 360.;
}
}
else if (p.get("south_pole", pole)) {
spole_ = PointLonLat(pole.data());
npole_ = PointLonLat(spole_.lon() - 180., spole_.lat() + 180.);
if (npole_.lat() > 90) {
npole_.lon() += 180.;
npole_ = PointLonLat(spole_.lon() - 180., -spole_.lat());
if (npole_.lon() < 0.) {
npole_.lon() += 360.;
}
}

Expand Down
2,563 changes: 1,689 additions & 874 deletions src/tests/grid/fctest_stretchedrotatedgaussiangrid.F90
View file

Large diffs are not rendered by default.

2,537 changes: 1,690 additions & 847 deletions src/tests/grid/test_stretchedrotatedgaussian.cc
View file

Large diffs are not rendered by default.

8 changes: 4 additions & 4 deletions src/tests/projection/fctest_projection.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ module projection_fixture
call config%set("north_pole", [2.0_dp,46.7_dp] )
projection = atlas_Projection(config)
FCTEST_CHECK_EQUAL( projection%type(), "rotated_schmidt" )
FCTEST_CHECK_EQUAL( projection%hash(), "148d7ceb58250c0f48dc6b590941a341" )
FCTEST_CHECK_EQUAL( projection%hash(), "3a39e0635b7d0a45f684696ca89825e6" )
call config%final()
call projection%final()
END_TEST
Expand All @@ -91,7 +91,7 @@ module projection_fixture
projection = atlas_RotatedSchmidtProjection( stretching_factor=2.0_dp, &
north_pole=[2.0_dp,46.7_dp], rotation_angle=180.0_dp)
FCTEST_CHECK_EQUAL( projection%type(), "rotated_schmidt" )
FCTEST_CHECK_EQUAL( projection%hash(), "148d7ceb58250c0f48dc6b590941a341" )
FCTEST_CHECK_EQUAL( projection%hash(), "3a39e0635b7d0a45f684696ca89825e6" )
call projection%final()
END_TEST

Expand Down Expand Up @@ -130,7 +130,7 @@ module projection_fixture
call config%set("rotation_angle", 180.0_dp)
projection = atlas_Projection(config)
FCTEST_CHECK_EQUAL( projection%type(), "rotated_lonlat" )
FCTEST_CHECK_EQUAL( projection%hash(), "2b6db0e1ccbe7c514dd726f408f92adb" )
FCTEST_CHECK_EQUAL( projection%hash(), "79586cfbc8145cdef1a25d075a9ae07e" )
call config%final()
call projection%final()
END_TEST
Expand All @@ -139,7 +139,7 @@ module projection_fixture
type(atlas_Projection) :: projection
projection = atlas_RotatedLonLatProjection([2.0_dp,46.7_dp],180._dp)
FCTEST_CHECK_EQUAL( projection%type(), "rotated_lonlat" )
FCTEST_CHECK_EQUAL( projection%hash(), "2b6db0e1ccbe7c514dd726f408f92adb" )
FCTEST_CHECK_EQUAL( projection%hash(), "79586cfbc8145cdef1a25d075a9ae07e" )
call projection%final()
END_TEST

Expand Down
49 changes: 25 additions & 24 deletions src/tests/projection/test_rotation.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,8 @@ bool equivalent(const PointLonLat& p1, const PointLonLat& p2) {
auto f = [=](double lon) { return 10. + std::cos(lon * d2r); };

return is_approximately_equal(p1.lat(), p2.lat(), eps) &&
(std::abs(p2.lat()) == 90. || is_approximately_equal(f(p1.lon()), f(p2.lon()), eps));
(is_approximately_equal(std::abs(p2.lat()), 90., eps) ||
is_approximately_equal(f(p1.lon()), f(p2.lon()), eps));
}

#define EXPECT_EQUIVALENT(p1, p2) EXPECT(equivalent(p1, p2))
Expand Down Expand Up @@ -127,27 +128,27 @@ CASE("test_rotation_angle") {
Rotation rot(Config("rotation_angle", 180.) | Config("north_pole", std::vector<double>{2.0, 46.7}));

const PointLonLat ref[] = {
{-178.00000, -46.70000}, {-178.00000, -16.70000}, {-178.00000, 13.30000}, {-178.00000, 43.30000},
{-178.00000, 73.30000}, {2.00000, 76.70000}, {2.00000, 46.70000}, {-178.00000, -46.70000},
{-162.62343, -19.46929}, {-152.02366, 8.65459}, {-139.57464, 36.43683}, {-113.10894, 61.43199},
{-39.88245, 68.00825}, {2.00000, 46.70000}, {-178.00000, -46.70000}, {-148.83443, -27.31067},
{-129.26346, -3.83700}, {-110.79116, 20.05422}, {-85.87917, 41.36507}, {-44.42496, 53.29508},
{2.00000, 46.70000}, {-178.00000, -46.70000}, {-137.90794, -39.07002}, {-109.60146, -21.33906},
{-88.00000, 0.00000}, {-66.39854, 21.33906}, {-38.09206, 39.07002}, {2.00000, 46.70000},
{-178.00000, -46.70000}, {-131.57504, -53.29508}, {-90.12083, -41.36507}, {-65.20884, -20.05422},
{-46.73654, 3.83700}, {-27.16557, 27.31067}, {2.00000, 46.70000}, {-178.00000, -46.70000},
{-136.11755, -68.00825}, {-62.89106, -61.43199}, {-36.42536, -36.43683}, {-23.97634, -8.65459},
{-13.37657, 19.46929}, {2.00000, 46.70000}, {-178.00000, -46.70000}, {-178.00000, -76.70000},
{2.00000, -73.30000}, {2.00000, -43.30000}, {2.00000, -13.30000}, {2.00000, 16.70000},
{2.00000, 46.70000}, {-178.00000, -46.70000}, {140.11755, -68.00825}, {66.89106, -61.43199},
{40.42536, -36.43683}, {27.97634, -8.65459}, {17.37657, 19.46929}, {2.00000, 46.70000},
{-178.00000, -46.70000}, {135.57504, -53.29508}, {94.12083, -41.36507}, {69.20884, -20.05422},
{50.73654, 3.83700}, {31.16557, 27.31067}, {2.00000, 46.70000}, {-178.00000, -46.70000},
{141.90794, -39.07002}, {113.60146, -21.33906}, {92.00000, 0.00000}, {70.39854, 21.33906},
{42.09206, 39.07002}, {2.00000, 46.70000}, {-178.00000, -46.70000}, {152.83443, -27.31067},
{133.26346, -3.83700}, {114.79116, 20.05422}, {89.87917, 41.36507}, {48.42496, 53.29508},
{2.00000, 46.70000}, {-178.00000, -46.70000}, {166.62343, -19.46929}, {156.02366, 8.65459},
{143.57464, 36.43683}, {117.10894, 61.43199}, {43.88245, 68.00825}, {2.00000, 46.70000},
{-178.00000,-46.70000}, {-178.00000,-76.70000}, {2.00000,-73.30000}, {2.00000,-43.30000},
{2.00000,-13.30000}, {2.00000,16.70000}, {2.00000,46.70000}, {-178.00000,-46.70000},
{140.11755,-68.00825}, {66.89106,-61.43199}, {40.42536,-36.43683}, {27.97634,-8.65459},
{17.37657,19.46929}, {2.00000,46.70000}, {-178.00000,-46.70000}, {135.57504,-53.29508},
{94.12083,-41.36507}, {69.20884,-20.05422}, {50.73654,3.83700}, {31.16557,27.31067},
{2.00000,46.70000}, {-178.00000,-46.70000}, {141.90794,-39.07002}, {113.60146,-21.33906},
{92.00000,0.00000}, {70.39854,21.33906}, {42.09206,39.07002}, {2.00000,46.70000},
{-178.00000,-46.70000}, {152.83443,-27.31067}, {133.26346,-3.83700}, {114.79116,20.05422},
{89.87917,41.36507}, {48.42496,53.29508}, {2.00000,46.70000}, {-178.00000,-46.70000},
{166.62343,-19.46929}, {156.02366,8.65459}, {143.57464,36.43683}, {117.10894,61.43199},
{43.88245,68.00825}, {2.00000,46.70000}, {-178.00000,-46.70000}, {-178.00000,-16.70000},
{-178.00000,13.30000}, {-178.00000,43.30000}, {-178.00000,73.30000}, {2.00000,76.70000},
{2.00000,46.70000}, {-178.00000,-46.70000}, {-162.62343,-19.46929}, {-152.02366,8.65459},
{-139.57464,36.43683}, {-113.10894,61.43199}, {-39.88245,68.00825}, {2.00000,46.70000},
{-178.00000,-46.70000}, {-148.83443,-27.31067}, {-129.26346,-3.83700}, {-110.79116,20.05422},
{-85.87917,41.36507}, {-44.42496,53.29508}, {2.00000,46.70000}, {-178.00000,-46.70000},
{-137.90794,-39.07002}, {-109.60146,-21.33906}, {-88.00000,0.00000}, {-66.39854,21.33906},
{-38.09206,39.07002}, {2.00000,46.70000}, {-178.00000,-46.70000}, {-131.57504,-53.29508},
{-90.12083,-41.36507}, {-65.20884,-20.05422}, {-46.73654,3.83700}, {-27.16557,27.31067},
{2.00000,46.70000}, {-178.00000,-46.70000}, {-136.11755,-68.00825}, {-62.89106,-61.43199},
{-36.42536,-36.43683}, {-23.97634,-8.65459}, {-13.37657,19.46929}, {2.00000,46.70000},
};


Expand Down Expand Up @@ -200,14 +201,14 @@ CASE("test_rotation") {
EXPECT_EQUIVALENT(magics.unrotate(r), p);

p = {0., 0.};
r = {-176., -50.};
r = {4., 50.};
EXPECT_EQUIVALENT(rotation.rotate(p), r);
EXPECT_EQUIVALENT(magics.rotate(p), r);
EXPECT_EQUIVALENT(rotation.unrotate(r), p);
EXPECT_EQUIVALENT(magics.unrotate(r), p);

p = {-180., 45.};
r = {-176., 85.};
r = {-176., -5.};
EXPECT_EQUIVALENT(rotation.rotate(p), r);
EXPECT_EQUIVALENT(magics.rotate(p), r);
EXPECT_EQUIVALENT(rotation.unrotate(r), p);
Expand Down

0 comments on commit 3f0055c

Please # to comment.