Custom unit factors in unitconvert.

Added the possibility to use custom unit factors. Similar to the
classic +to_meter parameter the conversion factor is related to meters,
i.e. the factor for conversion from kilometers meters is 1000. The
custom unit factors is given using the existing xy_in, xy_out, z_in and
z_out parameters, for example:

    proj=unitconvert xy_in=4.5 +xy_out=mm

3 files changed, 84 insertions, 71 deletions

diff --git a/src/PJ_unitconvert.c b/src/PJ_unitconvert.c
index 9cd2acf8..bf1fcd37 100644
--- a/src/PJ_unitconvert.c
+++ b/src/PJ_unitconvert.c
@@ -80,6 +80,13 @@ struct TIME_UNITS {
     char        *name;      /* comments */
 };
 
+struct pj_opaque_unitconvert {
+    int     t_in_id;        /* time unit id for the time input unit   */
+    int     t_out_id;       /* time unit id for the time output unit  */
+    double  xy_factor;      /* unit conversion factor for horizontal components */
+    double  z_factor;       /* unit conversion factor for vertical components */
+};
+
 
 /***********************************************************************/
 static int is_leap_year(int year) {
@@ -255,7 +262,6 @@ static double mjd_to_yyyymmdd(double mjd) {
     return year*10000.0 + month*100.0 + day;
 }
 
-
 static const struct TIME_UNITS time_units[] = {
     {"mjd",         mjd_to_mjd,         mjd_to_mjd,         "Modified julian date"},
     {"decimalyear", decimalyear_to_mjd, mjd_to_decimalyear, "Decimal year"},
@@ -264,15 +270,6 @@ static const struct TIME_UNITS time_units[] = {
     {NULL,          NULL,               NULL,               NULL}
 };
 
-struct pj_opaque_unitconvert {
-    int     t_in_id;    /* time unit id for the time input unit   */
-    int     t_out_id;   /* time unit id for the time output unit  */
-    int     xy_in_id;   /* unit id for the horizontal input unit  */
-    int     xy_out_id;  /* unit id for the horizontal output unit */
-    int     z_in_id;    /* unit id for the vertical input unit    */
-    int     z_out_id;   /* unit id for the vertical output unit   */
-};
-
 
 /***********************************************************************/
 static XY forward_2d(LP lp, PJ *P) {
@@ -283,8 +280,8 @@ static XY forward_2d(LP lp, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.lp = lp;
 
-    point.xy.x *= pj_units[Q->xy_in_id].factor / pj_units[Q->xy_out_id].factor;
-    point.xy.y *= pj_units[Q->xy_in_id].factor / pj_units[Q->xy_out_id].factor;
+    point.xy.x *= Q->xy_factor;
+    point.xy.y *= Q->xy_factor;
 
     return point.xy;
 }
@@ -299,8 +296,8 @@ static LP reverse_2d(XY xy, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.xy = xy;
 
-    point.xy.x *= pj_units[Q->xy_out_id].factor / pj_units[Q->xy_in_id].factor;
-    point.xy.y *= pj_units[Q->xy_out_id].factor / pj_units[Q->xy_in_id].factor;
+    point.xy.x /= Q->xy_factor;
+    point.xy.y /= Q->xy_factor;
 
     return point.lp;
 }
@@ -318,7 +315,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
     /* take care of the horizontal components in the 2D function */
     point.xy = forward_2d(point.lp, P);
 
-    point.xyz.z *= pj_units[Q->z_in_id].factor / pj_units[Q->z_out_id].factor;
+    point.xyz.z *= Q->z_factor;
 
     return point.xyz;
 }
@@ -335,7 +332,7 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
     /* take care of the horizontal components in the 2D function */
     point.lp = reverse_2d(point.xy, P);
 
-    point.xyz.z *= pj_units[Q->z_out_id].factor / pj_units[Q->z_in_id].factor;
+    point.xyz.z /= Q->z_factor;
 
     return point.lpz;
 }
@@ -387,6 +384,7 @@ PJ *CONVERSION(unitconvert,0) {
     struct pj_opaque_unitconvert *Q = pj_calloc (1, sizeof (struct pj_opaque_unitconvert));
     char *s, *name;
     int i;
+    double f;
 
     if (0==Q)
         return pj_default_destructor (P, ENOMEM);
@@ -406,40 +404,63 @@ PJ *CONVERSION(unitconvert,0) {
     Q->t_in_id = -1;
     Q->t_out_id = -1;
 
+    Q->xy_factor = 1.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) return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID);
-
-        Q->xy_in_id = i;
-        proj_log_debug(P, "xy_in unit: %s", pj_units[i].name);
+        if (s) {
+            f = pj_units[i].factor;
+            proj_log_debug(P, "xy_in unit: %s", pj_units[i].name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dxy_in").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID);
+        }
+        if (f != 0.0)
+            Q->xy_factor *= f;
     }
 
     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) return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID);
-
-        Q->xy_out_id = i;
-        proj_log_debug(P, "xy_out unit: %s", pj_units[i].name);
+        if (s) {
+            f = pj_units[i].factor;
+            proj_log_debug(P, "xy_out unit: %s", pj_units[i].name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dxy_out").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID);
+        }
+        if (f != 0.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 (!s) return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID); /* unknown unit conversion id */
-
-        Q->z_in_id = i;
-        proj_log_debug(P, "z_in unit: %s", pj_units[i].name);
+        if (s) {
+            f = pj_units[i].factor;
+            proj_log_debug(P, "z_in unit: %s", pj_units[i].name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dz_in").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID);
+        }
+        if (f != 0.0)
+            Q->z_factor *= f;
     }
 
     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) return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID); /* unknown unit conversion id */
-
-        Q->z_out_id = i;
-        proj_log_debug(P, "z_out unit: %s", pj_units[i].name);
+        if (s) {
+            f = pj_units[i].factor;
+            proj_log_debug(P, "z_out unit: %s", pj_units[i].name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dz_out").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOW_UNIT_ID);
+        }
+        if (f != 0.0)
+            Q->z_factor /= f;
     }
 
 
diff --git a/src/gie.c b/src/gie.c
index c6440a2d..e37acf62 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -1850,34 +1850,6 @@ static int test_time(const char* args, double tol, double t_in, double t_exp) {
     return ret;
 }
 
-static int test_xyz(const char* args, double tol, PJ_COORD in, PJ_COORD exp) {
-    PJ_COORD out = {{0,0,0,0}}, obs_in = {{0,0,0,0}};
-    PJ *P = proj_create(PJ_DEFAULT_CTX, args);
-    int ret = 0;
-
-    if (P == 0)
-        return 5;
-
-    obs_in.xyz = in.xyz;
-    out = proj_trans(P, PJ_FWD, obs_in);
-    if (proj_xyz_dist(out.xyz, exp.xyz) > tol) {
-        printf("exp: %10.10g, %10.10g, %10.10g\n", exp.xyz.x, exp.xyz.y, exp.xyz.z);
-        printf("out: %10.10g, %10.10g, %10.10g\n", out.xyz.x, out.xyz.y, out.xyz.z);
-        ret = 1;
-    }
-
-    out = proj_trans(P, PJ_INV, out);
-    if (proj_xyz_dist(out.xyz, in.xyz) > tol) {
-        printf("exp: %g, %g, %g\n", in.xyz.x, in.xyz.y, in.xyz.z);
-        printf("out: %g, %g, %g\n", out.xyz.x, out.xyz.y, out.xyz.z);
-        ret += 2;
-    }
-    proj_destroy(P);
-    proj_log_level(NULL, 0);
-    return ret;
-}
-
-
 static int pj_unitconvert_selftest (void) {
     int ret = 0;
     char args1[] = "+proj=unitconvert +t_in=decimalyear +t_out=decimalyear";
@@ -1892,22 +1864,14 @@ static int pj_unitconvert_selftest (void) {
     char args4[] = "+proj=unitconvert +t_in=gps_week +t_out=decimalyear";
     double in4 = 1877.71428, exp4 = 2016.0;
 
-    char args5[] = "+proj=unitconvert +xy_in=m +xy_out=dm +z_in=cm +z_out=mm";
-    PJ_COORD in5 = {{55.25, 23.23, 45.5, 0}}, exp5 = {{552.5, 232.3, 455.0, 0}};
-
-    char args6[] = "+proj=unitconvert +xy_in=m +xy_out=m +z_in=m +z_out=m";
-    PJ_COORD in6 = {{12.3, 45.6, 7.89, 0}};
-
-    char args7[] = "+proj=unitconvert +t_in=yyyymmdd +t_out=yyyymmdd";
-    double in7 = 20170131;
+    char args5[] = "+proj=unitconvert +t_in=yyyymmdd +t_out=yyyymmdd";
+    double in5 = 20170131;
 
     ret = test_time(args1, 1e-6, in1, in1);   if (ret) return ret + 10;
     ret = test_time(args2, 1e-6, in2, in2);   if (ret) return ret + 20;
     ret = test_time(args3, 1e-6, in3, in3);   if (ret) return ret + 30;
     ret = test_time(args4, 1e-6, in4, exp4);  if (ret) return ret + 40;
-    ret = test_xyz (args5, 1e-10, in5, exp5); if (ret) return ret + 50;
-    ret = test_xyz (args6, 1e-10, in6, in6);  if (ret) return ret + 60;
-    ret = test_time(args7, 1e-6, in7, in7);   if (ret) return ret + 70;
+    ret = test_time(args5, 1e-6, in5, in5);   if (ret) return ret + 50;
 
     return 0;
 
diff --git a/test/gie/unitconvert.gie b/test/gie/unitconvert.gie
new file mode 100644
index 00000000..30d9dc0b
--- /dev/null
+++ b/test/gie/unitconvert.gie
@@ -0,0 +1,28 @@
+-------------------------------------------------------------------------------
+                        Tests for the unitconvert operation
+-------------------------------------------------------------------------------
+
+<gie>
+
+operation   proj=unitconvert xy_in=m xy_out=dm z_in=cm z_out=mm
+tolerance   0.1
+accept       55.25   23.23   45.5
+expect      552.5   232.3   455.0
+
+operation   proj=unitconvert +xy_in=m +xy_out=m +z_in=m +z_out=m
+tolerance   0.1
+accept      12.3    45.6    7.89
+expect      12.3    45.6    7.89
+
+operation   proj=unitconvert xy_in=dm xy_out=dm
+tolerance   0.1
+accept      1 1 1 1
+expect      1 1 1 1
+
+
+operation   proj=unitconvert xy_in=2.0 xy_out=4.0
+tolerance   0.1
+accept      1       1       1   1
+expect      0.5     0.5     1   1
+
+</gie>