"Fossies" - the Fresh Open Source Software Archive

Member "proj-6.2.1/src/transform.cpp" (6 Mar 2019, 34433 Bytes) of package /linux/privat/proj-6.2.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "transform.cpp" see the Fossies "Dox" file reference documentation.

    1 /******************************************************************************
    2  * Project:  PROJ.4
    3  * Purpose:  Perform overall coordinate system to coordinate system
    4  *           transformations (pj_transform() function) including reprojection
    5  *           and datum shifting.
    6  * Author:   Frank Warmerdam, warmerdam@pobox.com
    7  *
    8  ******************************************************************************
    9  * Copyright (c) 2000, Frank Warmerdam
   10  *
   11  * Permission is hereby granted, free of charge, to any person obtaining a
   12  * copy of this software and associated documentation files (the "Software"),
   13  * to deal in the Software without restriction, including without limitation
   14  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   15  * and/or sell copies of the Software, and to permit persons to whom the
   16  * Software is furnished to do so, subject to the following conditions:
   17  *
   18  * The above copyright notice and this permission notice shall be included
   19  * in all copies or substantial portions of the Software.
   20  *
   21  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
   22  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   23  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
   24  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   25  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   26  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
   27  * DEALINGS IN THE SOFTWARE.
   28  *****************************************************************************/
   29 
   30 #include <errno.h>
   31 #include <math.h>
   32 #include <string.h>
   33 
   34 #include "proj.h"
   35 #include "proj_internal.h"
   36 #include "geocent.h"
   37 
   38 static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
   39                            long point_count, int point_offset,
   40                            double *x, double *y, double *z );
   41 
   42 #ifndef SRS_WGS84_SEMIMAJOR
   43 #define SRS_WGS84_SEMIMAJOR 6378137.0
   44 #endif
   45 
   46 #ifndef SRS_WGS84_ESQUARED
   47 #define SRS_WGS84_ESQUARED 0.0066943799901413165
   48 #endif
   49 
   50 #define Dx_BF (defn->datum_params[0])
   51 #define Dy_BF (defn->datum_params[1])
   52 #define Dz_BF (defn->datum_params[2])
   53 #define Rx_BF (defn->datum_params[3])
   54 #define Ry_BF (defn->datum_params[4])
   55 #define Rz_BF (defn->datum_params[5])
   56 #define M_BF  (defn->datum_params[6])
   57 
   58 /*
   59 ** This table is intended to indicate for any given error code
   60 ** whether that error will occur for all locations (ie.
   61 ** it is a problem with the coordinate system as a whole) in which case the
   62 ** value would be 0, or if the problem is with the point being transformed
   63 ** in which case the value is 1.
   64 **
   65 ** At some point we might want to move this array in with the error message
   66 ** list or something, but while experimenting with it this should be fine.
   67 **
   68 **
   69 ** NOTE (2017-10-01): Non-transient errors really should have resulted in a
   70 ** PJ==0 during initialization, and hence should be handled at the level
   71 ** before calling pj_transform. The only obvious example of the contrary
   72 ** appears to be the PJD_ERR_GRID_AREA case, which may also be taken to
   73 ** mean "no grids available"
   74 **
   75 **
   76 */
   77 
   78 static const int transient_error[70] = {
   79     /*             0  1  2  3  4  5  6  7  8  9   */
   80     /* 0 to 9 */   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   81     /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,
   82     /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
   83     /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   84     /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
   85     /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0,
   86     /* 60 to 69 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,};
   87 
   88 
   89 /* -------------------------------------------------------------------- */
   90 /*      Read transient_error[] in a safe way.                           */
   91 /* -------------------------------------------------------------------- */
   92 static int get_transient_error_value(int pos_index)
   93 {
   94     const int array_size =
   95         (int)(sizeof(transient_error) / sizeof(transient_error[0]));
   96     if( pos_index < 0 || pos_index >= array_size ) {
   97         return 0;
   98     }
   99     return transient_error[pos_index];
  100 }
  101 
  102 
  103 /* -------------------------------------------------------------------- */
  104 /*      Transform unusual input coordinate axis orientation to          */
  105 /*      standard form if needed.                                        */
  106 /* -------------------------------------------------------------------- */
  107 static int adjust_axes (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) {
  108     /* Nothing to do? */
  109     if (0==strcmp(P->axis,"enu"))
  110         return 0;
  111 
  112     return adjust_axis( P->ctx, P->axis,
  113                 dir==PJ_FWD ? 1: 0, n, dist, x, y, z );
  114 }
  115 
  116 
  117 
  118 /* ----------------------------------------------------------------------- */
  119 /*    Transform geographic (lat/long) source coordinates to                */
  120 /*    cartesian ("geocentric"), if needed                                  */
  121 /* ----------------------------------------------------------------------- */
  122 static int geographic_to_cartesian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) {
  123     int res;
  124     long i;
  125     double fac = P->to_meter;
  126 
  127     /* Nothing to do? */
  128     if (!P->is_geocent)
  129         return 0;
  130 
  131     if ( z == nullptr ) {
  132         pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
  133         return PJD_ERR_GEOCENTRIC;
  134     }
  135 
  136     if (PJ_FWD==dir) {
  137         fac = P->fr_meter;
  138         res = pj_geodetic_to_geocentric( P->a_orig, P->es_orig, n, dist, x, y, z );
  139         if (res)
  140             return res;
  141     }
  142 
  143     if (fac != 1.0) {
  144         for( i = 0; i < n; i++ ) {
  145             if( x[dist*i] != HUGE_VAL ) {
  146                 x[dist*i] *= fac;
  147                 y[dist*i] *= fac;
  148                 z[dist*i] *= fac;
  149             }
  150         }
  151     }
  152 
  153     if (PJ_FWD==dir)
  154         return 0;
  155     return pj_geocentric_to_geodetic(
  156         P->a_orig, P->es_orig,
  157         n, dist,
  158         x, y, z
  159     );
  160 }
  161 
  162 
  163 
  164 
  165 
  166 
  167 
  168 
  169 
  170 
  171 /* -------------------------------------------------------------------- */
  172 /*      Transform destination points to projection coordinates, if      */
  173 /*      desired.                                                        */
  174 /*                                                                      */
  175 /*      Ought to fold this into projected_to_geographic                 */
  176 /* -------------------------------------------------------------------- */
  177 static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *y, double *z) {
  178     long i;
  179 
  180     /* Nothing to do? */
  181     if (P->is_latlong && !P->geoc && P->vto_meter == 1.0)
  182         return 0;
  183     if (P->is_geocent)
  184         return 0;
  185 
  186     if(P->fwd3d != nullptr && !(z == nullptr && P->is_latlong))
  187     {
  188         /* Three dimensions must be defined */
  189         if ( z == nullptr)
  190         {
  191             pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
  192             return PJD_ERR_GEOCENTRIC;
  193         }
  194 
  195         for( i = 0; i < n; i++ )
  196         {
  197             PJ_XYZ projected_loc;
  198             PJ_LPZ geodetic_loc;
  199 
  200             geodetic_loc.lam = x[dist*i];
  201             geodetic_loc.phi = y[dist*i];
  202             geodetic_loc.z = z[dist*i];
  203 
  204             if (geodetic_loc.lam == HUGE_VAL)
  205                 continue;
  206 
  207             proj_errno_reset( P );
  208             projected_loc = pj_fwd3d( geodetic_loc, P);
  209             if( P->ctx->last_errno != 0 )
  210             {
  211                 if( (P->ctx->last_errno != EDOM
  212                         && P->ctx->last_errno != ERANGE)
  213                     && (P->ctx->last_errno > 0
  214                         || P->ctx->last_errno < -44 || n == 1
  215                         || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
  216                 {
  217                     return P->ctx->last_errno;
  218                 }
  219                 else
  220                 {
  221                     projected_loc.x = HUGE_VAL;
  222                     projected_loc.y = HUGE_VAL;
  223                     projected_loc.z = HUGE_VAL;
  224                 }
  225             }
  226 
  227             x[dist*i] = projected_loc.x;
  228             y[dist*i] = projected_loc.y;
  229             z[dist*i] = projected_loc.z;
  230         }
  231         return 0;
  232     }
  233 
  234     for( i = 0; i <n; i++ )
  235     {
  236         PJ_XY         projected_loc;
  237         PJ_LP          geodetic_loc;
  238 
  239         geodetic_loc.lam = x[dist*i];
  240         geodetic_loc.phi = y[dist*i];
  241 
  242         if( geodetic_loc.lam == HUGE_VAL )
  243             continue;
  244 
  245         proj_errno_reset( P );
  246         projected_loc = pj_fwd( geodetic_loc, P );
  247         if( P->ctx->last_errno != 0 )
  248         {
  249             if( (P->ctx->last_errno != EDOM
  250                     && P->ctx->last_errno != ERANGE)
  251                 && (P->ctx->last_errno > 0
  252                     || P->ctx->last_errno < -44 || n == 1
  253                     || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
  254             {
  255                 return P->ctx->last_errno;
  256             }
  257             else
  258             {
  259                 projected_loc.x = HUGE_VAL;
  260                 projected_loc.y = HUGE_VAL;
  261             }
  262         }
  263 
  264         x[dist*i] = projected_loc.x;
  265         y[dist*i] = projected_loc.y;
  266     }
  267     return 0;
  268 }
  269 
  270 
  271 
  272 
  273 
  274 /* ----------------------------------------------------------------------- */
  275 /*    Transform projected source coordinates to lat/long, if needed        */
  276 /* ----------------------------------------------------------------------- */
  277 static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *y, double *z) {
  278     long i;
  279 
  280     /* Nothing to do? */
  281     if (P->is_latlong && !P->geoc && P->vto_meter == 1.0)
  282         return 0;
  283     if (P->is_geocent)
  284         return 0;
  285 
  286     /* Check first if projection is invertible. */
  287     if( (P->inv3d == nullptr) && (P->inv == nullptr))
  288     {
  289         pj_ctx_set_errno(pj_get_ctx(P), PJD_ERR_NON_CONV_INV_MERI_DIST);
  290         pj_log( pj_get_ctx(P), PJ_LOG_ERROR,
  291                 "pj_transform(): source projection not invertable" );
  292         return PJD_ERR_NON_CONV_INV_MERI_DIST;
  293     }
  294 
  295     /* If invertible - First try inv3d if defined */
  296     if (P->inv3d != nullptr && !(z == nullptr && P->is_latlong))
  297     {
  298         /* Three dimensions must be defined */
  299         if ( z == nullptr)
  300         {
  301             pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
  302             return PJD_ERR_GEOCENTRIC;
  303         }
  304 
  305         for (i=0; i < n; i++)
  306         {
  307             PJ_XYZ projected_loc;
  308             PJ_LPZ geodetic_loc;
  309 
  310             projected_loc.x = x[dist*i];
  311             projected_loc.y = y[dist*i];
  312             projected_loc.z = z[dist*i];
  313 
  314             if (projected_loc.x == HUGE_VAL)
  315                 continue;
  316 
  317             proj_errno_reset( P );
  318             geodetic_loc = pj_inv3d(projected_loc, P);
  319             if( P->ctx->last_errno != 0 )
  320             {
  321                 if( (P->ctx->last_errno != EDOM
  322                         && P->ctx->last_errno != ERANGE)
  323                     && (P->ctx->last_errno > 0
  324                         || P->ctx->last_errno < -44 || n == 1
  325                         || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
  326                 {
  327                     return P->ctx->last_errno;
  328                 }
  329                 else
  330                 {
  331                     geodetic_loc.lam = HUGE_VAL;
  332                     geodetic_loc.phi = HUGE_VAL;
  333                     geodetic_loc.z = HUGE_VAL;
  334                 }
  335             }
  336 
  337             x[dist*i] = geodetic_loc.lam;
  338             y[dist*i] = geodetic_loc.phi;
  339             z[dist*i] = geodetic_loc.z;
  340 
  341         }
  342         return 0;
  343     }
  344 
  345     /* Fallback to the original PROJ.4 API 2d inversion - inv */
  346     for( i = 0; i < n; i++ ) {
  347         PJ_XY         projected_loc;
  348         PJ_LP          geodetic_loc;
  349 
  350         projected_loc.x = x[dist*i];
  351         projected_loc.y = y[dist*i];
  352 
  353         if( projected_loc.x == HUGE_VAL )
  354             continue;
  355 
  356         proj_errno_reset( P );
  357         geodetic_loc = pj_inv( projected_loc, P );
  358         if( P->ctx->last_errno != 0 )
  359         {
  360             if( (P->ctx->last_errno != EDOM
  361                     && P->ctx->last_errno != ERANGE)
  362                 && (P->ctx->last_errno > 0
  363                     || P->ctx->last_errno < -44 || n == 1
  364                     || get_transient_error_value(-P->ctx->last_errno) == 0 ) )
  365             {
  366                 return P->ctx->last_errno;
  367             }
  368             else
  369             {
  370                 geodetic_loc.lam = HUGE_VAL;
  371                 geodetic_loc.phi = HUGE_VAL;
  372             }
  373         }
  374 
  375         x[dist*i] = geodetic_loc.lam;
  376         y[dist*i] = geodetic_loc.phi;
  377     }
  378     return 0;
  379 }
  380 
  381 
  382 
  383 /* -------------------------------------------------------------------- */
  384 /*            Adjust for the prime meridian if needed.                  */
  385 /* -------------------------------------------------------------------- */
  386 static int prime_meridian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x) {
  387     int i;
  388     double pm = P->from_greenwich;
  389 
  390     /* Nothing to do? */
  391     if (pm==0.0)
  392         return 0;
  393     if (!(P->is_geocent || P->is_latlong))
  394         return 0;
  395 
  396     if (dir==PJ_FWD)
  397         pm = -pm;
  398 
  399     for (i = 0;  i < n;  i++)
  400         if (x[dist*i] != HUGE_VAL)
  401             x[dist*i] += pm;
  402 
  403     return 0;
  404 }
  405 
  406 
  407 
  408 /* -------------------------------------------------------------------- */
  409 /*            Adjust for vertical scale factor if needed                */
  410 /* -------------------------------------------------------------------- */
  411 static int height_unit (PJ *P, PJ_DIRECTION dir, long n, int dist, double *z) {
  412     int i;
  413     double fac = P->vto_meter;
  414 
  415     if (PJ_FWD==dir)
  416         fac = P->vfr_meter;
  417 
  418     /* Nothing to do? */
  419     if (fac==1.0)
  420         return 0;
  421     if (nullptr==z)
  422         return 0;
  423     if (P->is_latlong)
  424         return 0; /* done in pj_inv3d() / pj_fwd3d() */
  425 
  426     for (i = 0;  i < n;  i++)
  427         if (z[dist*i] != HUGE_VAL )
  428             z[dist*i] *= fac;
  429 
  430     return 0;
  431 }
  432 
  433 
  434 
  435 /* -------------------------------------------------------------------- */
  436 /*           Transform to ellipsoidal heights if needed                 */
  437 /* -------------------------------------------------------------------- */
  438 static int geometric_to_orthometric (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) {
  439     int err;
  440     if (0==P->has_geoid_vgrids)
  441         return 0;
  442     if (z==nullptr)
  443         return PJD_ERR_GEOCENTRIC;
  444     err = pj_apply_vgridshift (P, "sgeoidgrids",
  445               &(P->vgridlist_geoid),
  446               &(P->vgridlist_geoid_count),
  447               dir==PJ_FWD ? 1 : 0, n, dist, x, y, z );
  448     if (err)
  449         return pj_ctx_get_errno(P->ctx);
  450     return 0;
  451 }
  452 
  453 
  454 
  455 /* -------------------------------------------------------------------- */
  456 /*      Convert datums if needed, and possible.                         */
  457 /* -------------------------------------------------------------------- */
  458 static int datum_transform (PJ *P, PJ *Q, long n, int dist, double *x, double *y, double *z) {
  459     if (0==pj_datum_transform (P, Q, n, dist, x, y, z))
  460         return 0;
  461     if (P->ctx->last_errno)
  462         return P->ctx->last_errno;
  463     return Q->ctx->last_errno;
  464 }
  465 
  466 
  467 
  468 
  469 
  470 /* -------------------------------------------------------------------- */
  471 /*      If a wrapping center other than 0 is provided, rewrap around    */
  472 /*      the suggested center (for latlong coordinate systems only).     */
  473 /* -------------------------------------------------------------------- */
  474 static int long_wrap (PJ *P, long n, int dist, double *x) {
  475     long i;
  476 
  477     /* Nothing to do? */
  478     if (P->is_geocent)
  479         return 0;
  480     if (!P->is_long_wrap_set)
  481         return 0;
  482     if (!P->is_latlong)
  483         return 0;
  484 
  485     for (i = 0;  i < n;  i++ ) {
  486         double val = x[dist*i];
  487         if (val == HUGE_VAL)
  488             continue;
  489 
  490         /* Get fast in ] -2 PI, 2 PI [ range */
  491         val = fmod(val, M_TWOPI);
  492         while( val < P->long_wrap_center - M_PI )
  493             val += M_TWOPI;
  494         while( val > P->long_wrap_center + M_PI )
  495             val -= M_TWOPI;
  496         x[dist*i] = val;
  497     }
  498     return 0;
  499 }
  500 
  501 
  502 
  503 /************************************************************************/
  504 /*                            pj_transform()                            */
  505 /*                                                                      */
  506 /*      Currently this function doesn't recognise if two projections    */
  507 /*      are identical (to short circuit reprojection) because it is     */
  508 /*      difficult to compare PJ structures (since there are some        */
  509 /*      projection specific components).                                */
  510 /************************************************************************/
  511 
  512 int pj_transform(
  513     PJ *src, PJ *dst,
  514     long point_count, int point_offset,
  515     double *x, double *y, double *z
  516 ){
  517     int       err;
  518 
  519     src->ctx->last_errno = 0;
  520     dst->ctx->last_errno = 0;
  521 
  522     if( point_offset == 0 )
  523         point_offset = 1;
  524 
  525     /* Bring input to "normal form": longitude, latitude, ellipsoidal height */
  526 
  527     err = adjust_axes (src, PJ_INV, point_count, point_offset, x, y, z);
  528     if (err)
  529         return err;
  530     err = geographic_to_cartesian (src, PJ_INV, point_count, point_offset, x, y, z);
  531     if (err)
  532         return err;
  533     err = projected_to_geographic (src, point_count, point_offset, x, y, z);
  534     if (err)
  535         return err;
  536     err = prime_meridian (src, PJ_INV, point_count, point_offset, x);
  537     if (err)
  538         return err;
  539     err = height_unit (src, PJ_INV, point_count, point_offset, z);
  540     if (err)
  541         return err;
  542     err = geometric_to_orthometric (src, PJ_INV, point_count, point_offset, x, y, z);
  543     if (err)
  544         return err;
  545 
  546     /* At the center of the process we do the datum shift (if needed) */
  547 
  548     err = datum_transform(src, dst, point_count, point_offset, x, y, z );
  549     if (err)
  550         return err;
  551 
  552     /* Now get out on the other side: Bring "normal form" to output form */
  553 
  554     err = geometric_to_orthometric (dst, PJ_FWD, point_count, point_offset, x, y, z);
  555     if (err)
  556         return err;
  557     err = height_unit (dst, PJ_FWD, point_count, point_offset, z);
  558     if (err)
  559         return err;
  560     err = prime_meridian (dst, PJ_FWD, point_count, point_offset, x);
  561     if (err)
  562         return err;
  563     err = geographic_to_cartesian (dst, PJ_FWD, point_count, point_offset, x, y, z);
  564     if (err)
  565         return err;
  566     err = geographic_to_projected (dst, point_count, point_offset, x, y, z);
  567     if (err)
  568         return err;
  569     err = long_wrap (dst, point_count, point_offset, x);
  570     if (err)
  571         return err;
  572     err = adjust_axes (dst, PJ_FWD, point_count, point_offset, x, y, z);
  573     if (err)
  574         return err;
  575 
  576     return 0;
  577 }
  578 
  579 
  580 
  581 /************************************************************************/
  582 /*                     pj_geodetic_to_geocentric()                      */
  583 /************************************************************************/
  584 
  585 int pj_geodetic_to_geocentric( double a, double es,
  586                                long point_count, int point_offset,
  587                                double *x, double *y, double *z )
  588 
  589 {
  590     double b;
  591     int    i;
  592     GeocentricInfo gi;
  593     int ret_errno = 0;
  594 
  595     if( es == 0.0 )
  596         b = a;
  597     else
  598         b = a * sqrt(1-es);
  599 
  600     if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
  601     {
  602         return PJD_ERR_GEOCENTRIC;
  603     }
  604 
  605     for( i = 0; i < point_count; i++ )
  606     {
  607         long io = i * point_offset;
  608 
  609         if( x[io] == HUGE_VAL  )
  610             continue;
  611 
  612         if( pj_Convert_Geodetic_To_Geocentric( &gi, y[io], x[io], z[io],
  613                                                x+io, y+io, z+io ) != 0 )
  614         {
  615             ret_errno = PJD_ERR_LAT_OR_LON_EXCEED_LIMIT;
  616             x[io] = y[io] = HUGE_VAL;
  617             /* but keep processing points! */
  618         }
  619     }
  620 
  621     return ret_errno;
  622 }
  623 
  624 /************************************************************************/
  625 /*                     pj_geocentric_to_geodetic()                      */
  626 /************************************************************************/
  627 
  628 int pj_geocentric_to_geodetic( double a, double es,
  629                                long point_count, int point_offset,
  630                                double *x, double *y, double *z )
  631 
  632 {
  633     double b;
  634     int    i;
  635     GeocentricInfo gi;
  636 
  637     if( es == 0.0 )
  638         b = a;
  639     else
  640         b = a * sqrt(1-es);
  641 
  642     if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
  643     {
  644         return PJD_ERR_GEOCENTRIC;
  645     }
  646 
  647     for( i = 0; i < point_count; i++ )
  648     {
  649         long io = i * point_offset;
  650 
  651         if( x[io] == HUGE_VAL )
  652             continue;
  653 
  654         pj_Convert_Geocentric_To_Geodetic( &gi, x[io], y[io], z[io],
  655                                            y+io, x+io, z+io );
  656     }
  657 
  658     return 0;
  659 }
  660 
  661 /************************************************************************/
  662 /*                         pj_compare_datums()                          */
  663 /*                                                                      */
  664 /*      Returns TRUE if the two datums are identical, otherwise         */
  665 /*      FALSE.                                                          */
  666 /************************************************************************/
  667 
  668 int pj_compare_datums( PJ *srcdefn, PJ *dstdefn )
  669 
  670 {
  671     if( srcdefn->datum_type != dstdefn->datum_type )
  672     {
  673         return 0;
  674     }
  675     else if( srcdefn->a_orig != dstdefn->a_orig
  676              || ABS(srcdefn->es_orig - dstdefn->es_orig) > 0.000000000050 )
  677     {
  678         /* the tolerance for es is to ensure that GRS80 and WGS84 are
  679            considered identical */
  680         return 0;
  681     }
  682     else if( srcdefn->datum_type == PJD_3PARAM )
  683     {
  684         return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
  685                 && srcdefn->datum_params[1] == dstdefn->datum_params[1]
  686                 && srcdefn->datum_params[2] == dstdefn->datum_params[2]);
  687     }
  688     else if( srcdefn->datum_type == PJD_7PARAM )
  689     {
  690         return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
  691                 && srcdefn->datum_params[1] == dstdefn->datum_params[1]
  692                 && srcdefn->datum_params[2] == dstdefn->datum_params[2]
  693                 && srcdefn->datum_params[3] == dstdefn->datum_params[3]
  694                 && srcdefn->datum_params[4] == dstdefn->datum_params[4]
  695                 && srcdefn->datum_params[5] == dstdefn->datum_params[5]
  696                 && srcdefn->datum_params[6] == dstdefn->datum_params[6]);
  697     }
  698     else if( srcdefn->datum_type == PJD_GRIDSHIFT )
  699     {
  700         const char* srcnadgrids =
  701             pj_param(srcdefn->ctx, srcdefn->params,"snadgrids").s;
  702         const char* dstnadgrids =
  703             pj_param(dstdefn->ctx, dstdefn->params,"snadgrids").s;
  704         return srcnadgrids != nullptr && dstnadgrids != nullptr &&
  705                strcmp( srcnadgrids, dstnadgrids ) == 0;
  706     }
  707     else
  708         return 1;
  709 }
  710 
  711 /************************************************************************/
  712 /*                       pj_geocentic_to_wgs84()                        */
  713 /************************************************************************/
  714 
  715 static
  716 int pj_geocentric_to_wgs84( PJ *defn,
  717                             long point_count, int point_offset,
  718                             double *x, double *y, double *z )
  719 
  720 {
  721     int       i;
  722 
  723     if( defn->datum_type == PJD_3PARAM )
  724     {
  725         for( i = 0; i < point_count; i++ )
  726         {
  727             long io = i * point_offset;
  728 
  729             if( x[io] == HUGE_VAL )
  730                 continue;
  731 
  732             x[io] = x[io] + Dx_BF;
  733             y[io] = y[io] + Dy_BF;
  734             z[io] = z[io] + Dz_BF;
  735         }
  736     }
  737     else if( defn->datum_type == PJD_7PARAM )
  738     {
  739         for( i = 0; i < point_count; i++ )
  740         {
  741             long io = i * point_offset;
  742             double x_out, y_out, z_out;
  743 
  744             if( x[io] == HUGE_VAL )
  745                 continue;
  746 
  747             x_out = M_BF*(       x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
  748             y_out = M_BF*( Rz_BF*x[io] +       y[io] - Rx_BF*z[io]) + Dy_BF;
  749             z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] +       z[io]) + Dz_BF;
  750 
  751             x[io] = x_out;
  752             y[io] = y_out;
  753             z[io] = z_out;
  754         }
  755     }
  756 
  757     return 0;
  758 }
  759 
  760 /************************************************************************/
  761 /*                      pj_geocentic_from_wgs84()                       */
  762 /************************************************************************/
  763 
  764 static
  765 int pj_geocentric_from_wgs84( PJ *defn,
  766                               long point_count, int point_offset,
  767                               double *x, double *y, double *z )
  768 
  769 {
  770     int       i;
  771 
  772     if( defn->datum_type == PJD_3PARAM )
  773     {
  774         for( i = 0; i < point_count; i++ )
  775         {
  776             long io = i * point_offset;
  777 
  778             if( x[io] == HUGE_VAL )
  779                 continue;
  780 
  781             x[io] = x[io] - Dx_BF;
  782             y[io] = y[io] - Dy_BF;
  783             z[io] = z[io] - Dz_BF;
  784         }
  785     }
  786     else if( defn->datum_type == PJD_7PARAM )
  787     {
  788         for( i = 0; i < point_count; i++ )
  789         {
  790             long io = i * point_offset;
  791             double x_tmp, y_tmp, z_tmp;
  792 
  793             if( x[io] == HUGE_VAL )
  794                 continue;
  795 
  796             x_tmp = (x[io] - Dx_BF) / M_BF;
  797             y_tmp = (y[io] - Dy_BF) / M_BF;
  798             z_tmp = (z[io] - Dz_BF) / M_BF;
  799 
  800             x[io] =        x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp;
  801             y[io] = -Rz_BF*x_tmp +       y_tmp + Rx_BF*z_tmp;
  802             z[io] =  Ry_BF*x_tmp - Rx_BF*y_tmp +       z_tmp;
  803         }
  804     }
  805 
  806     return 0;
  807 }
  808 
  809 /************************************************************************/
  810 /*                         pj_datum_transform()                         */
  811 /*                                                                      */
  812 /*      The input should be long/lat/z coordinates in radians in the    */
  813 /*      source datum, and the output should be long/lat/z               */
  814 /*      coordinates in radians in the destination datum.                */
  815 /************************************************************************/
  816 
  817 int pj_datum_transform( PJ *src, PJ *dst,
  818                         long point_count, int point_offset,
  819                         double *x, double *y, double *z )
  820 
  821 {
  822     double      src_a, src_es, dst_a, dst_es;
  823     int         z_is_temp = FALSE;
  824 
  825 /* -------------------------------------------------------------------- */
  826 /*      We cannot do any meaningful datum transformation if either      */
  827 /*      the source or destination are of an unknown datum type          */
  828 /*      (ie. only a +ellps declaration, no +datum).  This is new        */
  829 /*      behavior for PROJ 4.6.0.                                        */
  830 /* -------------------------------------------------------------------- */
  831     if( src->datum_type == PJD_UNKNOWN
  832         || dst->datum_type == PJD_UNKNOWN )
  833         return 0;
  834 
  835 /* -------------------------------------------------------------------- */
  836 /*      Short cut if the datums are identical.                          */
  837 /* -------------------------------------------------------------------- */
  838     if( pj_compare_datums( src, dst ) )
  839         return 0;
  840 
  841     src_a = src->a_orig;
  842     src_es = src->es_orig;
  843 
  844     dst_a = dst->a_orig;
  845     dst_es = dst->es_orig;
  846 
  847 /* -------------------------------------------------------------------- */
  848 /*      Create a temporary Z array if one is not provided.              */
  849 /* -------------------------------------------------------------------- */
  850     if( z == nullptr )
  851     {
  852         size_t  bytes = sizeof(double) * point_count * point_offset;
  853         z = (double *) pj_malloc(bytes);
  854         memset( z, 0, bytes );
  855         z_is_temp = TRUE;
  856     }
  857 
  858 #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; }}
  859 
  860 /* -------------------------------------------------------------------- */
  861 /*  If this datum requires grid shifts, then apply it to geodetic   */
  862 /*      coordinates.                                                    */
  863 /* -------------------------------------------------------------------- */
  864     if( src->datum_type == PJD_GRIDSHIFT )
  865     {
  866         pj_apply_gridshift_2( src, 0, point_count, point_offset, x, y, z );
  867         CHECK_RETURN(src);
  868 
  869         src_a = SRS_WGS84_SEMIMAJOR;
  870         src_es = SRS_WGS84_ESQUARED;
  871     }
  872 
  873     if( dst->datum_type == PJD_GRIDSHIFT )
  874     {
  875         dst_a = SRS_WGS84_SEMIMAJOR;
  876         dst_es = SRS_WGS84_ESQUARED;
  877     }
  878 
  879 /* ==================================================================== */
  880 /*      Do we need to go through geocentric coordinates?                */
  881 /* ==================================================================== */
  882     if( src_es != dst_es || src_a != dst_a
  883         || src->datum_type == PJD_3PARAM
  884         || src->datum_type == PJD_7PARAM
  885         || dst->datum_type == PJD_3PARAM
  886         || dst->datum_type == PJD_7PARAM)
  887     {
  888 /* -------------------------------------------------------------------- */
  889 /*      Convert to geocentric coordinates.                              */
  890 /* -------------------------------------------------------------------- */
  891         src->ctx->last_errno =
  892             pj_geodetic_to_geocentric( src_a, src_es,
  893                                        point_count, point_offset, x, y, z );
  894         CHECK_RETURN(src);
  895 
  896 /* -------------------------------------------------------------------- */
  897 /*      Convert between datums.                                         */
  898 /* -------------------------------------------------------------------- */
  899         if( src->datum_type == PJD_3PARAM
  900             || src->datum_type == PJD_7PARAM )
  901         {
  902             pj_geocentric_to_wgs84( src, point_count, point_offset,x,y,z);
  903             CHECK_RETURN(src);
  904         }
  905 
  906         if( dst->datum_type == PJD_3PARAM
  907             || dst->datum_type == PJD_7PARAM )
  908         {
  909             pj_geocentric_from_wgs84( dst, point_count,point_offset,x,y,z);
  910             CHECK_RETURN(dst);
  911         }
  912 
  913 /* -------------------------------------------------------------------- */
  914 /*      Convert back to geodetic coordinates.                           */
  915 /* -------------------------------------------------------------------- */
  916         dst->ctx->last_errno =
  917             pj_geocentric_to_geodetic( dst_a, dst_es,
  918                                        point_count, point_offset, x, y, z );
  919         CHECK_RETURN(dst);
  920     }
  921 
  922 /* -------------------------------------------------------------------- */
  923 /*      Apply grid shift to destination if required.                    */
  924 /* -------------------------------------------------------------------- */
  925     if( dst->datum_type == PJD_GRIDSHIFT )
  926     {
  927         pj_apply_gridshift_2( dst, 1, point_count, point_offset, x, y, z );
  928         CHECK_RETURN(dst);
  929     }
  930 
  931     if( z_is_temp )
  932         pj_dalloc( z );
  933 
  934     return 0;
  935 }
  936 
  937 /************************************************************************/
  938 /*                           adjust_axis()                           */
  939 /*                                                                      */
  940 /*      Normalize or de-normalized the x/y/z axes.  The normal form     */
  941 /*      is "enu" (easting, northing, up).                               */
  942 /************************************************************************/
  943 static int adjust_axis( projCtx ctx,
  944                            const char *axis, int denormalize_flag,
  945                            long point_count, int point_offset,
  946                            double *x, double *y, double *z )
  947 
  948 {
  949     double x_in, y_in, z_in = 0.0;
  950     int i, i_axis;
  951 
  952     if( !denormalize_flag )
  953     {
  954         for( i = 0; i < point_count; i++ )
  955         {
  956             x_in = x[point_offset*i];
  957             y_in = y[point_offset*i];
  958             if( z )
  959                 z_in = z[point_offset*i];
  960 
  961             for( i_axis = 0; i_axis < 3; i_axis++ )
  962             {
  963                 double value;
  964 
  965                 if( i_axis == 0 )
  966                     value = x_in;
  967                 else if( i_axis == 1 )
  968                     value = y_in;
  969                 else
  970                     value = z_in;
  971 
  972                 switch( axis[i_axis] )
  973                 {
  974                   case 'e':
  975                     x[point_offset*i] = value;
  976                     break;
  977                   case 'w':
  978                     x[point_offset*i] = -value;
  979                     break;
  980                   case 'n':
  981                     y[point_offset*i] = value;
  982                     break;
  983                   case 's':
  984                     y[point_offset*i] = -value;
  985                     break;
  986                   case 'u':
  987                     if( z )
  988                         z[point_offset*i] = value;
  989                     break;
  990                   case 'd':
  991                     if( z )
  992                         z[point_offset*i] = -value;
  993                     break;
  994                   default:
  995                     pj_ctx_set_errno( ctx, PJD_ERR_AXIS );
  996                     return PJD_ERR_AXIS;
  997                 }
  998             } /* i_axis */
  999         } /* i (point) */
 1000     }
 1001 
 1002     else /* denormalize */
 1003     {
 1004         for( i = 0; i < point_count; i++ )
 1005         {
 1006             x_in = x[point_offset*i];
 1007             y_in = y[point_offset*i];
 1008             if( z )
 1009                 z_in = z[point_offset*i];
 1010 
 1011             for( i_axis = 0; i_axis < 3; i_axis++ )
 1012             {
 1013                 double *target;
 1014 
 1015                 if( i_axis == 2 && z == nullptr )
 1016                     continue;
 1017 
 1018                 if( i_axis == 0 )
 1019                     target = x;
 1020                 else if( i_axis == 1 )
 1021                     target = y;
 1022                 else
 1023                     target = z;
 1024 
 1025                 switch( axis[i_axis] )
 1026                 {
 1027                   case 'e':
 1028                     target[point_offset*i] = x_in; break;
 1029                   case 'w':
 1030                     target[point_offset*i] = -x_in; break;
 1031                   case 'n':
 1032                     target[point_offset*i] = y_in; break;
 1033                   case 's':
 1034                     target[point_offset*i] = -y_in; break;
 1035                   case 'u':
 1036                     target[point_offset*i] = z_in; break;
 1037                   case 'd':
 1038                     target[point_offset*i] = -z_in; break;
 1039                   default:
 1040                     pj_ctx_set_errno( ctx, PJD_ERR_AXIS );
 1041                     return PJD_ERR_AXIS;
 1042                 }
 1043             } /* i_axis */
 1044         } /* i (point) */
 1045     }
 1046 
 1047     return 0;
 1048 }