PageRenderTime 36ms CodeModel.GetById 2ms app.highlight 27ms RepoModel.GetById 1ms app.codeStats 1ms

/Proj4/pj_transform.c

http://github.com/route-me/route-me
C | 740 lines | 409 code | 100 blank | 231 comment | 188 complexity | 18919984cdbeb6f07b2e8eac21735a21 MD5 | raw file
  1/******************************************************************************
  2 * $Id: pj_transform.c,v 1.24 2007/12/03 15:48:20 fwarmerdam Exp $
  3 *
  4 * Project:  PROJ.4
  5 * Purpose:  Perform overall coordinate system to coordinate system 
  6 *           transformations (pj_transform() function) including reprojection
  7 *           and datum shifting.
  8 * Author:   Frank Warmerdam, warmerdam@pobox.com
  9 *
 10 ******************************************************************************
 11 * Copyright (c) 2000, Frank Warmerdam
 12 *
 13 * Permission is hereby granted, free of charge, to any person obtaining a
 14 * copy of this software and associated documentation files (the "Software"),
 15 * to deal in the Software without restriction, including without limitation
 16 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 17 * and/or sell copies of the Software, and to permit persons to whom the
 18 * Software is furnished to do so, subject to the following conditions:
 19 *
 20 * The above copyright notice and this permission notice shall be included
 21 * in all copies or substantial portions of the Software.
 22 *
 23 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 24 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 25 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 26 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 27 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 28 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 29 * DEALINGS IN THE SOFTWARE.
 30 ******************************************************************************
 31 *
 32 * $Log: pj_transform.c,v $
 33 * Revision 1.24  2007/12/03 15:48:20  fwarmerdam
 34 * Improve WGS84 ES precision to avoid unnecesary transformation (#1531)
 35 *
 36 * Revision 1.23  2007/11/26 00:21:59  fwarmerdam
 37 * Modified PJ structure to hold a_orig, es_orig, ellipsoid definition before
 38 * adjustment for spherical projections.
 39 * Modified pj_datum_transform() to use the original ellipsoid parameters,
 40 * not the ones adjusted for spherical projections.
 41 * Modified pj_datum_transform() to not attempt any datum shift via
 42 * geocentric coordinates if the source *or* destination are raw ellipsoids
 43 * (ie. PJD_UNKNOWN).  All per PROJ bug #1602, GDAL bug #2025.
 44 *
 45 * Revision 1.22  2007/09/11 20:32:25  fwarmerdam
 46 * mark the transient error array const
 47 *
 48 * Revision 1.21  2007/09/11 20:19:36  fwarmerdam
 49 * avoid use of static variables to make reentrant
 50 *
 51 * Revision 1.20  2006/10/12 21:04:39  fwarmerdam
 52 * Added experimental +lon_wrap argument to set a "center point" for
 53 * longitude wrapping of longitude values coming out of pj_transform().
 54 *
 55 * Revision 1.19  2006/05/10 19:23:47  fwarmerdam
 56 * Don't apply to_meter in pj_transform() if the value is HUGE_VAL.
 57 *
 58 * Revision 1.18  2006/05/01 21:13:54  fwarmerdam
 59 * make out of range errors in geodetic to geocentric a transient error
 60 *
 61 * Revision 1.17  2006/03/20 17:54:34  fwarmerdam
 62 * pj_geodetic_to_geocentric returns -14 now for lat out of range
 63 *
 64 * Revision 1.16  2006/02/17 02:26:14  fwarmerdam
 65 * ERANGE/EDOM treated as transient errors
 66 *
 67 * Revision 1.15  2005/12/04 14:47:37  fwarmerdam
 68 * use symbolic names as per patch from Martin Vermeer
 69 *
 70 * Revision 1.14  2004/11/05 06:05:11  fwarmerdam
 71 * Fixed pj_geocentric_to_geodetic() to not try and process HUGE_VAL values
 72 * (those that have failed some previous transform step).  Related to bug:5B
 73 *     http://bugzilla.remotesensing.org/show_bug.cgi?id=642
 74 *
 75 * Revision 1.13  2004/10/25 15:34:36  fwarmerdam
 76 * make names of geodetic funcs from geotrans unique
 77 *
 78 * Revision 1.12  2004/05/03 19:45:23  warmerda
 79 * Altered so that raw ellpses are treated as a essentially having a
 80 * +towgs84=0,0,0 specification so ellpisoid shifts are applied.
 81 * Fixed so that prime meridian shifts are applied if the coordinate system is
 82 * not lat/long (ie. if it is projected).  This fixes:
 83 * http://bugzilla.remotesensing.org/show_bug.cgi?id=510
 84 *
 85 * Revision 1.11  2004/01/24 09:37:19  warmerda
 86 * pj_transform() will no longer return an error code if any of the points are
 87 * transformable.  In this case the application is expected to check for
 88 * HUGE_VAL to identify failed points.
 89 * As part of the implementation, I added a list of pj_errno values that
 90 * are transient (ie per-point) rather than indicating global failure for the
 91 * coordinate system definition.  We use this in deciding which pj_fwd and
 92 * pj_inv error codes are really fatal and should be reported.
 93 *
 94 * Revision 1.10  2003/08/21 02:09:06  warmerda
 95 * added a bunch of HUGE_VAL checking
 96 *
 97 * Revision 1.9  2003/03/26 16:52:30  warmerda
 98 * added check that an inverse transformation func exists
 99 *
100 * Revision 1.8  2002/12/14 20:35:43  warmerda
101 * implement units support for geocentric coordinates
102 *
103 * Revision 1.7  2002/12/14 20:14:35  warmerda
104 * added geocentric support
105 *
106 * Revision 1.6  2002/12/09 16:01:02  warmerda
107 * added prime meridian support
108 *
109 * Revision 1.5  2002/12/01 19:25:26  warmerda
110 * applied fix for 7 param shift in pj_geocentric_from_wgs84, see bug 194
111 *
112 * Revision 1.4  2002/02/15 14:30:36  warmerda
113 * provide default Z array if none passed in in pj_datum_transform()
114 *
115 * Revision 1.3  2001/04/04 21:13:21  warmerda
116 * do arcsecond/radian and ppm datum parm transformation in pj_set_datum()
117 *
118 * Revision 1.2  2001/04/04 16:08:08  warmerda
119 * rewrote 7 param datum shift to match EPSG:9606, now works with example
120 *
121 * Revision 1.1  2000/07/06 23:32:27  warmerda
122 * New
123 *
124 */
125
126#include "projects.h"
127#include <string.h>
128#include <math.h>
129#include "geocent.h"
130
131PJ_CVSID("$Id: pj_transform.c,v 1.24 2007/12/03 15:48:20 fwarmerdam Exp $");
132
133#ifndef SRS_WGS84_SEMIMAJOR
134#define SRS_WGS84_SEMIMAJOR 6378137.0
135#endif
136
137#ifndef SRS_WGS84_ESQUARED
138#define SRS_WGS84_ESQUARED 0.0066943799901413165
139#endif
140
141#define Dx_BF (defn->datum_params[0])
142#define Dy_BF (defn->datum_params[1])
143#define Dz_BF (defn->datum_params[2])
144#define Rx_BF (defn->datum_params[3])
145#define Ry_BF (defn->datum_params[4])
146#define Rz_BF (defn->datum_params[5])
147#define M_BF  (defn->datum_params[6])
148
149/* 
150** This table is intended to indicate for any given error code in 
151** the range 0 to -44, whether that error will occur for all locations (ie.
152** it is a problem with the coordinate system as a whole) in which case the
153** value would be 0, or if the problem is with the point being transformed
154** in which case the value is 1. 
155**
156** At some point we might want to move this array in with the error message
157** list or something, but while experimenting with it this should be fine. 
158*/
159
160static const int transient_error[45] = {
161    /*             0  1  2  3  4  5  6  7  8  9   */
162    /* 0 to 9 */   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
163    /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,  
164    /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
165    /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
166    /* 40 to 44 */ 0, 0, 0, 0, 0 };
167
168/************************************************************************/
169/*                            pj_transform()                            */
170/*                                                                      */
171/*      Currently this function doesn't recognise if two projections    */
172/*      are identical (to short circuit reprojection) because it is     */
173/*      difficult to compare PJ structures (since there are some        */
174/*      projection specific components).                                */
175/************************************************************************/
176
177int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
178                  double *x, double *y, double *z )
179
180{
181    long      i;
182    int       need_datum_shift;
183
184    pj_errno = 0;
185
186    if( point_offset == 0 )
187        point_offset = 1;
188
189/* -------------------------------------------------------------------- */
190/*      Transform geocentric source coordinates to lat/long.            */
191/* -------------------------------------------------------------------- */
192    if( srcdefn->is_geocent )
193    {
194        if( z == NULL )
195        {
196            pj_errno = PJD_ERR_GEOCENTRIC;
197            return PJD_ERR_GEOCENTRIC;
198        }
199
200        if( srcdefn->to_meter != 1.0 )
201        {
202            for( i = 0; i < point_count; i++ )
203            {
204                if( x[point_offset*i] != HUGE_VAL )
205                {
206                    x[point_offset*i] *= srcdefn->to_meter;
207                    y[point_offset*i] *= srcdefn->to_meter;
208                }
209            }
210        }
211
212        if( pj_geocentric_to_geodetic( srcdefn->a_orig, srcdefn->es_orig,
213                                       point_count, point_offset, 
214                                       x, y, z ) != 0) 
215            return pj_errno;
216    }
217
218/* -------------------------------------------------------------------- */
219/*      Transform source points to lat/long, if they aren't             */
220/*      already.                                                        */
221/* -------------------------------------------------------------------- */
222    else if( !srcdefn->is_latlong )
223    {
224        if( srcdefn->inv == NULL )
225        {
226            pj_errno = -17; /* this isn't correct, we need a no inverse err */
227            if( getenv( "PROJ_DEBUG" ) != NULL )
228            {
229                fprintf( stderr, 
230                       "pj_transform(): source projection not invertable\n" );
231            }
232            return pj_errno;
233        }
234
235        for( i = 0; i < point_count; i++ )
236        {
237            XY         projected_loc;
238            LP	       geodetic_loc;
239
240            projected_loc.u = x[point_offset*i];
241            projected_loc.v = y[point_offset*i];
242
243            if( projected_loc.u == HUGE_VAL )
244                continue;
245
246            geodetic_loc = pj_inv( projected_loc, srcdefn );
247            if( pj_errno != 0 )
248            {
249                if( (pj_errno != 33 /*EDOM*/ && pj_errno != 34 /*ERANGE*/ )
250                    && (pj_errno > 0 || pj_errno < -44 || point_count == 1
251                        || transient_error[-pj_errno] == 0 ) )
252                    return pj_errno;
253                else
254                {
255                    geodetic_loc.u = HUGE_VAL;
256                    geodetic_loc.v = HUGE_VAL;
257                }
258            }
259
260            x[point_offset*i] = geodetic_loc.u;
261            y[point_offset*i] = geodetic_loc.v;
262        }
263    }
264/* -------------------------------------------------------------------- */
265/*      But if they are already lat long, adjust for the prime          */
266/*      meridian if there is one in effect.                             */
267/* -------------------------------------------------------------------- */
268    if( srcdefn->from_greenwich != 0.0 )
269    {
270        for( i = 0; i < point_count; i++ )
271        {
272            if( x[point_offset*i] != HUGE_VAL )
273                x[point_offset*i] += srcdefn->from_greenwich;
274        }
275    }
276
277/* -------------------------------------------------------------------- */
278/*      Convert datums if needed, and possible.                         */
279/* -------------------------------------------------------------------- */
280    if( pj_datum_transform( srcdefn, dstdefn, point_count, point_offset, 
281                            x, y, z ) != 0 )
282        return pj_errno;
283
284/* -------------------------------------------------------------------- */
285/*      But if they are staying lat long, adjust for the prime          */
286/*      meridian if there is one in effect.                             */
287/* -------------------------------------------------------------------- */
288    if( dstdefn->from_greenwich != 0.0 )
289    {
290        for( i = 0; i < point_count; i++ )
291        {
292            if( x[point_offset*i] != HUGE_VAL )
293                x[point_offset*i] -= dstdefn->from_greenwich;
294        }
295    }
296
297
298/* -------------------------------------------------------------------- */
299/*      Transform destination latlong to geocentric if required.        */
300/* -------------------------------------------------------------------- */
301    if( dstdefn->is_geocent )
302    {
303        if( z == NULL )
304        {
305            pj_errno = PJD_ERR_GEOCENTRIC;
306            return PJD_ERR_GEOCENTRIC;
307        }
308
309        pj_geodetic_to_geocentric( dstdefn->a_orig, dstdefn->es_orig,
310                                   point_count, point_offset, x, y, z );
311
312        if( dstdefn->fr_meter != 1.0 )
313        {
314            for( i = 0; i < point_count; i++ )
315            {
316                if( x[point_offset*i] != HUGE_VAL )
317                {
318                    x[point_offset*i] *= dstdefn->fr_meter;
319                    y[point_offset*i] *= dstdefn->fr_meter;
320                }
321            }
322        }
323    }
324
325/* -------------------------------------------------------------------- */
326/*      Transform destination points to projection coordinates, if      */
327/*      desired.                                                        */
328/* -------------------------------------------------------------------- */
329    else if( !dstdefn->is_latlong )
330    {
331        for( i = 0; i < point_count; i++ )
332        {
333            XY         projected_loc;
334            LP	       geodetic_loc;
335
336            geodetic_loc.u = x[point_offset*i];
337            geodetic_loc.v = y[point_offset*i];
338
339            if( geodetic_loc.u == HUGE_VAL )
340                continue;
341
342            projected_loc = pj_fwd( geodetic_loc, dstdefn );
343            if( pj_errno != 0 )
344            {
345                if( (pj_errno != 33 /*EDOM*/ && pj_errno != 34 /*ERANGE*/ )
346                    && (pj_errno > 0 || pj_errno < -44 || point_count == 1
347                        || transient_error[-pj_errno] == 0 ) )
348                    return pj_errno;
349                else
350                {
351                    projected_loc.u = HUGE_VAL;
352                    projected_loc.v = HUGE_VAL;
353                }
354            }
355
356            x[point_offset*i] = projected_loc.u;
357            y[point_offset*i] = projected_loc.v;
358        }
359    }
360
361/* -------------------------------------------------------------------- */
362/*      If a wrapping center other than 0 is provided, rewrap around    */
363/*      the suggested center (for latlong coordinate systems only).     */
364/* -------------------------------------------------------------------- */
365    else if( dstdefn->is_latlong && dstdefn->long_wrap_center != 0 )
366    {
367        for( i = 0; i < point_count; i++ )
368        {
369            if( x[point_offset*i] == HUGE_VAL )
370                continue;
371
372            while( x[point_offset*i] < dstdefn->long_wrap_center - HALFPI )
373                x[point_offset*i] += PI;
374            while( x[point_offset*i] > dstdefn->long_wrap_center + HALFPI )
375                x[point_offset*i] -= PI;
376        }
377    }
378
379    return 0;
380}
381
382/************************************************************************/
383/*                     pj_geodetic_to_geocentric()                      */
384/************************************************************************/
385
386int pj_geodetic_to_geocentric( double a, double es, 
387                               long point_count, int point_offset,
388                               double *x, double *y, double *z )
389
390{
391    double b;
392    int    i;
393    GeocentricInfo gi;
394
395    pj_errno = 0;
396
397    if( es == 0.0 )
398        b = a;
399    else
400        b = a * sqrt(1-es);
401
402    if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
403    {
404        pj_errno = PJD_ERR_GEOCENTRIC;
405        return pj_errno;
406    }
407
408    for( i = 0; i < point_count; i++ )
409    {
410        long io = i * point_offset;
411
412        if( x[io] == HUGE_VAL  )
413            continue;
414
415        if( pj_Convert_Geodetic_To_Geocentric( &gi, y[io], x[io], z[io], 
416                                               x+io, y+io, z+io ) != 0 )
417        {
418            pj_errno = -14;
419            x[io] = y[io] = HUGE_VAL;
420            /* but keep processing points! */
421        }
422    }
423
424    return pj_errno;
425}
426
427/************************************************************************/
428/*                     pj_geodetic_to_geocentric()                      */
429/************************************************************************/
430
431int pj_geocentric_to_geodetic( double a, double es, 
432                               long point_count, int point_offset,
433                               double *x, double *y, double *z )
434
435{
436    double b;
437    int    i;
438    GeocentricInfo gi;
439
440    if( es == 0.0 )
441        b = a;
442    else
443        b = a * sqrt(1-es);
444
445    if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
446    {
447        pj_errno = PJD_ERR_GEOCENTRIC;
448        return pj_errno;
449    }
450
451    for( i = 0; i < point_count; i++ )
452    {
453        long io = i * point_offset;
454
455        if( x[io] == HUGE_VAL )
456            continue;
457
458        pj_Convert_Geocentric_To_Geodetic( &gi, x[io], y[io], z[io], 
459                                           y+io, x+io, z+io );
460    }
461
462    return 0;
463}
464
465/************************************************************************/
466/*                         pj_compare_datums()                          */
467/*                                                                      */
468/*      Returns TRUE if the two datums are identical, otherwise         */
469/*      FALSE.                                                          */
470/************************************************************************/
471
472int pj_compare_datums( PJ *srcdefn, PJ *dstdefn )
473
474{
475    if( srcdefn->datum_type != dstdefn->datum_type )
476    {
477        return 0;
478    }
479    else if( srcdefn->a_orig != dstdefn->a_orig 
480             || ABS(srcdefn->es_orig - dstdefn->es_orig) > 0.000000000050 )
481    {
482        /* the tolerence for es is to ensure that GRS80 and WGS84 are
483           considered identical */
484        return 0;
485    }
486    else if( srcdefn->datum_type == PJD_3PARAM )
487    {
488        return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
489                && srcdefn->datum_params[1] == dstdefn->datum_params[1]
490                && srcdefn->datum_params[2] == dstdefn->datum_params[2]);
491    }
492    else if( srcdefn->datum_type == PJD_7PARAM )
493    {
494        return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
495                && srcdefn->datum_params[1] == dstdefn->datum_params[1]
496                && srcdefn->datum_params[2] == dstdefn->datum_params[2]
497                && srcdefn->datum_params[3] == dstdefn->datum_params[3]
498                && srcdefn->datum_params[4] == dstdefn->datum_params[4]
499                && srcdefn->datum_params[5] == dstdefn->datum_params[5]
500                && srcdefn->datum_params[6] == dstdefn->datum_params[6]);
501    }
502    else if( srcdefn->datum_type == PJD_GRIDSHIFT )
503    {
504        return strcmp( pj_param(srcdefn->params,"snadgrids").s,
505                       pj_param(dstdefn->params,"snadgrids").s ) == 0;
506    }
507    else
508        return 1;
509}
510
511/************************************************************************/
512/*                       pj_geocentic_to_wgs84()                        */
513/************************************************************************/
514
515int pj_geocentric_to_wgs84( PJ *defn, 
516                            long point_count, int point_offset,
517                            double *x, double *y, double *z )
518
519{
520    int       i;
521
522    pj_errno = 0;
523
524    if( defn->datum_type == PJD_3PARAM )
525    {
526        for( i = 0; i < point_count; i++ )
527        {
528            long io = i * point_offset;
529            
530            if( x[io] == HUGE_VAL )
531                continue;
532
533            x[io] = x[io] + Dx_BF;
534            y[io] = y[io] + Dy_BF;
535            z[io] = z[io] + Dz_BF;
536        }
537    }
538    else if( defn->datum_type == PJD_7PARAM )
539    {
540        for( i = 0; i < point_count; i++ )
541        {
542            long io = i * point_offset;
543            double x_out, y_out, z_out;
544
545            if( x[io] == HUGE_VAL )
546                continue;
547
548            x_out = M_BF*(       x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
549            y_out = M_BF*( Rz_BF*x[io] +       y[io] - Rx_BF*z[io]) + Dy_BF;
550            z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] +       z[io]) + Dz_BF;
551
552            x[io] = x_out;
553            y[io] = y_out;
554            z[io] = z_out;
555        }
556    }
557
558    return 0;
559}
560
561/************************************************************************/
562/*                      pj_geocentic_from_wgs84()                       */
563/************************************************************************/
564
565int pj_geocentric_from_wgs84( PJ *defn, 
566                              long point_count, int point_offset,
567                              double *x, double *y, double *z )
568
569{
570    int       i;
571
572    pj_errno = 0;
573
574    if( defn->datum_type == PJD_3PARAM )
575    {
576        for( i = 0; i < point_count; i++ )
577        {
578            long io = i * point_offset;
579
580            if( x[io] == HUGE_VAL )
581                continue;
582            
583            x[io] = x[io] - Dx_BF;
584            y[io] = y[io] - Dy_BF;
585            z[io] = z[io] - Dz_BF;
586        }
587    }
588    else if( defn->datum_type == PJD_7PARAM )
589    {
590        for( i = 0; i < point_count; i++ )
591        {
592            long io = i * point_offset;
593            double x_tmp, y_tmp, z_tmp;
594
595            if( x[io] == HUGE_VAL )
596                continue;
597
598            x_tmp = (x[io] - Dx_BF) / M_BF;
599            y_tmp = (y[io] - Dy_BF) / M_BF;
600            z_tmp = (z[io] - Dz_BF) / M_BF;
601
602            x[io] =        x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp;
603            y[io] = -Rz_BF*x_tmp +       y_tmp + Rx_BF*z_tmp;
604            z[io] =  Ry_BF*x_tmp - Rx_BF*y_tmp +       z_tmp;
605        }
606    }
607
608    return 0;
609}
610
611/************************************************************************/
612/*                         pj_datum_transform()                         */
613/*                                                                      */
614/*      The input should be long/lat/z coordinates in radians in the    */
615/*      source datum, and the output should be long/lat/z               */
616/*      coordinates in radians in the destination datum.                */
617/************************************************************************/
618
619int pj_datum_transform( PJ *srcdefn, PJ *dstdefn, 
620                        long point_count, int point_offset,
621                        double *x, double *y, double *z )
622
623{
624    double      src_a, src_es, dst_a, dst_es;
625    int         z_is_temp = FALSE;
626
627    pj_errno = 0;
628
629/* -------------------------------------------------------------------- */
630/*      We cannot do any meaningful datum transformation if either      */
631/*      the source or destination are of an unknown datum type          */
632/*      (ie. only a +ellps declaration, no +datum).  This is new        */
633/*      behavior for PROJ 4.6.0.                                        */
634/* -------------------------------------------------------------------- */
635    if( srcdefn->datum_type == PJD_UNKNOWN
636        || dstdefn->datum_type == PJD_UNKNOWN )
637        return 0;
638
639/* -------------------------------------------------------------------- */
640/*      Short cut if the datums are identical.                          */
641/* -------------------------------------------------------------------- */
642    if( pj_compare_datums( srcdefn, dstdefn ) )
643        return 0;
644
645    src_a = srcdefn->a_orig;
646    src_es = srcdefn->es_orig;
647
648    dst_a = dstdefn->a_orig;
649    dst_es = dstdefn->es_orig;
650
651/* -------------------------------------------------------------------- */
652/*      Create a temporary Z array if one is not provided.              */
653/* -------------------------------------------------------------------- */
654    if( z == NULL )
655    {
656        int	bytes = sizeof(double) * point_count * point_offset;
657        z = (double *) pj_malloc(bytes);
658        memset( z, 0, bytes );
659        z_is_temp = TRUE;
660    }
661
662#define CHECK_RETURN {if( pj_errno != 0 && (pj_errno > 0 || transient_error[-pj_errno] == 0) ) { if( z_is_temp ) pj_dalloc(z); return pj_errno; }}
663
664/* -------------------------------------------------------------------- */
665/*	If this datum requires grid shifts, then apply it to geodetic   */
666/*      coordinates.                                                    */
667/* -------------------------------------------------------------------- */
668    if( srcdefn->datum_type == PJD_GRIDSHIFT )
669    {
670        pj_apply_gridshift( pj_param(srcdefn->params,"snadgrids").s, 0, 
671                            point_count, point_offset, x, y, z );
672        CHECK_RETURN;
673
674        src_a = SRS_WGS84_SEMIMAJOR;
675        src_es = SRS_WGS84_ESQUARED;
676    }
677
678    if( dstdefn->datum_type == PJD_GRIDSHIFT )
679    {
680        dst_a = SRS_WGS84_SEMIMAJOR;
681        dst_es = SRS_WGS84_ESQUARED;
682    }
683        
684/* ==================================================================== */
685/*      Do we need to go through geocentric coordinates?                */
686/* ==================================================================== */
687    if( src_es != dst_es || src_a != dst_a
688        || srcdefn->datum_type == PJD_3PARAM 
689        || srcdefn->datum_type == PJD_7PARAM
690        || dstdefn->datum_type == PJD_3PARAM 
691        || dstdefn->datum_type == PJD_7PARAM)
692    {
693/* -------------------------------------------------------------------- */
694/*      Convert to geocentric coordinates.                              */
695/* -------------------------------------------------------------------- */
696        pj_geodetic_to_geocentric( src_a, src_es,
697                                   point_count, point_offset, x, y, z );
698        CHECK_RETURN;
699
700/* -------------------------------------------------------------------- */
701/*      Convert between datums.                                         */
702/* -------------------------------------------------------------------- */
703        if( srcdefn->datum_type == PJD_3PARAM 
704            || srcdefn->datum_type == PJD_7PARAM )
705        {
706            pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
707            CHECK_RETURN;
708        }
709
710        if( dstdefn->datum_type == PJD_3PARAM 
711            || dstdefn->datum_type == PJD_7PARAM )
712        {
713            pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
714            CHECK_RETURN;
715        }
716
717/* -------------------------------------------------------------------- */
718/*      Convert back to geodetic coordinates.                           */
719/* -------------------------------------------------------------------- */
720        pj_geocentric_to_geodetic( dst_a, dst_es,
721                                   point_count, point_offset, x, y, z );
722        CHECK_RETURN;
723    }
724
725/* -------------------------------------------------------------------- */
726/*      Apply grid shift to destination if required.                    */
727/* -------------------------------------------------------------------- */
728    if( dstdefn->datum_type == PJD_GRIDSHIFT )
729    {
730        pj_apply_gridshift( pj_param(dstdefn->params,"snadgrids").s, 1,
731                            point_count, point_offset, x, y, z );
732        CHECK_RETURN;
733    }
734
735    if( z_is_temp )
736        pj_dalloc( z );
737
738    return 0;
739}
740