Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Regression in rHealpix projection since v4.9.1 #1206

Closed
pelson opened this Issue Dec 21, 2018 · 7 comments

Comments

Projects
None yet
2 participants
@pelson
Copy link
Contributor

pelson commented Dec 21, 2018

I'm working on providing the Rectangular Healpix projection in cartopy. I started the process some time ago (around Christmas 馃巹last year as it happens) and was using v4.9.1. It produced some pretty pleasant results:

Cartopy healpix projection

(FWIW, cartopy has a general image transformation tool that can produce images like this for any projection that has a forward and inverse defined)

I didn't have the time to investigate at the time, but it turns out that this projection has had a regression in later versions of proj, particularly with regards to the lower rectangle. The following image is indicative of the issue:

figure_1

A quick demonstration of the issue:

$ proj
Rel. 4.9.1, 04 March 2015
usage: proj [ -beEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

echo 40 -50 | proj +proj=rhealpix  | cs2cs +proj=rhealpix  +to +proj=latlon -f %.4f
40.0000	-50.0000 0.0000

And with 5.2.0:

$ proj
Rel. 5.2.0, September 15th, 2018
usage: proj [ -bdeEfiIlmorsStTvVwW [args] ] [ +opts[=arg] ] [ files ]
$ echo 40 -50 | proj +proj=rhealpix  | cs2cs +proj=rhealpix  +to +proj=latlon -f %.4f
-130.0000	50.0000 0.0000

You can see that the round-trip transformation is no longer correct.

The image gives us another hint of what is going on - notice how cartopy is projecting the antarctic geometry to the north pole, which suggests that it is the lon/lat -> healpix meters conversion which is wrong.

@kbevers kbevers added the bug label Dec 21, 2018

@pelson

This comment has been minimized.

Copy link
Contributor Author

pelson commented Dec 21, 2018

Bisect results:

 * 5.2.0 - bad
 * 5.0.1 - bad
 * 4.9.1 - good
 * 95ffe331 - bad
 * e8fd3d30 - good
 * 57af0179 - bad
 * 509d54ee9e2977764ceb641a00747ec0b3ba090a - good
 * 82d0797275be73a50e72016e532efaffdba71bf3 - bad
 * edd1521df96ea5750b59f507fc4e8b9f8dbf3282 - good
 * 307914d451d87afd94b66aed872bba96d1968c85 - good
 * 6a921265db9ff12a263ff1a88118a69a65b7a2df - bad
 * 44e2818ac9dcb9e25d0057d106d3aa07d35c05bd - bad

And the cause: 196a1bb. cc @kbevers

@pelson

This comment has been minimized.

Copy link
Contributor Author

pelson commented Dec 21, 2018

A quick check of the commit, looks like inverting the following change gets some of the way (command line round-trip works again), but I'm still having issues with it within cartopy...:

diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index e1b48daa..6f6680b4 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -479,6 +479,7 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
             tmpRot = rot[get_rotate_index(capmap.cn - pole)];
         } else {
             pole = south_square;
+            a[1] = -M_HALFPI + pole*0;
             tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
         }
     } else {
@@ -492,6 +493,7 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
             tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
         } else {
             pole = south_square;
+            a[1] =  -M_HALFPI  + capmap.cn*0;
             tmpRot = rot[get_rotate_index(capmap.cn - pole)];
         }
     }

196a1bb#r31757442

pelson added a commit to pelson/proj.4 that referenced this issue Dec 21, 2018

@pelson

This comment has been minimized.

Copy link
Contributor Author

pelson commented Dec 21, 2018

Confirming that the full fix is below:

diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index e1b48daa..8aa23222 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -472,26 +472,29 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
     if (inverse == 0) {
         /* Rotate (x, y) about its polar cap tip and then translate it to
            north_square or south_square. */
-        a[0] =  -3.0*M_FORTPI + pole*M_HALFPI;
-        a[1] =  M_HALFPI;
         if (capmap.region == north) {
             pole = north_square;
+            a[0] =  -3.0*M_FORTPI + pole*M_HALFPI;
+            a[1] =  M_HALFPI;
             tmpRot = rot[get_rotate_index(capmap.cn - pole)];
         } else {
             pole = south_square;
+            a[0] =  -3.0*M_FORTPI + pole*M_HALFPI;
+            a[1] =  -M_HALFPI;
             tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
         }
     } else {
         /* Inverse function.
          Unrotate (x, y) and then translate it back. */
         a[0] = -3.0*M_FORTPI + capmap.cn*M_HALFPI;
-        a[1] = M_HALFPI;
         /* disassemble */
         if (capmap.region == north) {
             pole = north_square;
+            a[1] = M_HALFPI;
             tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
         } else {
             pole = south_square;
+            a[1] =  -M_HALFPI;
             tmpRot = rot[get_rotate_index(capmap.cn - pole)];
         }
     }

Unfortunately, applying this to the 5.2 branch is a bit harder, as it has been subsequently factored out.

@kbevers

This comment has been minimized.

Copy link
Member

kbevers commented Dec 21, 2018

What is the fix above applied to? The master branch?

Unfortunately, applying this to the 5.2 branch is a bit harder, as it has been subsequently factored out.

Don't worry, we haven't scheduled any patch releases for the 5.2 branch.

@pelson

This comment has been minimized.

Copy link
Contributor Author

pelson commented Dec 21, 2018

It was applied on top of the bad commit (to try avoiding too many moving parts).
I've put up #1206 against 5.2 which fixes the regression. I haven't added a test yet though.
Let's move over to the PR for implementation discussion.

pelson added a commit to pelson/proj.4 that referenced this issue Dec 21, 2018

@pelson

This comment has been minimized.

Copy link
Contributor Author

pelson commented Dec 24, 2018

Another key change:

diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index e1b48daa..3c6a331a 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -162,8 +162,8 @@ int in_image(double x, double y, int proj, int north_square, int south_square) {
         double healpixVertsJit[][2] = {
             {-1.0*M_PI - EPS, M_FORTPI},
             {-3.0*M_FORTPI, M_HALFPI + EPS},
-            {M_HALFPI, M_FORTPI + EPS},
-            {M_FORTPI, M_HALFPI + EPS},
+            {-M_HALFPI, M_FORTPI + EPS},
+            {-M_FORTPI, M_HALFPI + EPS},
             {0.0, M_FORTPI + EPS},
             {M_FORTPI, M_HALFPI + EPS},
             {M_HALFPI, M_FORTPI + EPS},
@@ -192,7 +192,7 @@ int in_image(double x, double y, int proj, int north_square, int south_square) {
         rhealpixVertsJit[0][0]  = -1.0*M_PI - EPS;
         rhealpixVertsJit[0][1]  = M_FORTPI + EPS;
         rhealpixVertsJit[1][0]  = -1.0*M_PI + north_square*M_HALFPI- EPS;
-        rhealpixVertsJit[1][2]  = M_FORTPI + EPS;
+        rhealpixVertsJit[1][1]  = M_FORTPI + EPS;
         rhealpixVertsJit[2][0]  = -1.0*M_PI + north_square*M_HALFPI- EPS;
         rhealpixVertsJit[2][1]  = 3*M_FORTPI + EPS;
         rhealpixVertsJit[3][0]  = -1.0*M_PI + (north_square + 1.0)*M_HALFPI + EPS;

pelson added a commit to pelson/proj.4 that referenced this issue Dec 24, 2018

Fixed rHealpix projection, closing OSGeo#1206.
Numerical test results based on function output, rather than mathematical derivation. Verified global coverage with graphical eyeballing through cartopy.

kbevers added a commit that referenced this issue Dec 26, 2018

Merge pull request #1207 from pelson/healpix_fix
Fixed rHealpix projection, closing #1206.
@pelson

This comment has been minimized.

Copy link
Contributor Author

pelson commented Dec 26, 2018

Solved by #1207.

@pelson pelson closed this Dec 26, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can鈥檛 perform that action at this time.