Merge pull request #1054 from rouault/add_angular_units_to_unitconvert

Add support for deg, rad and grad in unitconvert (fixes #1052), and document that it supports numeric factors (refs #1053)

6 files changed, 164 insertions, 38 deletions

diff --git a/docs/source/operations/conversions/unitconvert.rst b/docs/source/operations/conversions/unitconvert.rst
index 6c8605fe..fab8026f 100644
--- a/docs/source/operations/conversions/unitconvert.rst
+++ b/docs/source/operations/conversions/unitconvert.rst
@@ -6,7 +6,7 @@ Unit conversion
 
 .. versionadded:: 5.0.0
 
-Convert between various distance and time units.
+Convert between various distance, angular and time units.
 
 +---------------------+--------------------------------------------------------+
 | **Alias**           | unitconvert                                            |
@@ -41,24 +41,32 @@ expected to be in units of decimalyears. This can be fixed with `unitconvert`::
 Parameters
 ################################################################################
 
-.. option:: +xy_in=<unit>
+.. option:: +xy_in=<unit> or <conversion_factor>
 
-    Horizontal input units. See :ref:`distance_units` for a list of available
+    Horizontal input units. See :ref:`distance_units` and :ref:`angular_units`
+    for a list of available units. `<conversion_factor>` is the conversion factor
+    from the input unit to metre for linear units, or to radian for angular
     units.
 
-.. option:: +xy_out=<unit>
+.. option:: +xy_out=<unit> or <conversion_factor>
 
-    Horizontal output units. See :ref:`distance_units` for a list of available
+    Horizontal output units. See :ref:`distance_units` and :ref:`angular_units`
+    for a list of available units. `<conversion_factor>` is the conversion factor
+    from the output unit to metre for linear units, or to radian for angular
     units.
 
-.. option:: +z_in=<unit>
+.. option:: +z_in=<unit> or <conversion_factor>
 
-    Vertical output units. See :ref:`distance_units` for a list of available
+    Vertical output units. See :ref:`distance_units` and :ref:`angular_units`
+    for a list of available units. `<conversion_factor>` is the conversion factor
+    from the input unit to metre for linear units, or to radian for angular
     units.
 
-.. option:: +z_out=<unit>
+.. option:: +z_out=<unit> or <conversion_factor>
 
-    Vertical output units. See :ref:`distance_units` for a list of available
+    Vertical output units. See :ref:`distance_units` and :ref:`angular_units`
+    for a list of available units.  `<conversion_factor>` is the conversion factor
+    from the output unit to metre for linear units, or to radian for angular
     units.
 
 .. option:: +t_in=<unit>
@@ -74,7 +82,7 @@ Parameters
 Distance units
 ###############################################################################
 
-In the table below all distance units supported by PROJ is listed.
+In the table below all distance units supported by PROJ are listed.
 The same list can also be produced on the command line with :program:`proj` or
 :program:`cs2cs`, by adding the `-lu` flag when calling the utility.
 
@@ -124,12 +132,31 @@ The same list can also be produced on the command line with :program:`proj` or
 | ind-ch   | Indian Chain                    |
 +----------+---------------------------------+
 
+.. _angular_units:
+
+Angular units
+###############################################################################
+
+.. versionadded:: 5.2.0
+
+In the table below all angular units supported by PROJ `unitconvert` are listed.
+
++----------+---------------------------------+
+| **Label**| **Name**                        |
++----------+---------------------------------+
+|    deg   | Degree                          |
++----------+---------------------------------+
+|   grad   | Grad                            |
++----------+---------------------------------+
+|    rad   | Radian                          |
++----------+---------------------------------+
+
 .. _time_units:
 
 Time units
 ###############################################################################
 
-In the table below all time units supported by PROJ is listed.
+In the table below all time units supported by PROJ are listed.
 
 +--------------+-----------------------------+
 | **label**    | **Name**                    |
diff --git a/src/PJ_unitconvert.c b/src/PJ_unitconvert.c
index 25bdc2e7..6e53f6df 100644
--- a/src/PJ_unitconvert.c
+++ b/src/PJ_unitconvert.c
@@ -382,6 +382,59 @@ static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
     return out;
 }
 
+/* M_PI / 200 */
+#define GRAD_TO_RAD 0.015707963267948967
+
+static const struct PJ_UNITS
+pj_angular_units[] = {
+    {"rad",         "1.0",                   "Radian", 1.0},
+    {"deg",         "0.017453292519943296",  "Degree", DEG_TO_RAD},
+    {"grad",        "0.015707963267948967",  "Grad",   GRAD_TO_RAD},
+    {NULL,          NULL,                     NULL,    0.0}
+};
+
+
+/***********************************************************************/
+static double get_unit_conversion_factor(const char* name,
+                                         int* p_is_linear,
+                                         const char** p_normalized_name) {
+/***********************************************************************/
+    int i;
+    const char* s;
+
+    /* Try first with linear units */
+    for (i = 0; (s = pj_units[i].id) ; ++i) {
+        if ( strcmp(s, name) == 0 ) {
+            if( p_normalized_name ) {
+                *p_normalized_name = pj_units[i].name;
+            }
+            if( p_is_linear ) {
+                *p_is_linear = 1;
+            }
+            return pj_units[i].factor;
+        }
+    }
+
+    /* And then angular units */
+    for (i = 0; (s = pj_angular_units[i].id) ; ++i) {
+        if ( strcmp(s, name) == 0 ) {
+            if( p_normalized_name ) {
+                *p_normalized_name = pj_angular_units[i].name;
+            }
+            if( p_is_linear ) {
+                *p_is_linear = 0;
+            }
+            return pj_angular_units[i].factor;
+        }
+    }
+    if( p_normalized_name ) {
+        *p_normalized_name = NULL;
+    }
+    if( p_is_linear ) {
+        *p_is_linear = -1;
+    }
+    return 0.0;
+}
 
 /***********************************************************************/
 PJ *CONVERSION(unitconvert,0) {
@@ -390,6 +443,10 @@ PJ *CONVERSION(unitconvert,0) {
     char *s, *name;
     int i;
     double f;
+    int xy_in_is_linear = -1; /* unknown */
+    int xy_out_is_linear = -1; /* unknown */
+    int z_in_is_linear = -1; /* unknown */
+    int z_out_is_linear = -1; /* unknown */
 
     if (0==Q)
         return pj_default_destructor (P, ENOMEM);
@@ -413,11 +470,10 @@ PJ *CONVERSION(unitconvert,0) {
     Q->z_factor  = 1.0;
 
     if ((name = pj_param (P->ctx, P->params, "sxy_in").s) != NULL) {
-        for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i);
-
-        if (s) {
-            f = pj_units[i].factor;
-            proj_log_debug(P, "xy_in unit: %s", pj_units[i].name);
+        const char* normalized_name = NULL;
+        f = get_unit_conversion_factor(name, &xy_in_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "xy_in unit: %s", normalized_name);
         } else {
             if ( (f = pj_param (P->ctx, P->params, "dxy_in").f) == 0.0)
                 return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
@@ -427,11 +483,10 @@ PJ *CONVERSION(unitconvert,0) {
     }
 
     if ((name = pj_param (P->ctx, P->params, "sxy_out").s) != NULL) {
-        for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i);
-
-        if (s) {
-            f = pj_units[i].factor;
-            proj_log_debug(P, "xy_out unit: %s", pj_units[i].name);
+        const char* normalized_name = NULL;
+        f = get_unit_conversion_factor(name, &xy_out_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "xy_out unit: %s", normalized_name);
         } else {
             if ( (f = pj_param (P->ctx, P->params, "dxy_out").f) == 0.0)
                 return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
@@ -440,12 +495,17 @@ PJ *CONVERSION(unitconvert,0) {
             Q->xy_factor /= f;
     }
 
-    if ((name = pj_param (P->ctx, P->params, "sz_in").s) != NULL) {
-        for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i);
+    if( xy_in_is_linear >= 0 && xy_out_is_linear >= 0 &&
+        xy_in_is_linear != xy_out_is_linear ) {
+        proj_log_debug(P, "inconsistent unit type between xy_in and xy_out");
+        return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT);
+    }
 
-        if (s) {
-            f = pj_units[i].factor;
-            proj_log_debug(P, "z_in unit: %s", pj_units[i].name);
+    if ((name = pj_param (P->ctx, P->params, "sz_in").s) != NULL) {
+        const char* normalized_name = NULL;
+        f = get_unit_conversion_factor(name, &z_in_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "z_in unit: %s", normalized_name);
         } else {
             if ( (f = pj_param (P->ctx, P->params, "dz_in").f) == 0.0)
                 return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
@@ -455,11 +515,10 @@ PJ *CONVERSION(unitconvert,0) {
     }
 
     if ((name = pj_param (P->ctx, P->params, "sz_out").s) != NULL) {
-        for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i);
-
-        if (s) {
-            f = pj_units[i].factor;
-            proj_log_debug(P, "z_out unit: %s", pj_units[i].name);
+        const char* normalized_name = NULL;
+        f = get_unit_conversion_factor(name, &z_out_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "z_out unit: %s", normalized_name);
         } else {
             if ( (f = pj_param (P->ctx, P->params, "dz_out").f) == 0.0)
                 return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
@@ -468,6 +527,11 @@ PJ *CONVERSION(unitconvert,0) {
             Q->z_factor /= f;
     }
 
+    if( z_in_is_linear >= 0 && z_out_is_linear >= 0 &&
+        z_in_is_linear != z_out_is_linear ) {
+        proj_log_debug(P, "inconsistent unit type between z_in and z_out");
+        return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT);
+    }
 
     if ((name = pj_param (P->ctx, P->params, "st_in").s) != NULL) {
         for (i = 0; (s = time_units[i].id) && strcmp(name, s) ; ++i);
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index e06ea58b..d37e6d5f 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -67,6 +67,7 @@ pj_err_list[] = {
     "ellipsoidal usage unsupported",                                   /* -56 */
     "only one +init allowed for non-pipeline operations",              /* -57 */
     "argument not numerical or out of range",                          /* -58 */
+    "inconsistent unit type between input and output",                 /* -59 */
 
     /* When adding error messages, remember to update ID defines in
        projects.h, and transient_error array in pj_transform                  */
diff --git a/src/pj_transform.c b/src/pj_transform.c
index e67cc525..b459e4a2 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -69,8 +69,8 @@ static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
 #define M_BF  (defn->datum_params[6])
 
 /*
-** This table is intended to indicate for any given error code in
-** the range 0 to -56, whether that error will occur for all locations (ie.
+** This table is intended to indicate for any given error code
+** whether that error will occur for all locations (ie.
 ** it is a problem with the coordinate system as a whole) in which case the
 ** value would be 0, or if the problem is with the point being transformed
 ** in which case the value is 1.
@@ -98,6 +98,19 @@ static const int transient_error[60] = {
     /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0 };
 
 
+/* -------------------------------------------------------------------- */
+/*      Read transient_error[] in a safe way.                           */
+/* -------------------------------------------------------------------- */
+static int get_transient_error_value(int pos_index)
+{
+    const int array_size =
+        (int)(sizeof(transient_error) / sizeof(transient_error[0]));
+    if( pos_index < 0 || pos_index >= array_size ) {
+        return 0;
+    }
+    return transient_error[pos_index];
+}
+
 
 /* -------------------------------------------------------------------- */
 /*      Transform unusual input coordinate axis orientation to          */
@@ -211,7 +224,7 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *
                         && P->ctx->last_errno != ERANGE)
                     && (P->ctx->last_errno > 0
                         || P->ctx->last_errno < -44 || n == 1
-                        || transient_error[-P->ctx->last_errno] == 0 ) )
+                        || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
                 {
                     return P->ctx->last_errno;
                 }
@@ -249,7 +262,7 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *
                     && P->ctx->last_errno != ERANGE)
                 && (P->ctx->last_errno > 0
                     || P->ctx->last_errno < -44 || n == 1
-                    || transient_error[-P->ctx->last_errno] == 0 ) )
+                    || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
             {
                 return P->ctx->last_errno;
             }
@@ -319,7 +332,7 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *
                         && P->ctx->last_errno != ERANGE)
                     && (P->ctx->last_errno > 0
                         || P->ctx->last_errno < -44 || n == 1
-                        || transient_error[-P->ctx->last_errno] == 0 ) )
+                        || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
                 {
                     return P->ctx->last_errno;
                 }
@@ -358,7 +371,7 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *
                     && P->ctx->last_errno != ERANGE)
                 && (P->ctx->last_errno > 0
                     || P->ctx->last_errno < -44 || n == 1
-                    || transient_error[-P->ctx->last_errno] == 0 ) )
+                    || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
             {
                 return P->ctx->last_errno;
             }
@@ -850,7 +863,7 @@ int pj_datum_transform( PJ *src, PJ *dst,
         z_is_temp = TRUE;
     }
 
-#define CHECK_RETURN(defn) {if( defn->ctx->last_errno != 0 && (defn->ctx->last_errno > 0 || transient_error[-defn->ctx->last_errno] == 0) ) { if( z_is_temp ) pj_dalloc(z); return defn->ctx->last_errno; }}
+#define CHECK_RETURN(defn) {if( defn->ctx->last_errno != 0 && (defn->ctx->last_errno > 0 || get_transient_error_value(-defn->ctx->last_errno) == 0) ) { if( z_is_temp ) pj_dalloc(z); return defn->ctx->last_errno; }}
 
 /* -------------------------------------------------------------------- */
 /*	If this datum requires grid shifts, then apply it to geodetic   */
diff --git a/src/projects.h b/src/projects.h
index 5b923813..2dc073fa 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -583,6 +583,9 @@ enum deprecated_constants_for_now_dropped_analytical_factors {
 #define PJD_ERR_ELLIPSOIDAL_UNSUPPORTED -56
 #define PJD_ERR_TOO_MANY_INITS          -57
 #define PJD_ERR_INVALID_ARG             -58
+#define PJD_ERR_INCONSISTENT_UNIT       -59
+/* NOTE: Remember to update pj_strerrno.c and transient_error in */
+/* pj_transform.c when adding new value */
 
 struct projFileAPI_t;
 
diff --git a/test/gie/unitconvert.gie b/test/gie/unitconvert.gie
index 30d9dc0b..7bbea49d 100644
--- a/test/gie/unitconvert.gie
+++ b/test/gie/unitconvert.gie
@@ -25,4 +25,22 @@ tolerance   0.1
 accept      1       1       1   1
 expect      0.5     0.5     1   1
 
+operation   proj=unitconvert xy_in=deg xy_out=rad
+tolerance   0.000000000001
+accept      1       1       1   1
+expect      0.017453292519943296 0.017453292519943296 1 1
+
+operation   proj=unitconvert xy_in=grad xy_out=deg
+tolerance   0.000000000001
+accept      50 50 1 1
+expect      45 45 1 1
+
+operation   proj=unitconvert xy_in=m xy_out=rad
+accept      1 1 1 1
+expect      failure
+
+operation   proj=unitconvert z_in=rad z_out=m
+accept      1 1 1 1
+expect      failure
+
 </gie>