Add rotation support to the HEALPix projection (#1638)

3 files changed, 54 insertions, 4 deletions

diff --git a/docs/source/operations/projections/healpix.rst b/docs/source/operations/projections/healpix.rst
index e13c5456..d28afa8d 100644
--- a/docs/source/operations/projections/healpix.rst
+++ b/docs/source/operations/projections/healpix.rst
@@ -51,6 +51,15 @@ Parameters
 
 .. note:: All parameters for the projection are optional.
 
+.. option:: +rot_xy
+
+    .. versionadded:: 7.0.0
+
+    Rotation of the HEALPix map in degrees. A positive value results in a clockwise
+    rotation around (x_0, y_0) in the cartesian / projected coordinate space.
+
+    *Defaults to 0.0.*
+
 .. include:: ../options/lon_0.rst
 
 .. include:: ../options/x_0.rst
diff --git a/src/projections/healpix.cpp b/src/projections/healpix.cpp
index 515f4f6f..4acb4f69 100644
--- a/src/projections/healpix.cpp
+++ b/src/projections/healpix.cpp
@@ -36,7 +36,7 @@
 #include "proj.h"
 #include "proj_internal.h"
 
-PROJ_HEAD(healpix, "HEAPJ_LPix") "\n\tSph&Ell";
+PROJ_HEAD(healpix, "HEAPJ_LPix") "\n\tSph&Ell\n\trot_xy=";
 PROJ_HEAD(rhealpix, "rHEAPJ_LPix") "\n\tSph&Ell\n\tnorth_square= south_square=";
 
 /* Matrix for counterclockwise rotation by pi/2: */
@@ -56,6 +56,7 @@ namespace { // anonymous namespace
 struct pj_opaque {
     int north_square;
     int south_square;
+    double rot_xy;
     double qp;
     double *apa;
 };
@@ -78,6 +79,12 @@ static double sign (double v) {
     return v > 0 ? 1 : (v < 0 ? -1 : 0);
 }
 
+static PJ_XY rotate(PJ_XY p, double angle) {
+    PJ_XY result;
+    result.x = p.x * cos(angle) - p.y * sin(angle);
+    result.y = p.y * cos(angle) + p.x * sin(angle);
+    return result;
+}
 
 /**
  * Return the index of the matrix in ROT.
@@ -178,7 +185,8 @@ static int in_image(double x, double y, int proj, int north_square,
             {-M_FORTPI,   -M_HALFPI - EPS},
             {-M_HALFPI,   -M_FORTPI - EPS},
             {-3*M_FORTPI, -M_HALFPI - EPS},
-            {-M_PI - EPS, -M_FORTPI}
+            {-M_PI - EPS, -M_FORTPI},
+            {-M_PI - EPS,  M_FORTPI}
         };
         return pnpoly((int)sizeof(healpixVertsJit)/
                       sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
@@ -513,17 +521,22 @@ static PJ_XY combine_caps(double x, double y, int north_square, int south_square
 
 static PJ_XY s_healpix_forward(PJ_LP lp, PJ *P) { /* sphere  */
     (void) P;
-    return healpix_sphere(lp);
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    return rotate(healpix_sphere(lp), -Q->rot_xy);
 }
 
 
 static PJ_XY e_healpix_forward(PJ_LP lp, PJ *P) { /* ellipsoid  */
     lp.phi = auth_lat(P, lp.phi, 0);
-    return healpix_sphere(lp);
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    return rotate(healpix_sphere(lp), -Q->rot_xy);
 }
 
 
 static PJ_LP s_healpix_inverse(PJ_XY xy, PJ *P) { /* sphere */
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    xy = rotate(xy, Q->rot_xy);
+
     /* Check whether (x, y) lies in the HEAPJ_LPix image */
     if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
         PJ_LP lp;
@@ -538,6 +551,8 @@ static PJ_LP s_healpix_inverse(PJ_XY xy, PJ *P) { /* sphere */
 
 static PJ_LP e_healpix_inverse(PJ_XY xy, PJ *P) { /* ellipsoid */
     PJ_LP lp = {0.0,0.0};
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    xy = rotate(xy, Q->rot_xy);
 
     /* Check whether (x, y) lies in the HEAPJ_LPix image. */
     if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
@@ -622,6 +637,9 @@ PJ *PROJECTION(healpix) {
     P->opaque = Q;
     P->destructor = destructor;
 
+    double angle = pj_param(P->ctx, P->params,"drot_xy").f;
+    Q->rot_xy = PJ_TORAD(angle);
+
     if (P->es != 0.0) {
         Q->apa = pj_authset(P->es);             /* For auth_lat(). */
         if (nullptr==Q->apa)
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index a70baba5..f86b7e79 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -2021,6 +2021,29 @@ expect  -0.001790493 0.000759909
 accept  -200 -100
 expect  -0.001790493 -0.000759909
 
+-------------------------------------------------------------------------------
+operation +proj=healpix   +R=6400000   +lat_1=0.5 +lat_2=2  +rot_xy=42
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  2 1
+expect  254069.735470912856 -51696.237925639456
+accept  2 -1
+expect  77970.559536809917 -247274.186569161975
+accept  -2 1
+expect  -77970.559536809917 247274.186569161975
+accept  -2 -1
+expect  -254069.735470912856 51696.237925639456
+
+direction inverse
+accept  254069.735470912856 -51696.237925639456
+expect  2 1
+accept  77970.559536809917 -247274.186569161975
+expect  2 -1
+accept  -77970.559536809917 247274.186569161975
+expect  -2 1
+accept  -254069.735470912856 51696.237925639456
+expect  -2 -1
+
 
 ===============================================================================
 rHEALPix