diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/PJ_qsc.c | 2 | ||||
| -rw-r--r-- | src/pj_inv.c | 4 | ||||
| -rw-r--r-- | src/pj_transform.c | 684 | ||||
| -rw-r--r-- | src/proj.h | 4 | ||||
| -rw-r--r-- | src/proj_4D_api.c | 30 | ||||
| -rw-r--r-- | src/proj_internal.h | 4 |
6 files changed, 410 insertions, 318 deletions
diff --git a/src/PJ_qsc.c b/src/PJ_qsc.c index 36dabeb9..0f5d2299 100644 --- a/src/PJ_qsc.c +++ b/src/PJ_qsc.c @@ -14,7 +14,7 @@ * is described in * [LK12] * M. Lambers and A. Kolb, "Ellipsoidal Cube Maps for Accurate Rendering of - * Planetary-Scale Terrain Data", Proc. Pacfic Graphics (Short Papers), Sep. + * Planetary-Scale Terrain Data", Proc. Pacific Graphics (Short Papers), Sep. * 2012 * * You have to choose one of the following projection centers, diff --git a/src/pj_inv.c b/src/pj_inv.c index 4ea88b69..9918c59d 100644 --- a/src/pj_inv.c +++ b/src/pj_inv.c @@ -115,9 +115,9 @@ static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) { return coo; /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */ - /* Multiplying by ra, rather than dividing by a because the CALCOFI projection */ + /* Multiplying by ra, rather than dividing by a because the CalCOFI projection */ /* stomps on a and hence (apparently) depends on this to roundtrip correctly */ - /* (CALCOFI avoids further scaling by stomping - but a better solution is possible) */ + /* (CalCOFI avoids further scaling by stomping - but a better solution is possible) */ coo.xyz.x *= P->ra; coo.xyz.y *= P->ra; return coo; diff --git a/src/pj_transform.c b/src/pj_transform.c index c72df2e7..06a2dfb0 100644 --- a/src/pj_transform.c +++ b/src/pj_transform.c @@ -27,11 +27,23 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ -#include <projects.h> +#include "projects.h" #include <string.h> #include <math.h> #include "geocent.h" + +/* Apply transformation to observation - in forward or inverse direction */ +/* Copied from proj.h */ +enum PJ_DIRECTION { + PJ_FWD = 1, /* Forward */ + PJ_IDENT = 0, /* Do nothing */ + PJ_INV = -1 /* Inverse */ +}; +typedef enum PJ_DIRECTION PJ_DIRECTION; + + + static int pj_adjust_axis( projCtx ctx, const char *axis, int denormalize_flag, long point_count, int point_offset, double *x, double *y, double *z ); @@ -81,388 +93,470 @@ static const int transient_error[60] = { /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0 }; -/************************************************************************/ -/* pj_transform() */ -/* */ -/* Currently this function doesn't recognise if two projections */ -/* are identical (to short circuit reprojection) because it is */ -/* difficult to compare PJ structures (since there are some */ -/* projection specific components). */ -/************************************************************************/ - -int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - long i; - srcdefn->ctx->last_errno = 0; - dstdefn->ctx->last_errno = 0; - - 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; +static int adjust_axes (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) { + /* Nothing to do? */ + if (0==strcmp(P->axis,"enu")) + return 0; + + return pj_adjust_axis( P->ctx, P->axis, + dir==PJ_FWD ? 1: 0, n, dist, x, y, z ); +} + - err = pj_adjust_axis( srcdefn->ctx, srcdefn->axis, - 0, point_count, point_offset, x, y, z ); - if( err != 0 ) - return err; + +/* ----------------------------------------------------------------------- */ +/* Transform cartesian ("geocentric") source coordinates to lat/long, */ +/* if needed */ +/* ----------------------------------------------------------------------- */ +static int geographic_to_cartesian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) { + int res; + long i; + double fac = P->to_meter; + + /* Nothing to do? */ + if (!P->is_geocent) + return 0; + + if ( z == NULL ) { + pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC); + return PJD_ERR_GEOCENTRIC; + } + + if (PJ_FWD==dir) { + fac = P->fr_meter; + res = pj_geodetic_to_geocentric( P->a_orig, P->es_orig, n, dist, x, y, z ); + if (res) + return res; + } + + if (fac != 1.0) { + for( i = 0; i < n; i++ ) { + if( x[dist*i] != HUGE_VAL ) { + x[dist*i] *= fac; + y[dist*i] *= fac; + z[dist*i] *= fac; + } + } } + if (PJ_FWD==dir) + return 0; + return pj_geocentric_to_geodetic( + P->a_orig, P->es_orig, + n, dist, + x, y, z + ); +} + + + + + + + + + + /* -------------------------------------------------------------------- */ -/* Transform geocentric source coordinates to lat/long. */ +/* Transform destination points to projection coordinates, if */ +/* desired. */ +/* */ +/* Ought to fold this into projected_to_geographic */ /* -------------------------------------------------------------------- */ - if( srcdefn->is_geocent ) +static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *y, double *z) { + long i; + + /* Nothing to do? */ + if (P->is_latlong) + return 0; + if (P->is_geocent) + return 0; + + if(P->fwd3d != NULL) { - int err; - if( z == NULL ) + /* Three dimensions must be defined */ + if ( z == NULL) { - pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC); + pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC); return PJD_ERR_GEOCENTRIC; } - if( srcdefn->to_meter != 1.0 ) + for( i = 0; i < n; i++ ) { - for( i = 0; i < point_count; i++ ) + XYZ projected_loc; + LPZ geodetic_loc; + + geodetic_loc.u = x[dist*i]; + geodetic_loc.v = y[dist*i]; + geodetic_loc.w = z[dist*i]; + + if (geodetic_loc.u == HUGE_VAL) + continue; + + projected_loc = pj_fwd3d( geodetic_loc, P); + if( P->ctx->last_errno != 0 ) { - if( x[point_offset*i] != HUGE_VAL ) + if( (P->ctx->last_errno != 33 /*EDOM*/ + && P->ctx->last_errno != 34 /*ERANGE*/ ) + && (P->ctx->last_errno > 0 + || P->ctx->last_errno < -44 || n == 1 + || transient_error[-P->ctx->last_errno] == 0 ) ) + return P->ctx->last_errno; + else { - x[point_offset*i] *= srcdefn->to_meter; - y[point_offset*i] *= srcdefn->to_meter; - z[point_offset*i] *= srcdefn->to_meter; + projected_loc.u = HUGE_VAL; + projected_loc.v = HUGE_VAL; + projected_loc.w = HUGE_VAL; } } - } - err = pj_geocentric_to_geodetic( srcdefn->a_orig, srcdefn->es_orig, - point_count, point_offset, - x, y, z ); - if( err != 0 ) - return err; + x[dist*i] = projected_loc.u; + y[dist*i] = projected_loc.v; + z[dist*i] = projected_loc.w; + } + return 0; } -/* -------------------------------------------------------------------- */ -/* Transform source points to lat/long, if they aren't */ -/* already. */ -/* -------------------------------------------------------------------- */ - else if( !srcdefn->is_latlong ) + for( i = 0; i <n; i++ ) { + XY projected_loc; + LP geodetic_loc; - /* Check first if projection is invertible. */ - if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL)) - { - pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 ); - pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR, - "pj_transform(): source projection not invertable" ); - return -17; - } + geodetic_loc.u = x[dist*i]; + geodetic_loc.v = y[dist*i]; + + if( geodetic_loc.u == HUGE_VAL ) + continue; - /* If invertible - First try inv3d if defined */ - if (srcdefn->inv3d != NULL) + projected_loc = pj_fwd( geodetic_loc, P ); + if( P->ctx->last_errno != 0 ) { - /* Three dimensions must be defined */ - if ( z == NULL) + if( (P->ctx->last_errno != 33 /*EDOM*/ + && P->ctx->last_errno != 34 /*ERANGE*/ ) + && (P->ctx->last_errno > 0 + || P->ctx->last_errno < -44 || n == 1 + || transient_error[-P->ctx->last_errno] == 0 ) ) + return P->ctx->last_errno; + else { - pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC); - return PJD_ERR_GEOCENTRIC; + projected_loc.u = HUGE_VAL; + projected_loc.v = HUGE_VAL; } + } - for (i=0; i < point_count; i++) - { - XYZ projected_loc; - XYZ geodetic_loc; + x[dist*i] = projected_loc.u; + y[dist*i] = projected_loc.v; + } + return 0; +} - projected_loc.u = x[point_offset*i]; - projected_loc.v = y[point_offset*i]; - projected_loc.w = z[point_offset*i]; - if (projected_loc.u == HUGE_VAL) - continue; - geodetic_loc = pj_inv3d(projected_loc, srcdefn); - if( srcdefn->ctx->last_errno != 0 ) - { - if( (srcdefn->ctx->last_errno != 33 /*EDOM*/ - && srcdefn->ctx->last_errno != 34 /*ERANGE*/ ) - && (srcdefn->ctx->last_errno > 0 - || srcdefn->ctx->last_errno < -44 || point_count == 1 - || transient_error[-srcdefn->ctx->last_errno] == 0 ) ) - return srcdefn->ctx->last_errno; - else - { - geodetic_loc.u = HUGE_VAL; - geodetic_loc.v = HUGE_VAL; - geodetic_loc.w = HUGE_VAL; - } - } - x[point_offset*i] = geodetic_loc.u; - y[point_offset*i] = geodetic_loc.v; - z[point_offset*i] = geodetic_loc.w; - } +/* ----------------------------------------------------------------------- */ +/* Transform projected source coordinates to lat/long, if needed */ +/* ----------------------------------------------------------------------- */ +static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *y, double *z) { + long i; + + /* Nothing to do? */ + if (P->is_latlong) + return 0; + + /* Check first if projection is invertible. */ + if( (P->inv3d == NULL) && (P->inv == NULL)) + { + pj_ctx_set_errno( pj_get_ctx(P), -17 ); + pj_log( pj_get_ctx(P), PJ_LOG_ERROR, + "pj_transform(): source projection not invertable" ); + return -17; + } + /* If invertible - First try inv3d if defined */ + if (P->inv3d != NULL) + { + /* Three dimensions must be defined */ + if ( z == NULL) + { + pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC); + return PJD_ERR_GEOCENTRIC; } - else + + for (i=0; i < n; i++) { - /* Fallback to the original PROJ.4 API 2d inversion - inv */ - for( i = 0; i < point_count; i++ ) - { - XY projected_loc; - LP geodetic_loc; + XYZ projected_loc; + XYZ geodetic_loc; - projected_loc.u = x[point_offset*i]; - projected_loc.v = y[point_offset*i]; + projected_loc.u = x[dist*i]; + projected_loc.v = y[dist*i]; + projected_loc.w = z[dist*i]; - if( projected_loc.u == HUGE_VAL ) - continue; + if (projected_loc.u == HUGE_VAL) + continue; - geodetic_loc = pj_inv( projected_loc, srcdefn ); - if( srcdefn->ctx->last_errno != 0 ) + geodetic_loc = pj_inv3d(projected_loc, P); + if( P->ctx->last_errno != 0 ) + { + if( (P->ctx->last_errno != 33 /*EDOM*/ + && P->ctx->last_errno != 34 /*ERANGE*/ ) + && (P->ctx->last_errno > 0 + || P->ctx->last_errno < -44 || n == 1 + || transient_error[-P->ctx->last_errno] == 0 ) ) + return P->ctx->last_errno; + else { - if( (srcdefn->ctx->last_errno != 33 /*EDOM*/ - && srcdefn->ctx->last_errno != 34 /*ERANGE*/ ) - && (srcdefn->ctx->last_errno > 0 - || srcdefn->ctx->last_errno < -44 || point_count == 1 - || transient_error[-srcdefn->ctx->last_errno] == 0 ) ) - return srcdefn->ctx->last_errno; - else - { - geodetic_loc.u = HUGE_VAL; - geodetic_loc.v = HUGE_VAL; - } + geodetic_loc.u = HUGE_VAL; + geodetic_loc.v = HUGE_VAL; + geodetic_loc.w = HUGE_VAL; } - - x[point_offset*i] = geodetic_loc.u; - y[point_offset*i] = geodetic_loc.v; } + + x[dist*i] = geodetic_loc.u; + y[dist*i] = geodetic_loc.v; + z[dist*i] = geodetic_loc.w; + } + return 0; } -/* -------------------------------------------------------------------- */ -/* But if they are already lat long, adjust for the prime */ -/* meridian if there is one in effect. */ -/* -------------------------------------------------------------------- */ - if ((srcdefn->is_geocent || srcdefn->is_latlong) && ( srcdefn->from_greenwich != 0.0 )) - { - for( i = 0; i < point_count; i++ ) + /* Fallback to the original PROJ.4 API 2d inversion - inv */ + for( i = 0; i < n; i++ ) { + XY projected_loc; + LP geodetic_loc; + + projected_loc.u = x[dist*i]; + projected_loc.v = y[dist*i]; + + if( projected_loc.u == HUGE_VAL ) + continue; + + geodetic_loc = pj_inv( projected_loc, P ); + if( P->ctx->last_errno != 0 ) { - if( x[point_offset*i] != HUGE_VAL ) - x[point_offset*i] += srcdefn->from_greenwich; + if( (P->ctx->last_errno != 33 /*EDOM*/ + && P->ctx->last_errno != 34 /*ERANGE*/ ) + && (P->ctx->last_errno > 0 + || P->ctx->last_errno < -44 || n == 1 + || transient_error[-P->ctx->last_errno] == 0 ) ) + return P->ctx->last_errno; + else + { + geodetic_loc.u = HUGE_VAL; + geodetic_loc.v = HUGE_VAL; + } } - } -/* -------------------------------------------------------------------- */ -/* Do we need to translate from geoid to ellipsoidal vertical */ -/* datum? */ -/* -------------------------------------------------------------------- */ - if( srcdefn->has_geoid_vgrids && z != NULL ) - { - if( pj_apply_vgridshift( srcdefn, "sgeoidgrids", - &(srcdefn->vgridlist_geoid), - &(srcdefn->vgridlist_geoid_count), - 0, point_count, point_offset, x, y, z ) != 0 ) - return pj_ctx_get_errno(srcdefn->ctx); + x[dist*i] = geodetic_loc.u; + y[dist*i] = geodetic_loc.v; } + return 0; +} -/* -------------------------------------------------------------------- */ -/* Convert datums if needed, and possible. */ -/* -------------------------------------------------------------------- */ - if( pj_datum_transform( srcdefn, dstdefn, point_count, point_offset, - x, y, z ) != 0 ) - { - if( srcdefn->ctx->last_errno != 0 ) - return srcdefn->ctx->last_errno; - else - return dstdefn->ctx->last_errno; - } -/* -------------------------------------------------------------------- */ -/* Do we need to translate from ellipsoidal to geoid vertical */ -/* datum? */ -/* -------------------------------------------------------------------- */ - if( dstdefn->has_geoid_vgrids && z != NULL ) - { - if( pj_apply_vgridshift( dstdefn, "sgeoidgrids", - &(dstdefn->vgridlist_geoid), - &(dstdefn->vgridlist_geoid_count), - 1, point_count, point_offset, x, y, z ) != 0 ) - return dstdefn->ctx->last_errno; - } /* -------------------------------------------------------------------- */ -/* But if they are staying lat long, adjust for the prime */ -/* meridian if there is one in effect. */ +/* Adjust for the prime meridian if needed. */ /* -------------------------------------------------------------------- */ - if ((dstdefn->is_geocent || dstdefn->is_latlong) && ( dstdefn->from_greenwich != 0.0 )) - { - for( i = 0; i < point_count; i++ ) - { - if( x[point_offset*i] != HUGE_VAL ) - x[point_offset*i] -= dstdefn->from_greenwich; - } +static int prime_meridian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x) { + int i; + double pm = P->from_greenwich; + + /* Nothing to do? */ + if (pm==0.0) + return 0; + if (!(P->is_geocent || P->is_latlong)) + return 0; + + if (dir==PJ_FWD) + pm = -pm; + + for (i = 0; i < n; i++) + if (x[dist*i] != HUGE_VAL) + x[dist*i] += pm; + + return 0; } + + /* -------------------------------------------------------------------- */ -/* Transform destination latlong to geocentric if required. */ +/* Adjust for vertical scale factor if needed */ /* -------------------------------------------------------------------- */ - if( dstdefn->is_geocent ) - { - if( z == NULL ) - { - pj_ctx_set_errno( dstdefn->ctx, PJD_ERR_GEOCENTRIC ); - return PJD_ERR_GEOCENTRIC; - } +static int height_unit (PJ *P, PJ_DIRECTION dir, long n, int dist, double *z) { + int i; + double fac = P->vto_meter; - pj_geodetic_to_geocentric( dstdefn->a_orig, dstdefn->es_orig, - point_count, point_offset, x, y, z ); + if (PJ_FWD==dir) + fac = P->vfr_meter; - if( dstdefn->fr_meter != 1.0 ) - { - for( i = 0; i < point_count; i++ ) - { - if( x[point_offset*i] != HUGE_VAL ) - { - x[point_offset*i] *= dstdefn->fr_meter; - y[point_offset*i] *= dstdefn->fr_meter; - z[point_offset*i] *= srcdefn->fr_meter; - } - } - } - } + /* Nothing to do? */ + if (fac==1.0) + return 0; + if (0==z) + return 0; -/* -------------------------------------------------------------------- */ -/* Transform destination points to projection coordinates, if */ -/* desired. */ -/* -------------------------------------------------------------------- */ - else if( !dstdefn->is_latlong ) - { + for (i = 0; i < n; i++) + if (z[dist*i] != HUGE_VAL ) + z[dist*i] *= fac; - if( dstdefn->fwd3d != NULL) - { - /* Three dimensions must be defined */ - if ( z == NULL) - { - pj_ctx_set_errno( pj_get_ctx(dstdefn), PJD_ERR_GEOCENTRIC); - return PJD_ERR_GEOCENTRIC; - } + return 0; +} - for( i = 0; i < point_count; i++ ) - { - XYZ projected_loc; - LPZ geodetic_loc; - geodetic_loc.u = x[point_offset*i]; - geodetic_loc.v = y[point_offset*i]; - geodetic_loc.w = z[point_offset*i]; - if (geodetic_loc.u == HUGE_VAL) - continue; +/* -------------------------------------------------------------------- */ +/* Transform to ellipsoidal heights if needed */ +/* -------------------------------------------------------------------- */ +static int geometric_to_orthometric (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) { + int err; + if (0==P->has_geoid_vgrids) + return 0; + if (z==0) + return PJD_ERR_GEOCENTRIC; + err = pj_apply_vgridshift (P, "sgeoidgrids", + &(P->vgridlist_geoid), + &(P->vgridlist_geoid_count), + dir==PJ_FWD ? 1 : 0, n, dist, x, y, z ); + if (err) + return pj_ctx_get_errno(P->ctx); + return 0; +} - projected_loc = pj_fwd3d( geodetic_loc, dstdefn); - if( dstdefn->ctx->last_errno != 0 ) - { - if( (dstdefn->ctx->last_errno != 33 /*EDOM*/ - && dstdefn->ctx->last_errno != 34 /*ERANGE*/ ) - && (dstdefn->ctx->last_errno > 0 - || dstdefn->ctx->last_errno < -44 || point_count == 1 - || transient_error[-dstdefn->ctx->last_errno] == 0 ) ) - return dstdefn->ctx->last_errno; - else - { - projected_loc.u = HUGE_VAL; - projected_loc.v = HUGE_VAL; - projected_loc.w = HUGE_VAL; - } - } - x[point_offset*i] = projected_loc.u; - y[point_offset*i] = projected_loc.v; - z[point_offset*i] = projected_loc.w; - } - } - else - { - for( i = 0; i < point_count; i++ ) - { - XY projected_loc; - LP geodetic_loc; +/* -------------------------------------------------------------------- */ +/* Convert datums if needed, and possible. */ +/* -------------------------------------------------------------------- */ +static int datum_transform (PJ *P, PJ *Q, long n, int dist, double *x, double *y, double *z) { + if (0==pj_datum_transform (P, Q, n, dist, x, y, z)) + return 0; + if (P->ctx->last_errno) + return P->ctx->last_errno; + return Q->ctx->last_errno; +} - geodetic_loc.u = x[point_offset*i]; - geodetic_loc.v = y[point_offset*i]; - if( geodetic_loc.u == HUGE_VAL ) - continue; - projected_loc = pj_fwd( geodetic_loc, dstdefn ); - if( dstdefn->ctx->last_errno != 0 ) - { - if( (dstdefn->ctx->last_errno != 33 /*EDOM*/ - && dstdefn->ctx->last_errno != 34 /*ERANGE*/ ) - && (dstdefn->ctx->last_errno > 0 - || dstdefn->ctx->last_errno < -44 || point_count == 1 - || transient_error[-dstdefn->ctx->last_errno] == 0 ) ) - return dstdefn->ctx->last_errno; - else - { - projected_loc.u = HUGE_VAL; - projected_loc.v = HUGE_VAL; - } - } - x[point_offset*i] = projected_loc.u; - y[point_offset*i] = projected_loc.v; - } - } - } /* -------------------------------------------------------------------- */ /* If a wrapping center other than 0 is provided, rewrap around */ /* the suggested center (for latlong coordinate systems only). */ /* -------------------------------------------------------------------- */ - else if( dstdefn->is_latlong && dstdefn->is_long_wrap_set ) - { - for( i = 0; i < point_count; i++ ) - { - double val = x[point_offset*i]; - if( val == HUGE_VAL ) - continue; +static int long_wrap (PJ *P, long n, int dist, double *x) { + long i; - /* Get fast in ] -2 PI, 2 PI [ range */ - val = fmod(val, M_TWOPI); - while( val < dstdefn->long_wrap_center - M_PI ) - val += M_TWOPI; - while( val > dstdefn->long_wrap_center + M_PI ) - val -= M_TWOPI; - x[point_offset*i] = val; - } - } + /* Nothing to do? */ + if (P->is_geocent) + return 0; + if (!P->is_long_wrap_set) + return 0; + if (!P->is_latlong) + return 0; -/* -------------------------------------------------------------------- */ -/* Transform normalized axes into unusual output coordinate axis */ -/* orientation if needed. */ -/* -------------------------------------------------------------------- */ - if( strcmp(dstdefn->axis,"enu") != 0 ) - { - int err; + for (i = 0; i < n; i++ ) { + double val = x[dist*i]; + if (val == HUGE_VAL) + continue; - err = pj_adjust_axis( dstdefn->ctx, dstdefn->axis, - 1, point_count, point_offset, x, y, z ); - if( err != 0 ) - return err; + /* Get fast in ] -2 PI, 2 PI [ range */ + val = fmod(val, M_TWOPI); + while( val < P->long_wrap_center - M_PI ) + val += M_TWOPI; + while( val > P->long_wrap_center + M_PI ) + val -= M_TWOPI; + x[dist*i] = val; } + return 0; +} + + + +/************************************************************************/ +/* pj_transform() */ +/* */ +/* Currently this function doesn't recognise if two projections */ +/* are identical (to short circuit reprojection) because it is */ +/* difficult to compare PJ structures (since there are some */ +/* projection specific components). */ +/************************************************************************/ + +int pj_transform( + PJ *srcdefn, PJ *dstdefn, + long point_count, int point_offset, + double *x, double *y, double *z +){ + int err; + + srcdefn->ctx->last_errno = 0; + dstdefn->ctx->last_errno = 0; + + if( point_offset == 0 ) + point_offset = 1; + + /* Bring input to "normal form": longitude, latitude, ellipsoidal height */ + + err = adjust_axes (srcdefn, PJ_INV, point_count, point_offset, x, y, z); + if (err) + return err; + err = geographic_to_cartesian (srcdefn, PJ_INV, point_count, point_offset, x, y, z); + if (err) + return err; + err = projected_to_geographic (srcdefn, point_count, point_offset, x, y, z); + if (err) + return err; + err = prime_meridian (srcdefn, PJ_INV, point_count, point_offset, x); + if (err) + return err; + err = height_unit (srcdefn, PJ_INV, point_count, point_offset, z); + if (err) + return err; + err = geometric_to_orthometric (srcdefn, PJ_INV, point_count, point_offset, x, y, z); + if (err) + return err; + + /* At the center of the process we do the datum shift (if needed) */ + + err = datum_transform(srcdefn, dstdefn, point_count, point_offset, x, y, z ); + if (err) + return err; + + /* Now get out on the other side: Bring "normal form" to output form */ + + err = geometric_to_orthometric (dstdefn, PJ_FWD, point_count, point_offset, x, y, z); + if (err) + return err; + err = height_unit (dstdefn, PJ_FWD, point_count, point_offset, z); + if (err) + return err; + err = prime_meridian (dstdefn, PJ_FWD, point_count, point_offset, x); + if (err) + return err; + err = geographic_to_cartesian (dstdefn, PJ_FWD, point_count, point_offset, x, y, z); + if (err) + return err; + err = geographic_to_projected (dstdefn, point_count, point_offset, x, y, z); + if (err) + return err; + err = long_wrap (dstdefn, point_count, point_offset, x); + if (err) + return err; + err = adjust_axes (dstdefn, PJ_FWD, point_count, point_offset, x, y, z); + if (err) + return err; return 0; } + + /************************************************************************/ /* pj_geodetic_to_geocentric() */ /************************************************************************/ @@ -324,7 +324,7 @@ size_t proj_trans_generic ( PJ_COORD proj_coord (double x, double y, double z, double t); /* Measure internal consistency - in forward or inverse direction */ -double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo); +double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coord); /* Geodesic distance between two points with angular 2D coordinates */ double proj_lp_dist (const PJ *P, PJ_COORD a, PJ_COORD b); @@ -374,7 +374,7 @@ double proj_torad (double angle_in_degrees); double proj_todeg (double angle_in_radians); /* Geographical to geocentric latitude - another of the "simple, but useful" */ -PJ_COORD proj_geocentric_latitude (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo); +PJ_COORD proj_geocentric_latitude (const PJ *P, PJ_DIRECTION direction, PJ_COORD coord); double proj_dmstor(const char *is, char **rs); char* proj_rtodms(char *s, double r, int pos, int neg); diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index fb20978b..f0ea7083 100644 --- a/src/proj_4D_api.c +++ b/src/proj_4D_api.c @@ -113,7 +113,7 @@ double proj_xyz_dist (PJ_COORD a, PJ_COORD b) { /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */ -double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) { +double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coord) { int i; PJ_COORD t, org; @@ -126,9 +126,9 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) { } /* in the first half-step, we generate the output value */ - org = *coo; - *coo = proj_trans (P, direction, org); - t = *coo; + org = *coord; + *coord = proj_trans (P, direction, org); + t = *coord; /* now we take n-1 full steps in inverse direction: We are */ /* out of phase due to the half step already taken */ @@ -148,26 +148,26 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) { /**************************************************************************************/ -PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) { +PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coord) { /*************************************************************************************** -Apply the transformation P to the coordinate coo, preferring the 4D interfaces if +Apply the transformation P to the coordinate coord, preferring the 4D interfaces if available. See also pj_approx_2D_trans and pj_approx_3D_trans in pj_internal.c, which work similarly, but prefers the 2D resp. 3D interfaces if available. ***************************************************************************************/ if (0==P) - return coo; + return coord; if (P->inverted) direction = -direction; switch (direction) { case PJ_FWD: - return pj_fwd4d (coo, P); + return pj_fwd4d (coord, P); case PJ_INV: - return pj_inv4d (coo, P); + return pj_inv4d (coord, P); case PJ_IDENT: - return coo; + return coord; default: break; } @@ -357,7 +357,7 @@ size_t proj_trans_generic ( /*************************************************************************************/ -PJ_COORD proj_geocentric_latitude (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo) { +PJ_COORD proj_geocentric_latitude (const PJ *P, PJ_DIRECTION direction, PJ_COORD coord) { /************************************************************************************** Convert geographical latitude to geocentric (or the other way round if direction = PJ_INV) @@ -373,13 +373,13 @@ PJ_COORD proj_geocentric_latitude (const PJ *P, PJ_DIRECTION direction, PJ_COORD consequently, the input is copied directly to the output. **************************************************************************************/ const double limit = M_HALFPI - 1e-9; - PJ_COORD res = coo; - if ((coo.lp.phi > limit) || (coo.lp.phi < -limit) || (P->es==0)) + PJ_COORD res = coord; + if ((coord.lp.phi > limit) || (coord.lp.phi < -limit) || (P->es==0)) return res; if (direction==PJ_FWD) - res.lp.phi = atan (P->one_es * tan (coo.lp.phi) ); + res.lp.phi = atan (P->one_es * tan (coord.lp.phi) ); else - res.lp.phi = atan (P->rone_es * tan (coo.lp.phi) ); + res.lp.phi = atan (P->rone_es * tan (coord.lp.phi) ); return res; } diff --git a/src/proj_internal.h b/src/proj_internal.h index 071510c7..e72ce578 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -33,7 +33,7 @@ #endif #include <math.h> /* For M_PI */ -#include <proj.h> +#include "proj.h" #ifdef PROJECTS_H #error proj_internal.h must be included before projects.h @@ -73,8 +73,6 @@ enum pj_io_units { enum pj_io_units pj_left (PJ *P); enum pj_io_units pj_right (PJ *P); -PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD obs); - PJ_COORD proj_coord_error (void); void proj_context_errno_set (PJ_CONTEXT *ctx, int err); |