From 3405bf76e70d12df6b9dc0c9c129c7dc96918e8d Mon Sep 17 00:00:00 2001 From: Frank Warmerdam Date: Mon, 1 Mar 2010 04:32:18 +0000 Subject: added preliminary support for +axis (#18) git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1825 4e78687f-474d-0410-85f9-8d5e500ac6b2 --- src/pj_init.c | 28 +++++++++++ src/pj_strerrno.c | 1 + src/pj_transform.c | 140 ++++++++++++++++++++++++++++++++++++++++++++++++++++- src/projects.h | 6 ++- 4 files changed, 171 insertions(+), 4 deletions(-) (limited to 'src') 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; @@ -304,6 +305,33 @@ pj_init(int argc, char **argv) { PIN->over = pj_param(start, "bover").i; /* longitude center for wrapping */ + 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; + + /* 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; 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() */ @@ -92,6 +96,20 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset, if( point_offset == 0 ) 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. */ /* -------------------------------------------------------------------- */ @@ -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__ -- cgit v1.2.3