added preliminary support for +axis (#18)

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1825 4e78687f-474d-0410-85f9-8d5e500ac6b2

4 files changed, 171 insertions, 4 deletions

diff --git a/src/pj_init.c b/src/pj_init.c
index f40d2ad9..eb8cf268 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -270,6 +270,7 @@ pj_init(int argc, char **argv) {
         PIN->is_geocent = 0;
         PIN->is_long_wrap_set = 0;
         PIN->long_wrap_center = 0.0;
+        strcpy( PIN->axis, "enu" );
 
         /* set datum parameters */
         if (pj_datum_set(start, PIN)) goto bum_call;
@@ -308,6 +309,33 @@ pj_init(int argc, char **argv) {
  	if (PIN->is_long_wrap_set)
  		PIN->long_wrap_center = pj_param(start, "rlon_wrap").f;
 
+	/* axis orientation */
+        if( (pj_param(start,"saxis").s) != NULL )
+        {
+            static const char *axis_legal = "ewnsud";
+            const char *axis_arg = pj_param(start,"saxis").s;
+            if( strlen(axis_arg) != 3 )
+            {
+                pj_errno = PJD_ERR_AXIS;
+                goto bum_call;
+            }
+
+            if( strchr( axis_legal, axis_arg[0] ) == NULL
+                || strchr( axis_legal, axis_arg[1] ) == NULL
+                || (axis_arg[2] && strchr( axis_legal, axis_arg[1] ) == NULL))
+            {
+                pj_errno = PJD_ERR_AXIS;
+                goto bum_call;
+            }
+
+            /* it would be nice to validate we don't have on axis repeated */
+            strcpy( PIN->axis, axis_arg );
+        }
+
+ 	PIN->is_long_wrap_set = pj_param(start, "tlon_wrap").i;
+ 	if (PIN->is_long_wrap_set)
+ 		PIN->long_wrap_center = pj_param(start, "rlon_wrap").f;
+
 	/* central meridian */
 	PIN->lam0=pj_param(start, "rlon_0").f;
 
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index 4db1fb9a..320aaafc 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -52,6 +52,7 @@ pj_err_list[] = {
 	"unparseable coordinate system definition",	/* -44 */
 	"geocentric transformation missing z or ellps",	/* -45 */
 	"unknown prime meridian conversion id",		/* -46 */
+	"illegal axis orientation combination",		/* -47 */
 };
 	char *
 pj_strerrno(int err) 
diff --git a/src/pj_transform.c b/src/pj_transform.c
index c64b410f..ba904db2 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -36,6 +36,10 @@
 
 PJ_CVSID("$Id$");
 
+static int pj_adjust_axis( const char *axis, int denormalize_flag,
+                           long point_count, int point_offset, 
+                           double *x, double *y, double *z );
+
 #ifndef SRS_WGS84_SEMIMAJOR
 #define SRS_WGS84_SEMIMAJOR 6378137.0
 #endif
@@ -63,13 +67,13 @@ PJ_CVSID("$Id$");
 ** list or something, but while experimenting with it this should be fine. 
 */
 
-static const int transient_error[45] = {
+static const int transient_error[50] = {
     /*             0  1  2  3  4  5  6  7  8  9   */
     /* 0 to 9 */   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,  
     /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
     /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
-    /* 40 to 44 */ 0, 0, 0, 0, 0 };
+    /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
 
 /************************************************************************/
 /*                            pj_transform()                            */
@@ -93,6 +97,20 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
         point_offset = 1;
 
 /* -------------------------------------------------------------------- */
+/*      Transform unusual input coordinate axis orientation to          */
+/*      standard form if needed.                                        */
+/* -------------------------------------------------------------------- */
+    if( strcmp(srcdefn->axis,"enu") != 0 )
+    {
+        int err;
+
+        err = pj_adjust_axis( srcdefn->axis, 0, point_count, point_offset, 
+                              x, y, z );
+        if( err != 0 )
+            return err;
+    }
+
+/* -------------------------------------------------------------------- */
 /*      Transform geocentric source coordinates to lat/long.            */
 /* -------------------------------------------------------------------- */
     if( srcdefn->is_geocent )
@@ -282,6 +300,20 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
         }
     }
 
+/* -------------------------------------------------------------------- */
+/*      Transform normalized axes into unusual output coordinate axis   */
+/*      orientation if needed.                                          */
+/* -------------------------------------------------------------------- */
+    if( strcmp(dstdefn->axis,"enu") != 0 )
+    {
+        int err;
+
+        err = pj_adjust_axis( dstdefn->axis, 1, point_count, point_offset, 
+                              x, y, z );
+        if( err != 0 )
+            return err;
+    }
+
     return 0;
 }
 
@@ -644,3 +676,107 @@ int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
     return 0;
 }
 
+/************************************************************************/
+/*                           pj_adjust_axis()                           */
+/*                                                                      */
+/*      Normalize or de-normalized the x/y/z axes.  The normal form     */
+/*      is "enu" (easting, northing, up).                               */
+/************************************************************************/
+static int pj_adjust_axis( const char *axis, int denormalize_flag,
+                           long point_count, int point_offset, 
+                           double *x, double *y, double *z )
+
+{
+    double x_in, y_in, z_in = 0.0;
+    int i, i_axis;
+
+    if( !denormalize_flag )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            x_in = x[point_offset*i];
+            y_in = y[point_offset*i];
+            if( z )
+                z_in = z[point_offset*i];
+     
+            for( i_axis = 0; i_axis < 3; i_axis++ )
+            {
+                double value;
+
+                if( i_axis == 0 )
+                    value = x_in;
+                else if( i_axis == 1 )
+                    value = y_in;
+                else
+                    value = z_in;
+                
+                switch( axis[i_axis] )
+                {
+                  case 'e':
+                    x[point_offset*i] = value; break;
+                  case 'w':
+                    x[point_offset*i] = -value; break;
+                  case 'n':
+                    y[point_offset*i] = value; break;
+                  case 's':
+                    y[point_offset*i] = -value; break;
+                  case 'u':
+                    if( z ) z[point_offset*i] = value; break;
+                  case 'd':
+                    if( z ) z[point_offset*i] = -value; break;
+                  default:
+                    pj_errno = PJD_ERR_AXIS;
+                    return PJD_ERR_AXIS;
+                }
+            } /* i_axis */
+        } /* i (point) */
+    }
+
+    else /* denormalize */
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            x_in = x[point_offset*i];
+            y_in = y[point_offset*i];
+            if( z )
+                z_in = z[point_offset*i];
+     
+            for( i_axis = 0; i_axis < 3; i_axis++ )
+            {
+                double *target;
+
+                if( i_axis == 2 && z == NULL )
+                    continue;
+
+                if( i_axis == 0 )
+                    target = x;
+                else if( i_axis == 1 )
+                    target = y;
+                else
+                    target = z;
+                
+                switch( axis[i_axis] )
+                {
+                  case 'e':
+                    target[point_offset*i] = x_in; break;
+                  case 'w':
+                    target[point_offset*i] = -x_in; break;
+                  case 'n':
+                    target[point_offset*i] = y_in; break;
+                  case 's':
+                    target[point_offset*i] = -y_in; break;
+                  case 'u':
+                    target[point_offset*i] = z_in; break;
+                  case 'd':
+                    target[point_offset*i] = -z_in; break;
+                  default:
+                    pj_errno = PJD_ERR_AXIS;
+                    return PJD_ERR_AXIS;
+                }
+            } /* i_axis */
+        } /* i (point) */
+    }
+    
+    return 0;
+}
+
diff --git a/src/projects.h b/src/projects.h
index 7fec1ce5..a526cfe3 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -131,8 +131,9 @@ extern double hypot(double, double);
 #define PJD_GRIDSHIFT 3
 #define PJD_WGS84     4   /* WGS84 (or anything considered equivelent) */
 
-/* datum system errors */
-#define PJD_ERR_GEOCENTRIC -45
+/* library errors */
+#define PJD_ERR_GEOCENTRIC          -45
+#define PJD_ERR_AXIS                -47
 
 #define USE_PROJUV 
 
@@ -235,6 +236,7 @@ typedef struct PJconsts {
         double  from_greenwich; /* prime meridian offset (in radians) */
         double  long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/
         int     is_long_wrap_set;
+        char    axis[4];
         
 #ifdef PROJ_PARMS__
 PROJ_PARMS__