Regression in rHealpix projection since v4.9.1

[pelson]

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:

(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:

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.

[pelson]

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]

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]

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]

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]

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]

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]

Solved by #1207.