Merge pull request #1205 from kbevers/master

Better handling of pipelines including step with PJ_IO_UNITS_WHATEVER operations

16 files changed, 86 insertions, 42 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index d4ca51fa..e0d15773 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -68,8 +68,8 @@ int proj_angular_input (PJ *P, enum PJ_DIRECTION dir) {
     dir: {PJ_FWD, PJ_INV}
 ******************************************************************************/
     if (PJ_FWD==dir)
-        return pj_left (P)==PJ_IO_UNITS_ANGULAR;
-    return pj_right (P)==PJ_IO_UNITS_ANGULAR;
+        return pj_left (P)==PJ_IO_UNITS_RADIANS;
+    return pj_right (P)==PJ_IO_UNITS_RADIANS;
 }
 
 /*****************************************************************************/
diff --git a/src/conversions/axisswap.cpp b/src/conversions/axisswap.cpp
index 18d13f8c..97c8899a 100644
--- a/src/conversions/axisswap.cpp
+++ b/src/conversions/axisswap.cpp
@@ -282,8 +282,8 @@ PJ *CONVERSION(axisswap,0) {
     }
 
     if (pj_param(P->ctx, P->params, "tangularunits").i) {
-        P->left  = PJ_IO_UNITS_ANGULAR;
-        P->right = PJ_IO_UNITS_ANGULAR;
+        P->left  = PJ_IO_UNITS_RADIANS;
+        P->right = PJ_IO_UNITS_RADIANS;
     } else {
         P->left  = PJ_IO_UNITS_WHATEVER;
         P->right = PJ_IO_UNITS_WHATEVER;
diff --git a/src/conversions/cart.cpp b/src/conversions/cart.cpp
index f8ec256b..d9aea9b8 100644
--- a/src/conversions/cart.cpp
+++ b/src/conversions/cart.cpp
@@ -213,7 +213,7 @@ PJ *CONVERSION(cart,1) {
     P->inv3d  =  geodetic;
     P->fwd    =  cart_forward;
     P->inv    =  cart_reverse;
-    P->left   =  PJ_IO_UNITS_ANGULAR;
+    P->left   =  PJ_IO_UNITS_RADIANS;
     P->right  =  PJ_IO_UNITS_CARTESIAN;
     return P;
 }
diff --git a/src/conversions/geoc.cpp b/src/conversions/geoc.cpp
index 66ca6199..e0ca3df3 100644
--- a/src/conversions/geoc.cpp
+++ b/src/conversions/geoc.cpp
@@ -51,8 +51,8 @@ static PJ *CONVERSION(geoc, 1) {
     P->inv4d = inverse;
     P->fwd4d = forward;
 
-    P->left   =  PJ_IO_UNITS_ANGULAR;
-    P->right  =  PJ_IO_UNITS_ANGULAR;
+    P->left   =  PJ_IO_UNITS_RADIANS;
+    P->right  =  PJ_IO_UNITS_RADIANS;
 
     P->is_latlong = 1;
     return P;
diff --git a/src/conversions/geocent.cpp b/src/conversions/geocent.cpp
index c8d2d486..31f1a30b 100644
--- a/src/conversions/geocent.cpp
+++ b/src/conversions/geocent.cpp
@@ -56,7 +56,7 @@ PJ *CONVERSION (geocent, 0) {
     P->y0 = 0.0;
     P->inv = inverse;
     P->fwd = forward;
-    P->left = PJ_IO_UNITS_ANGULAR;
+    P->left = PJ_IO_UNITS_RADIANS;
     P->right = PJ_IO_UNITS_CARTESIAN;
 
     return P;
diff --git a/src/conversions/unitconvert.cpp b/src/conversions/unitconvert.cpp
index cbd81a7d..1e3372d6 100644
--- a/src/conversions/unitconvert.cpp
+++ b/src/conversions/unitconvert.cpp
@@ -457,6 +457,8 @@ PJ *CONVERSION(unitconvert,0) {
 
     P->left  = PJ_IO_UNITS_WHATEVER;
     P->right = PJ_IO_UNITS_WHATEVER;
+    P->skip_fwd_prepare = 1;
+    P->skip_inv_prepare = 1;
 
     /* if no time input/output unit is specified we can skip them */
     Q->t_in_id = -1;
@@ -476,6 +478,8 @@ PJ *CONVERSION(unitconvert,0) {
         }
         if (f != 0.0)
             Q->xy_factor *= f;
+        if (normalized_name != nullptr && strcmp(normalized_name, "Radian") == 0)
+            P->left = PJ_IO_UNITS_RADIANS;
     }
 
     if ((name = pj_param (P->ctx, P->params, "sxy_out").s) != nullptr) {
@@ -489,6 +493,8 @@ PJ *CONVERSION(unitconvert,0) {
         }
         if (f != 0.0)
             Q->xy_factor /= f;
+        if (normalized_name != nullptr && strcmp(normalized_name, "Radian") == 0)
+            P->right= PJ_IO_UNITS_RADIANS;
     }
 
     if( xy_in_is_linear >= 0 && xy_out_is_linear >= 0 &&
diff --git a/src/fwd.cpp b/src/fwd.cpp
index 02f31f87..a8c51934 100644
--- a/src/fwd.cpp
+++ b/src/fwd.cpp
@@ -48,7 +48,7 @@ static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
     if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
 
     /* Check validity of angular input coordinates */
-    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
+    if (INPUT_UNITS==PJ_IO_UNITS_RADIANS) {
         double t;
 
         /* check for latitude or longitude over-range */
@@ -135,7 +135,7 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
     case PJ_IO_UNITS_WHATEVER:
         break;
 
-    case PJ_IO_UNITS_ANGULAR:
+    case PJ_IO_UNITS_RADIANS:
         coo.lpz.z = P->vfr_meter * (coo.lpz.z + P->z0);
 
         if( P->is_long_wrap_set ) {
diff --git a/src/inv.cpp b/src/inv.cpp
index 83c5498c..a05f8376 100644
--- a/src/inv.cpp
+++ b/src/inv.cpp
@@ -82,7 +82,7 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
         coo.xyz.y *= P->ra;
         return coo;
 
-    case PJ_IO_UNITS_ANGULAR:
+    case PJ_IO_UNITS_RADIANS:
         coo.lpz.z = P->vto_meter * coo.lpz.z - P->z0;
         break;
     }
@@ -99,7 +99,7 @@ static PJ_COORD inv_finalize (PJ *P, PJ_COORD coo) {
         return proj_coord_error ();
     }
 
-    if (OUTPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
+    if (OUTPUT_UNITS==PJ_IO_UNITS_RADIANS) {
 
         /* Distance from central meridian, taking system zero meridian into account */
         coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
diff --git a/src/pipeline.cpp b/src/pipeline.cpp
index 5f7efe08..b24f71e7 100644
--- a/src/pipeline.cpp
+++ b/src/pipeline.cpp
@@ -329,6 +329,33 @@ static void set_ellipsoid(PJ *P) {
 }
 
 
+static enum pj_io_units get_next_non_whatever_unit(void *pipeline_data, int step, PJ_DIRECTION dir) {
+    PJ **pipeline = static_cast<struct pj_opaque*>(pipeline_data)->pipeline;
+    int nsteps = static_cast<struct pj_opaque*>(pipeline_data)->steps;
+    int i;
+
+    if (dir == PJ_FWD) {
+        for (i = step+1; i<=nsteps; i++) {
+            if (pj_left(pipeline[i]) != pj_right(pipeline[i]))
+                return pj_left(pipeline[i]);
+            if (pj_left(pipeline[i]) != PJ_IO_UNITS_WHATEVER)
+                return pj_left(pipeline[i]);
+            if (pj_right(pipeline[i]) != PJ_IO_UNITS_WHATEVER)
+                return pj_right(pipeline[i]);
+        }
+    } else {
+        for (i=step; i>1; i--) {
+            if (pj_right(pipeline[i]) != pj_left(pipeline[i]))
+                return pj_right(pipeline[i]);
+            if (pj_right(pipeline[i]) != PJ_IO_UNITS_WHATEVER)
+                return pj_right(pipeline[i]);
+            if (pj_left(pipeline[i]) != PJ_IO_UNITS_WHATEVER)
+                return pj_left(pipeline[i]);
+        }
+    }
+    return PJ_IO_UNITS_WHATEVER;
+}
+
 
 
 PJ *OPERATION(pipeline,0) {
@@ -480,14 +507,36 @@ PJ *OPERATION(pipeline,0) {
         }
     }
 
-    /* Check that output units from step i are compatible with expected units in step i+1 */
+
+    /* Replace PJ_IO_UNITS_WHATEVER with input/output units of neighbouring steps where */
+    /* it make sense. It does in most cases but not always, for instance                */
+    /*      proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step ...           */
+    /* where the left-hand side units of the first step shouldn't be changed to RADIANS */
+    /* as it will result in deg->rad conversions in cs2cs and other applications.       */
+    PJ **pipeline = static_cast<struct pj_opaque*>(P->opaque)->pipeline;
+    for (i=1; i<=nsteps; i++) {
+        if (pj_left(pipeline[i]) == PJ_IO_UNITS_WHATEVER && pj_right(pipeline[i]) == PJ_IO_UNITS_WHATEVER) {
+            pipeline[i]->left = get_next_non_whatever_unit(P->opaque, i, PJ_FWD);
+            pipeline[i]->right = get_next_non_whatever_unit(P->opaque, i, PJ_FWD);
+        }
+    }
+
+    for (i=nsteps; i>0; i--) {
+        if (pj_left(pipeline[i]) == PJ_IO_UNITS_WHATEVER && pj_right(pipeline[i]) == PJ_IO_UNITS_WHATEVER) {
+            pipeline[i]->right = get_next_non_whatever_unit(P->opaque, i, PJ_INV);
+            pipeline[i]->left = get_next_non_whatever_unit(P->opaque, i, PJ_INV);
+        }
+    }
+
+    /* Check that units between each steps match each other, fail if they don't */
     for (i = 1; i < nsteps; i++) {
-        enum pj_io_units unit_returned = pj_right (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i]);
-        enum pj_io_units unit_expected = pj_left  (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i+1]);
+        enum pj_io_units curr_step_output = pj_right (pipeline[i]);
+        enum pj_io_units next_step_input  = pj_left  (pipeline[i+1]);
 
-        if ( unit_returned == PJ_IO_UNITS_WHATEVER || unit_expected == PJ_IO_UNITS_WHATEVER )
+        if ( curr_step_output == PJ_IO_UNITS_WHATEVER || next_step_input == PJ_IO_UNITS_WHATEVER )
             continue;
-        if ( unit_returned != unit_expected ) {
+
+        if ( curr_step_output != next_step_input ) {
             proj_log_error (P, "Pipeline: Mismatched units between step %d and %d", i, i+1);
             return destructor (P, PJD_ERR_MALFORMED_PIPELINE);
         }
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 76b95130..63e10407 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -172,15 +172,12 @@ typedef long pj_int32;
 #define DIR_CHAR '/'
 #endif
 
-
-/* This enum is also conditionally defined in projects.h - but enums cannot */
-/* be forward declared and we need it here for the pj_left/right prototypes */
 enum pj_io_units {
     PJ_IO_UNITS_WHATEVER  = 0,  /* Doesn't matter (or depends on pipeline neighbours) */
     PJ_IO_UNITS_CLASSIC   = 1,  /* Scaled meters (right), projected system */
     PJ_IO_UNITS_PROJECTED = 2,  /* Meters, projected system */
     PJ_IO_UNITS_CARTESIAN = 3,  /* Meters, 3D cartesian system */
-    PJ_IO_UNITS_ANGULAR   = 4   /* Radians */
+    PJ_IO_UNITS_RADIANS   = 4   /* Radians */
 };
 enum pj_io_units pj_left (PJ *P);
 enum pj_io_units pj_right (PJ *P);
@@ -245,15 +242,7 @@ struct ARG_list;
 struct PJ_REGION_S;
 typedef struct PJ_REGION_S  PJ_Region;
 typedef struct ARG_list paralist;   /* parameter list */
-#ifndef PROJ_INTERNAL_H
-enum pj_io_units {
-    PJ_IO_UNITS_WHATEVER  = 0,  /* Doesn't matter (or depends on pipeline neighbours) */
-    PJ_IO_UNITS_CLASSIC   = 1,  /* Scaled meters (right), projected system */
-    PJ_IO_UNITS_PROJECTED = 2,  /* Meters, projected system */
-    PJ_IO_UNITS_CARTESIAN = 3,  /* Meters, 3D cartesian system */
-    PJ_IO_UNITS_ANGULAR   = 4   /* Radians */
-};
-#endif
+
 #ifndef PROJ_H
 typedef struct PJconsts PJ;         /* the PJ object herself */
 typedef union  PJ_COORD PJ_COORD;
@@ -692,7 +681,7 @@ C_NAMESPACE PJ *pj_##name (PJ *P) {                          \
         return nullptr;                                      \
     P->descr = des_##name;                                   \
     P->need_ellps = NEED_ELPJ_LPS;                              \
-    P->left  = PJ_IO_UNITS_ANGULAR;                          \
+    P->left  = PJ_IO_UNITS_RADIANS;                          \
     P->right = PJ_IO_UNITS_CLASSIC;                          \
     return P;                                                \
 }                                                            \
diff --git a/src/projections/latlong.cpp b/src/projections/latlong.cpp
index f114504b..970c4893 100644
--- a/src/projections/latlong.cpp
+++ b/src/projections/latlong.cpp
@@ -98,8 +98,8 @@ static PJ *latlong_setup (PJ *P) {
     P->fwd3d = latlong_forward_3d;
     P->inv4d = latlong_inverse_4d;
     P->fwd4d = latlong_forward_4d;
-    P->left = PJ_IO_UNITS_ANGULAR;
-    P->right = PJ_IO_UNITS_ANGULAR;
+    P->left = PJ_IO_UNITS_RADIANS;
+    P->right = PJ_IO_UNITS_RADIANS;
     return P;
 }
 
diff --git a/src/projections/ob_tran.cpp b/src/projections/ob_tran.cpp
index 1726b622..6daae394 100644
--- a/src/projections/ob_tran.cpp
+++ b/src/projections/ob_tran.cpp
@@ -237,7 +237,7 @@ PJ *PROJECTION(ob_tran) {
 
     /* Support some rather speculative test cases, where the rotated projection */
     /* is actually latlong. We do not want scaling in that case... */
-    if (Q->link->right==PJ_IO_UNITS_ANGULAR)
+    if (Q->link->right==PJ_IO_UNITS_RADIANS)
         P->right = PJ_IO_UNITS_PROJECTED;
 
 
diff --git a/src/transformations/affine.cpp b/src/transformations/affine.cpp
index 21529a20..bda54f1e 100644
--- a/src/transformations/affine.cpp
+++ b/src/transformations/affine.cpp
@@ -238,8 +238,8 @@ PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) {
     P->fwd    = forward_2d;
     P->inv    = reverse_2d;
 
-    P->left   = PJ_IO_UNITS_ANGULAR;
-    P->right  = PJ_IO_UNITS_ANGULAR;
+    P->left   = PJ_IO_UNITS_RADIANS;
+    P->right  = PJ_IO_UNITS_RADIANS;
 
     /* read args */
     Q->xoff = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD;
diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp
index 788c2ebb..2e2294cb 100644
--- a/src/transformations/hgridshift.cpp
+++ b/src/transformations/hgridshift.cpp
@@ -93,8 +93,8 @@ PJ *TRANSFORMATION(hgridshift,0) {
     P->fwd    = nullptr;
     P->inv    = nullptr;
 
-    P->left  = PJ_IO_UNITS_ANGULAR;
-    P->right = PJ_IO_UNITS_ANGULAR;
+    P->left  = PJ_IO_UNITS_RADIANS;
+    P->right = PJ_IO_UNITS_RADIANS;
 
     if (0==pj_param(P->ctx, P->params, "tgrids").i) {
         proj_log_error(P, "hgridshift: +grids parameter missing.");
diff --git a/src/transformations/molodensky.cpp b/src/transformations/molodensky.cpp
index 002c5286..c389fd32 100644
--- a/src/transformations/molodensky.cpp
+++ b/src/transformations/molodensky.cpp
@@ -287,8 +287,8 @@ PJ *TRANSFORMATION(molodensky,1) {
     P->fwd    = forward_2d;
     P->inv    = reverse_2d;
 
-    P->left   = PJ_IO_UNITS_ANGULAR;
-    P->right  = PJ_IO_UNITS_ANGULAR;
+    P->left   = PJ_IO_UNITS_RADIANS;
+    P->right  = PJ_IO_UNITS_RADIANS;
 
     /* read args */
     if (pj_param(P->ctx, P->params, "tdx").i) {
diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp
index 95537cda..fda38ec3 100644
--- a/src/transformations/vgridshift.cpp
+++ b/src/transformations/vgridshift.cpp
@@ -137,8 +137,8 @@ PJ *TRANSFORMATION(vgridshift,0) {
     P->fwd    = nullptr;
     P->inv    = nullptr;
 
-    P->left  = PJ_IO_UNITS_ANGULAR;
-    P->right = PJ_IO_UNITS_ANGULAR;
+    P->left  = PJ_IO_UNITS_RADIANS;
+    P->right = PJ_IO_UNITS_RADIANS;
 
     return P;
 }