aboutsummaryrefslogtreecommitdiff
path: root/src/transform.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/transform.cpp')
-rw-r--r--src/transform.cpp1047
1 files changed, 1047 insertions, 0 deletions
diff --git a/src/transform.cpp b/src/transform.cpp
new file mode 100644
index 00000000..433fc017
--- /dev/null
+++ b/src/transform.cpp
@@ -0,0 +1,1047 @@
+/******************************************************************************
+ * Project: PROJ.4
+ * Purpose: Perform overall coordinate system to coordinate system
+ * transformations (pj_transform() function) including reprojection
+ * and datum shifting.
+ * Author: Frank Warmerdam, warmerdam@pobox.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#include <errno.h>
+#include <math.h>
+#include <string.h>
+
+#include "proj.h"
+#include "projects.h"
+#include "geocent.h"
+
+static int adjust_axis( projCtx ctx, 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
+
+#ifndef SRS_WGS84_ESQUARED
+#define SRS_WGS84_ESQUARED 0.0066943799901413165
+#endif
+
+#define Dx_BF (defn->datum_params[0])
+#define Dy_BF (defn->datum_params[1])
+#define Dz_BF (defn->datum_params[2])
+#define Rx_BF (defn->datum_params[3])
+#define Ry_BF (defn->datum_params[4])
+#define Rz_BF (defn->datum_params[5])
+#define M_BF (defn->datum_params[6])
+
+/*
+** 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.
+**
+** At some point we might want to move this array in with the error message
+** list or something, but while experimenting with it this should be fine.
+**
+**
+** NOTE (2017-10-01): Non-transient errors really should have resulted in a
+** PJ==0 during initialization, and hence should be handled at the level
+** before calling pj_transform. The only obvious example of the contrary
+** appears to be the PJD_ERR_GRID_AREA case, which may also be taken to
+** mean "no grids available"
+**
+**
+*/
+
+static const int transient_error[60] = {
+ /* 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, 0, 0,
+ /* 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 };
+
+
+/* -------------------------------------------------------------------- */
+/* 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 */
+/* standard form if needed. */
+/* -------------------------------------------------------------------- */
+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 adjust_axis( P->ctx, P->axis,
+ dir==PJ_FWD ? 1: 0, n, dist, x, y, z );
+}
+
+
+
+/* ----------------------------------------------------------------------- */
+/* Transform geographic (lat/long) source coordinates to */
+/* cartesian ("geocentric"), 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 == nullptr ) {
+ 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 destination points to projection coordinates, if */
+/* desired. */
+/* */
+/* Ought to fold this into projected_to_geographic */
+/* -------------------------------------------------------------------- */
+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 && !P->geoc && P->vto_meter == 1.0)
+ return 0;
+ if (P->is_geocent)
+ return 0;
+
+ if(P->fwd3d != nullptr && !(z == nullptr && P->is_latlong))
+ {
+ /* Three dimensions must be defined */
+ if ( z == nullptr)
+ {
+ pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
+ return PJD_ERR_GEOCENTRIC;
+ }
+
+ for( i = 0; i < n; i++ )
+ {
+ XYZ projected_loc;
+ LPZ geodetic_loc;
+
+ geodetic_loc.lam = x[dist*i];
+ geodetic_loc.phi = y[dist*i];
+ geodetic_loc.z = z[dist*i];
+
+ if (geodetic_loc.lam == HUGE_VAL)
+ continue;
+
+ proj_errno_reset( P );
+ projected_loc = pj_fwd3d( geodetic_loc, P);
+ if( P->ctx->last_errno != 0 )
+ {
+ if( (P->ctx->last_errno != EDOM
+ && P->ctx->last_errno != ERANGE)
+ && (P->ctx->last_errno > 0
+ || P->ctx->last_errno < -44 || n == 1
+ || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
+ {
+ return P->ctx->last_errno;
+ }
+ else
+ {
+ projected_loc.x = HUGE_VAL;
+ projected_loc.y = HUGE_VAL;
+ projected_loc.z = HUGE_VAL;
+ }
+ }
+
+ x[dist*i] = projected_loc.x;
+ y[dist*i] = projected_loc.y;
+ z[dist*i] = projected_loc.z;
+ }
+ return 0;
+ }
+
+ for( i = 0; i <n; i++ )
+ {
+ XY projected_loc;
+ LP geodetic_loc;
+
+ geodetic_loc.lam = x[dist*i];
+ geodetic_loc.phi = y[dist*i];
+
+ if( geodetic_loc.lam == HUGE_VAL )
+ continue;
+
+ proj_errno_reset( P );
+ projected_loc = pj_fwd( geodetic_loc, P );
+ if( P->ctx->last_errno != 0 )
+ {
+ if( (P->ctx->last_errno != EDOM
+ && P->ctx->last_errno != ERANGE)
+ && (P->ctx->last_errno > 0
+ || P->ctx->last_errno < -44 || n == 1
+ || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
+ {
+ return P->ctx->last_errno;
+ }
+ else
+ {
+ projected_loc.x = HUGE_VAL;
+ projected_loc.y = HUGE_VAL;
+ }
+ }
+
+ x[dist*i] = projected_loc.x;
+ y[dist*i] = projected_loc.y;
+ }
+ return 0;
+}
+
+
+
+
+
+/* ----------------------------------------------------------------------- */
+/* 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 && !P->geoc && P->vto_meter == 1.0)
+ return 0;
+ if (P->is_geocent)
+ return 0;
+
+ /* Check first if projection is invertible. */
+ if( (P->inv3d == nullptr) && (P->inv == nullptr))
+ {
+ pj_ctx_set_errno(pj_get_ctx(P), PJD_ERR_NON_CONV_INV_MERI_DIST);
+ pj_log( pj_get_ctx(P), PJ_LOG_ERROR,
+ "pj_transform(): source projection not invertable" );
+ return PJD_ERR_NON_CONV_INV_MERI_DIST;
+ }
+
+ /* If invertible - First try inv3d if defined */
+ if (P->inv3d != nullptr && !(z == nullptr && P->is_latlong))
+ {
+ /* Three dimensions must be defined */
+ if ( z == nullptr)
+ {
+ pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
+ return PJD_ERR_GEOCENTRIC;
+ }
+
+ for (i=0; i < n; i++)
+ {
+ XYZ projected_loc;
+ LPZ geodetic_loc;
+
+ projected_loc.x = x[dist*i];
+ projected_loc.y = y[dist*i];
+ projected_loc.z = z[dist*i];
+
+ if (projected_loc.x == HUGE_VAL)
+ continue;
+
+ proj_errno_reset( P );
+ geodetic_loc = pj_inv3d(projected_loc, P);
+ if( P->ctx->last_errno != 0 )
+ {
+ if( (P->ctx->last_errno != EDOM
+ && P->ctx->last_errno != ERANGE)
+ && (P->ctx->last_errno > 0
+ || P->ctx->last_errno < -44 || n == 1
+ || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
+ {
+ return P->ctx->last_errno;
+ }
+ else
+ {
+ geodetic_loc.lam = HUGE_VAL;
+ geodetic_loc.phi = HUGE_VAL;
+ geodetic_loc.z = HUGE_VAL;
+ }
+ }
+
+ x[dist*i] = geodetic_loc.lam;
+ y[dist*i] = geodetic_loc.phi;
+ z[dist*i] = geodetic_loc.z;
+
+ }
+ return 0;
+ }
+
+ /* Fallback to the original PROJ.4 API 2d inversion - inv */
+ for( i = 0; i < n; i++ ) {
+ XY projected_loc;
+ LP geodetic_loc;
+
+ projected_loc.x = x[dist*i];
+ projected_loc.y = y[dist*i];
+
+ if( projected_loc.x == HUGE_VAL )
+ continue;
+
+ proj_errno_reset( P );
+ geodetic_loc = pj_inv( projected_loc, P );
+ if( P->ctx->last_errno != 0 )
+ {
+ if( (P->ctx->last_errno != EDOM
+ && P->ctx->last_errno != ERANGE)
+ && (P->ctx->last_errno > 0
+ || P->ctx->last_errno < -44 || n == 1
+ || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
+ {
+ return P->ctx->last_errno;
+ }
+ else
+ {
+ geodetic_loc.lam = HUGE_VAL;
+ geodetic_loc.phi = HUGE_VAL;
+ }
+ }
+
+ x[dist*i] = geodetic_loc.lam;
+ y[dist*i] = geodetic_loc.phi;
+ }
+ return 0;
+}
+
+
+
+/* -------------------------------------------------------------------- */
+/* Adjust for the prime meridian if needed. */
+/* -------------------------------------------------------------------- */
+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;
+}
+
+
+
+/* -------------------------------------------------------------------- */
+/* Adjust for vertical scale factor if needed */
+/* -------------------------------------------------------------------- */
+static int height_unit (PJ *P, PJ_DIRECTION dir, long n, int dist, double *z) {
+ int i;
+ double fac = P->vto_meter;
+
+ if (PJ_FWD==dir)
+ fac = P->vfr_meter;
+
+ /* Nothing to do? */
+ if (fac==1.0)
+ return 0;
+ if (nullptr==z)
+ return 0;
+ if (P->is_latlong)
+ return 0; /* done in pj_inv3d() / pj_fwd3d() */
+
+ for (i = 0; i < n; i++)
+ if (z[dist*i] != HUGE_VAL )
+ z[dist*i] *= fac;
+
+ return 0;
+}
+
+
+
+/* -------------------------------------------------------------------- */
+/* 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==nullptr)
+ 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;
+}
+
+
+
+/* -------------------------------------------------------------------- */
+/* 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;
+}
+
+
+
+
+
+/* -------------------------------------------------------------------- */
+/* If a wrapping center other than 0 is provided, rewrap around */
+/* the suggested center (for latlong coordinate systems only). */
+/* -------------------------------------------------------------------- */
+static int long_wrap (PJ *P, long n, int dist, double *x) {
+ long i;
+
+ /* Nothing to do? */
+ if (P->is_geocent)
+ return 0;
+ if (!P->is_long_wrap_set)
+ return 0;
+ if (!P->is_latlong)
+ return 0;
+
+ for (i = 0; i < n; i++ ) {
+ double val = x[dist*i];
+ if (val == HUGE_VAL)
+ continue;
+
+ /* 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 *src, PJ *dst,
+ long point_count, int point_offset,
+ double *x, double *y, double *z
+){
+ int err;
+
+ src->ctx->last_errno = 0;
+ dst->ctx->last_errno = 0;
+
+ if( point_offset == 0 )
+ point_offset = 1;
+
+ /* Bring input to "normal form": longitude, latitude, ellipsoidal height */
+
+ err = adjust_axes (src, PJ_INV, point_count, point_offset, x, y, z);
+ if (err)
+ return err;
+ err = geographic_to_cartesian (src, PJ_INV, point_count, point_offset, x, y, z);
+ if (err)
+ return err;
+ err = projected_to_geographic (src, point_count, point_offset, x, y, z);
+ if (err)
+ return err;
+ err = prime_meridian (src, PJ_INV, point_count, point_offset, x);
+ if (err)
+ return err;
+ err = height_unit (src, PJ_INV, point_count, point_offset, z);
+ if (err)
+ return err;
+ err = geometric_to_orthometric (src, 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(src, dst, 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 (dst, PJ_FWD, point_count, point_offset, x, y, z);
+ if (err)
+ return err;
+ err = height_unit (dst, PJ_FWD, point_count, point_offset, z);
+ if (err)
+ return err;
+ err = prime_meridian (dst, PJ_FWD, point_count, point_offset, x);
+ if (err)
+ return err;
+ err = geographic_to_cartesian (dst, PJ_FWD, point_count, point_offset, x, y, z);
+ if (err)
+ return err;
+ err = geographic_to_projected (dst, point_count, point_offset, x, y, z);
+ if (err)
+ return err;
+ err = long_wrap (dst, point_count, point_offset, x);
+ if (err)
+ return err;
+ err = adjust_axes (dst, PJ_FWD, point_count, point_offset, x, y, z);
+ if (err)
+ return err;
+
+ return 0;
+}
+
+
+
+/************************************************************************/
+/* pj_geodetic_to_geocentric() */
+/************************************************************************/
+
+int pj_geodetic_to_geocentric( double a, double es,
+ long point_count, int point_offset,
+ double *x, double *y, double *z )
+
+{
+ double b;
+ int i;
+ GeocentricInfo gi;
+ int ret_errno = 0;
+
+ if( es == 0.0 )
+ b = a;
+ else
+ b = a * sqrt(1-es);
+
+ if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
+ {
+ return PJD_ERR_GEOCENTRIC;
+ }
+
+ for( i = 0; i < point_count; i++ )
+ {
+ long io = i * point_offset;
+
+ if( x[io] == HUGE_VAL )
+ continue;
+
+ if( pj_Convert_Geodetic_To_Geocentric( &gi, y[io], x[io], z[io],
+ x+io, y+io, z+io ) != 0 )
+ {
+ ret_errno = PJD_ERR_LAT_OR_LON_EXCEED_LIMIT;
+ x[io] = y[io] = HUGE_VAL;
+ /* but keep processing points! */
+ }
+ }
+
+ return ret_errno;
+}
+
+/************************************************************************/
+/* pj_geocentric_to_geodetic() */
+/************************************************************************/
+
+int pj_geocentric_to_geodetic( double a, double es,
+ long point_count, int point_offset,
+ double *x, double *y, double *z )
+
+{
+ double b;
+ int i;
+ GeocentricInfo gi;
+
+ if( es == 0.0 )
+ b = a;
+ else
+ b = a * sqrt(1-es);
+
+ if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
+ {
+ return PJD_ERR_GEOCENTRIC;
+ }
+
+ for( i = 0; i < point_count; i++ )
+ {
+ long io = i * point_offset;
+
+ if( x[io] == HUGE_VAL )
+ continue;
+
+ pj_Convert_Geocentric_To_Geodetic( &gi, x[io], y[io], z[io],
+ y+io, x+io, z+io );
+ }
+
+ return 0;
+}
+
+/************************************************************************/
+/* pj_compare_datums() */
+/* */
+/* Returns TRUE if the two datums are identical, otherwise */
+/* FALSE. */
+/************************************************************************/
+
+int pj_compare_datums( PJ *srcdefn, PJ *dstdefn )
+
+{
+ if( srcdefn->datum_type != dstdefn->datum_type )
+ {
+ return 0;
+ }
+ else if( srcdefn->a_orig != dstdefn->a_orig
+ || ABS(srcdefn->es_orig - dstdefn->es_orig) > 0.000000000050 )
+ {
+ /* the tolerance for es is to ensure that GRS80 and WGS84 are
+ considered identical */
+ return 0;
+ }
+ else if( srcdefn->datum_type == PJD_3PARAM )
+ {
+ return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
+ && srcdefn->datum_params[1] == dstdefn->datum_params[1]
+ && srcdefn->datum_params[2] == dstdefn->datum_params[2]);
+ }
+ else if( srcdefn->datum_type == PJD_7PARAM )
+ {
+ return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
+ && srcdefn->datum_params[1] == dstdefn->datum_params[1]
+ && srcdefn->datum_params[2] == dstdefn->datum_params[2]
+ && srcdefn->datum_params[3] == dstdefn->datum_params[3]
+ && srcdefn->datum_params[4] == dstdefn->datum_params[4]
+ && srcdefn->datum_params[5] == dstdefn->datum_params[5]
+ && srcdefn->datum_params[6] == dstdefn->datum_params[6]);
+ }
+ else if( srcdefn->datum_type == PJD_GRIDSHIFT )
+ {
+ const char* srcnadgrids =
+ pj_param(srcdefn->ctx, srcdefn->params,"snadgrids").s;
+ const char* dstnadgrids =
+ pj_param(dstdefn->ctx, dstdefn->params,"snadgrids").s;
+ return srcnadgrids != nullptr && dstnadgrids != nullptr &&
+ strcmp( srcnadgrids, dstnadgrids ) == 0;
+ }
+ else
+ return 1;
+}
+
+/************************************************************************/
+/* pj_geocentic_to_wgs84() */
+/************************************************************************/
+
+static
+int pj_geocentric_to_wgs84( PJ *defn,
+ long point_count, int point_offset,
+ double *x, double *y, double *z )
+
+{
+ int i;
+
+ if( defn->datum_type == PJD_3PARAM )
+ {
+ for( i = 0; i < point_count; i++ )
+ {
+ long io = i * point_offset;
+
+ if( x[io] == HUGE_VAL )
+ continue;
+
+ x[io] = x[io] + Dx_BF;
+ y[io] = y[io] + Dy_BF;
+ z[io] = z[io] + Dz_BF;
+ }
+ }
+ else if( defn->datum_type == PJD_7PARAM )
+ {
+ for( i = 0; i < point_count; i++ )
+ {
+ long io = i * point_offset;
+ double x_out, y_out, z_out;
+
+ if( x[io] == HUGE_VAL )
+ continue;
+
+ x_out = M_BF*( x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
+ y_out = M_BF*( Rz_BF*x[io] + y[io] - Rx_BF*z[io]) + Dy_BF;
+ z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] + z[io]) + Dz_BF;
+
+ x[io] = x_out;
+ y[io] = y_out;
+ z[io] = z_out;
+ }
+ }
+
+ return 0;
+}
+
+/************************************************************************/
+/* pj_geocentic_from_wgs84() */
+/************************************************************************/
+
+static
+int pj_geocentric_from_wgs84( PJ *defn,
+ long point_count, int point_offset,
+ double *x, double *y, double *z )
+
+{
+ int i;
+
+ if( defn->datum_type == PJD_3PARAM )
+ {
+ for( i = 0; i < point_count; i++ )
+ {
+ long io = i * point_offset;
+
+ if( x[io] == HUGE_VAL )
+ continue;
+
+ x[io] = x[io] - Dx_BF;
+ y[io] = y[io] - Dy_BF;
+ z[io] = z[io] - Dz_BF;
+ }
+ }
+ else if( defn->datum_type == PJD_7PARAM )
+ {
+ for( i = 0; i < point_count; i++ )
+ {
+ long io = i * point_offset;
+ double x_tmp, y_tmp, z_tmp;
+
+ if( x[io] == HUGE_VAL )
+ continue;
+
+ x_tmp = (x[io] - Dx_BF) / M_BF;
+ y_tmp = (y[io] - Dy_BF) / M_BF;
+ z_tmp = (z[io] - Dz_BF) / M_BF;
+
+ x[io] = x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp;
+ y[io] = -Rz_BF*x_tmp + y_tmp + Rx_BF*z_tmp;
+ z[io] = Ry_BF*x_tmp - Rx_BF*y_tmp + z_tmp;
+ }
+ }
+
+ return 0;
+}
+
+/************************************************************************/
+/* pj_datum_transform() */
+/* */
+/* The input should be long/lat/z coordinates in radians in the */
+/* source datum, and the output should be long/lat/z */
+/* coordinates in radians in the destination datum. */
+/************************************************************************/
+
+int pj_datum_transform( PJ *src, PJ *dst,
+ long point_count, int point_offset,
+ double *x, double *y, double *z )
+
+{
+ double src_a, src_es, dst_a, dst_es;
+ int z_is_temp = FALSE;
+
+/* -------------------------------------------------------------------- */
+/* We cannot do any meaningful datum transformation if either */
+/* the source or destination are of an unknown datum type */
+/* (ie. only a +ellps declaration, no +datum). This is new */
+/* behavior for PROJ 4.6.0. */
+/* -------------------------------------------------------------------- */
+ if( src->datum_type == PJD_UNKNOWN
+ || dst->datum_type == PJD_UNKNOWN )
+ return 0;
+
+/* -------------------------------------------------------------------- */
+/* Short cut if the datums are identical. */
+/* -------------------------------------------------------------------- */
+ if( pj_compare_datums( src, dst ) )
+ return 0;
+
+ src_a = src->a_orig;
+ src_es = src->es_orig;
+
+ dst_a = dst->a_orig;
+ dst_es = dst->es_orig;
+
+/* -------------------------------------------------------------------- */
+/* Create a temporary Z array if one is not provided. */
+/* -------------------------------------------------------------------- */
+ if( z == nullptr )
+ {
+ size_t bytes = sizeof(double) * point_count * point_offset;
+ z = (double *) pj_malloc(bytes);
+ memset( z, 0, bytes );
+ z_is_temp = TRUE;
+ }
+
+#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 */
+/* coordinates. */
+/* -------------------------------------------------------------------- */
+ if( src->datum_type == PJD_GRIDSHIFT )
+ {
+ pj_apply_gridshift_2( src, 0, point_count, point_offset, x, y, z );
+ CHECK_RETURN(src);
+
+ src_a = SRS_WGS84_SEMIMAJOR;
+ src_es = SRS_WGS84_ESQUARED;
+ }
+
+ if( dst->datum_type == PJD_GRIDSHIFT )
+ {
+ dst_a = SRS_WGS84_SEMIMAJOR;
+ dst_es = SRS_WGS84_ESQUARED;
+ }
+
+/* ==================================================================== */
+/* Do we need to go through geocentric coordinates? */
+/* ==================================================================== */
+ if( src_es != dst_es || src_a != dst_a
+ || src->datum_type == PJD_3PARAM
+ || src->datum_type == PJD_7PARAM
+ || dst->datum_type == PJD_3PARAM
+ || dst->datum_type == PJD_7PARAM)
+ {
+/* -------------------------------------------------------------------- */
+/* Convert to geocentric coordinates. */
+/* -------------------------------------------------------------------- */
+ src->ctx->last_errno =
+ pj_geodetic_to_geocentric( src_a, src_es,
+ point_count, point_offset, x, y, z );
+ CHECK_RETURN(src);
+
+/* -------------------------------------------------------------------- */
+/* Convert between datums. */
+/* -------------------------------------------------------------------- */
+ if( src->datum_type == PJD_3PARAM
+ || src->datum_type == PJD_7PARAM )
+ {
+ pj_geocentric_to_wgs84( src, point_count, point_offset,x,y,z);
+ CHECK_RETURN(src);
+ }
+
+ if( dst->datum_type == PJD_3PARAM
+ || dst->datum_type == PJD_7PARAM )
+ {
+ pj_geocentric_from_wgs84( dst, point_count,point_offset,x,y,z);
+ CHECK_RETURN(dst);
+ }
+
+/* -------------------------------------------------------------------- */
+/* Convert back to geodetic coordinates. */
+/* -------------------------------------------------------------------- */
+ dst->ctx->last_errno =
+ pj_geocentric_to_geodetic( dst_a, dst_es,
+ point_count, point_offset, x, y, z );
+ CHECK_RETURN(dst);
+ }
+
+/* -------------------------------------------------------------------- */
+/* Apply grid shift to destination if required. */
+/* -------------------------------------------------------------------- */
+ if( dst->datum_type == PJD_GRIDSHIFT )
+ {
+ pj_apply_gridshift_2( dst, 1, point_count, point_offset, x, y, z );
+ CHECK_RETURN(dst);
+ }
+
+ if( z_is_temp )
+ pj_dalloc( z );
+
+ return 0;
+}
+
+/************************************************************************/
+/* adjust_axis() */
+/* */
+/* Normalize or de-normalized the x/y/z axes. The normal form */
+/* is "enu" (easting, northing, up). */
+/************************************************************************/
+static int adjust_axis( projCtx ctx,
+ 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_ctx_set_errno( ctx, 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 == nullptr )
+ 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_ctx_set_errno( ctx, PJD_ERR_AXIS );
+ return PJD_ERR_AXIS;
+ }
+ } /* i_axis */
+ } /* i (point) */
+ }
+
+ return 0;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 18:28:34 +0000