← Back to branch summary

~ubuntu-branches/ubuntu/breezy/proj/breezy

« back to all changes in this revision

Viewing changes to src/pj_transform.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter S Galbraith
  • Date: 2004-11-06 19:44:53 UTC
  • mto: This revision was merged to the branch mainline in revision 3.
  • Revision ID: james.westby@ubuntu.com-20041106194453-axnsmkh1zplal8mz
Tags: upstream-4.4.9
Import upstream version 4.4.9

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/pj_transform.c
1
1
/******************************************************************************
2
 
 * $Id: pj_transform.c,v 1.3 2001/04/04 21:13:21 warmerda Exp $
 
2
 * $Id: pj_transform.c,v 1.13 2004/10/25 15:34:36 fwarmerdam Exp $
3
3
 *
4
4
 * Project:  PROJ.4
5
5
 * Purpose:  Perform overall coordinate system to coordinate system 
6
6
 *           transformations (pj_transform() function) including reprojection
7
7
 *           and datum shifting.
8
 
 * Author:   Frank Warmerdam, warmerda@home.com
 
8
 * Author:   Frank Warmerdam, warmerdam@pobox.com
9
9
 *
10
10
 ******************************************************************************
11
11
 * Copyright (c) 2000, Frank Warmerdam
30
30
 ******************************************************************************
31
31
 *
32
32
 * $Log: pj_transform.c,v $
 
33
 * Revision 1.13  2004/10/25 15:34:36  fwarmerdam
 
34
 * make names of geodetic funcs from geotrans unique
 
35
 *
 
36
 * Revision 1.12  2004/05/03 19:45:23  warmerda
 
37
 * Altered so that raw ellpses are treated as a essentially having a
 
38
 * +towgs84=0,0,0 specification so ellpisoid shifts are applied.
 
39
 * Fixed so that prime meridian shifts are applied if the coordinate system is
 
40
 * not lat/long (ie. if it is projected).  This fixes:
 
41
 * http://bugzilla.remotesensing.org/show_bug.cgi?id=510
 
42
 *
 
43
 * Revision 1.11  2004/01/24 09:37:19  warmerda
 
44
 * pj_transform() will no longer return an error code if any of the points are
 
45
 * transformable.  In this case the application is expected to check for
 
46
 * HUGE_VAL to identify failed points.
 
47
 * As part of the implementation, I added a list of pj_errno values that
 
48
 * are transient (ie per-point) rather than indicating global failure for the
 
49
 * coordinate system definition.  We use this in deciding which pj_fwd and
 
50
 * pj_inv error codes are really fatal and should be reported.
 
51
 *
 
52
 * Revision 1.10  2003/08/21 02:09:06  warmerda
 
53
 * added a bunch of HUGE_VAL checking
 
54
 *
 
55
 * Revision 1.9  2003/03/26 16:52:30  warmerda
 
56
 * added check that an inverse transformation func exists
 
57
 *
 
58
 * Revision 1.8  2002/12/14 20:35:43  warmerda
 
59
 * implement units support for geocentric coordinates
 
60
 *
 
61
 * Revision 1.7  2002/12/14 20:14:35  warmerda
 
62
 * added geocentric support
 
63
 *
 
64
 * Revision 1.6  2002/12/09 16:01:02  warmerda
 
65
 * added prime meridian support
 
66
 *
 
67
 * Revision 1.5  2002/12/01 19:25:26  warmerda
 
68
 * applied fix for 7 param shift in pj_geocentric_from_wgs84, see bug 194
 
69
 *
 
70
 * Revision 1.4  2002/02/15 14:30:36  warmerda
 
71
 * provide default Z array if none passed in in pj_datum_transform()
 
72
 *
33
73
 * Revision 1.3  2001/04/04 21:13:21  warmerda
34
74
 * do arcsecond/radian and ppm datum parm transformation in pj_set_datum()
35
75
 *
46
86
#include <math.h>
47
87
#include "geocent.h"
48
88
 
 
89
PJ_CVSID("$Id: pj_transform.c,v 1.13 2004/10/25 15:34:36 fwarmerdam Exp $");
 
90
 
49
91
#ifndef SRS_WGS84_SEMIMAJOR
50
92
#define SRS_WGS84_SEMIMAJOR 6378137.0
51
93
#endif
62
104
#define Rz_BF (defn->datum_params[5])
63
105
#define M_BF  (defn->datum_params[6])
64
106
 
 
107
/* 
 
108
** This table is intended to indicate for any given error code in 
 
109
** the range 0 to -44, whether that error will occur for all locations (ie.
 
110
** it is a problem with the coordinate system as a whole) in which case the
 
111
** value would be 0, or if the problem is with the point being transformed
 
112
** in which case the value is 1. 
 
113
**
 
114
** At some point we might want to move this array in with the error message
 
115
** list or something, but while experimenting with it this should be fine. 
 
116
*/
 
117
 
 
118
static int transient_error[45] = {
 
119
    /*             0  1  2  3  4  5  6  7  8  9   */
 
120
    /* 0 to 9 */   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 
121
    /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,  
 
122
    /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
 
123
    /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
 
124
    /* 40 to 44 */ 0, 0, 0, 0, 0 };
 
125
 
65
126
/************************************************************************/
66
127
/*                            pj_transform()                            */
67
128
/*                                                                      */
84
145
        point_offset = 1;
85
146
 
86
147
/* -------------------------------------------------------------------- */
 
148
/*      Transform geocentric source coordinates to lat/long.            */
 
149
/* -------------------------------------------------------------------- */
 
150
    if( srcdefn->is_geocent )
 
151
    {
 
152
        if( z == NULL )
 
153
        {
 
154
            pj_errno = PJD_ERR_GEOCENTRIC;
 
155
            return PJD_ERR_GEOCENTRIC;
 
156
        }
 
157
 
 
158
        if( srcdefn->to_meter != 1.0 )
 
159
        {
 
160
            for( i = 0; i < point_count; i++ )
 
161
            {
 
162
                x[point_offset*i] *= srcdefn->to_meter;
 
163
                y[point_offset*i] *= srcdefn->to_meter;
 
164
            }
 
165
        }
 
166
 
 
167
        if( pj_geocentric_to_geodetic( srcdefn->a, srcdefn->es,
 
168
                                       point_count, point_offset, 
 
169
                                       x, y, z ) != 0) 
 
170
            return pj_errno;
 
171
    }
 
172
 
 
173
/* -------------------------------------------------------------------- */
87
174
/*      Transform source points to lat/long, if they aren't             */
88
175
/*      already.                                                        */
89
176
/* -------------------------------------------------------------------- */
90
 
    if( !srcdefn->is_latlong )
 
177
    else if( !srcdefn->is_latlong )
91
178
    {
 
179
        if( srcdefn->inv == NULL )
 
180
        {
 
181
            pj_errno = -17; /* this isn't correct, we need a no inverse err */
 
182
            if( getenv( "PROJ_DEBUG" ) != NULL )
 
183
            {
 
184
                fprintf( stderr, 
 
185
                       "pj_transform(): source projection not invertable\n" );
 
186
            }
 
187
            return pj_errno;
 
188
        }
 
189
 
92
190
        for( i = 0; i < point_count; i++ )
93
191
        {
94
192
            XY         projected_loc;
97
195
            projected_loc.u = x[point_offset*i];
98
196
            projected_loc.v = y[point_offset*i];
99
197
 
 
198
            if( projected_loc.u == HUGE_VAL )
 
199
                continue;
 
200
 
100
201
            geodetic_loc = pj_inv( projected_loc, srcdefn );
101
202
            if( pj_errno != 0 )
102
 
                return pj_errno;
 
203
            {
 
204
                if( pj_errno > 0 || pj_errno < -44 || point_count == 1
 
205
                    || transient_error[-pj_errno] == 0 )
 
206
                    return pj_errno;
 
207
                else
 
208
                {
 
209
                    geodetic_loc.u = HUGE_VAL;
 
210
                    geodetic_loc.v = HUGE_VAL;
 
211
                }
 
212
            }
103
213
 
104
214
            x[point_offset*i] = geodetic_loc.u;
105
215
            y[point_offset*i] = geodetic_loc.v;
106
216
        }
107
217
    }
 
218
/* -------------------------------------------------------------------- */
 
219
/*      But if they are already lat long, adjust for the prime          */
 
220
/*      meridian if there is one in effect.                             */
 
221
/* -------------------------------------------------------------------- */
 
222
    if( srcdefn->from_greenwich != 0.0 )
 
223
    {
 
224
        for( i = 0; i < point_count; i++ )
 
225
        {
 
226
            if( x[point_offset*i] != HUGE_VAL )
 
227
                x[point_offset*i] += srcdefn->from_greenwich;
 
228
        }
 
229
    }
108
230
 
109
231
/* -------------------------------------------------------------------- */
110
232
/*      Convert datums if needed, and possible.                         */
114
236
        return pj_errno;
115
237
 
116
238
/* -------------------------------------------------------------------- */
 
239
/*      But if they are staying lat long, adjust for the prime          */
 
240
/*      meridian if there is one in effect.                             */
 
241
/* -------------------------------------------------------------------- */
 
242
    if( dstdefn->from_greenwich != 0.0 )
 
243
    {
 
244
        for( i = 0; i < point_count; i++ )
 
245
        {
 
246
            if( x[point_offset*i] != HUGE_VAL )
 
247
                x[point_offset*i] -= dstdefn->from_greenwich;
 
248
        }
 
249
    }
 
250
 
 
251
 
 
252
/* -------------------------------------------------------------------- */
 
253
/*      Transform destination latlong to geocentric if required.        */
 
254
/* -------------------------------------------------------------------- */
 
255
    if( dstdefn->is_geocent )
 
256
    {
 
257
        if( z == NULL )
 
258
        {
 
259
            pj_errno = PJD_ERR_GEOCENTRIC;
 
260
            return PJD_ERR_GEOCENTRIC;
 
261
        }
 
262
 
 
263
        pj_geodetic_to_geocentric( dstdefn->a, dstdefn->es,
 
264
                                   point_count, point_offset, x, y, z );
 
265
 
 
266
        if( dstdefn->fr_meter != 1.0 )
 
267
        {
 
268
            for( i = 0; i < point_count; i++ )
 
269
            {
 
270
                if( x[point_offset*i] != HUGE_VAL )
 
271
                {
 
272
                    x[point_offset*i] *= dstdefn->fr_meter;
 
273
                    y[point_offset*i] *= dstdefn->fr_meter;
 
274
                }
 
275
            }
 
276
        }
 
277
    }
 
278
 
 
279
/* -------------------------------------------------------------------- */
117
280
/*      Transform destination points to projection coordinates, if      */
118
281
/*      desired.                                                        */
119
282
/* -------------------------------------------------------------------- */
120
 
    if( !dstdefn->is_latlong )
 
283
    else if( !dstdefn->is_latlong )
121
284
    {
122
285
        for( i = 0; i < point_count; i++ )
123
286
        {
127
290
            geodetic_loc.u = x[point_offset*i];
128
291
            geodetic_loc.v = y[point_offset*i];
129
292
 
 
293
            if( geodetic_loc.u == HUGE_VAL )
 
294
                continue;
 
295
 
130
296
            projected_loc = pj_fwd( geodetic_loc, dstdefn );
131
297
            if( pj_errno != 0 )
132
 
                return pj_errno;
 
298
            {
 
299
                if( pj_errno > 0 || pj_errno < -44 || point_count == 1
 
300
                    || transient_error[-pj_errno] == 0 )
 
301
                    return pj_errno;
 
302
                else
 
303
                {
 
304
                    projected_loc.u = HUGE_VAL;
 
305
                    projected_loc.v = HUGE_VAL;
 
306
                }
 
307
            }
133
308
 
134
309
            x[point_offset*i] = projected_loc.u;
135
310
            y[point_offset*i] = projected_loc.v;
156
331
    else
157
332
        b = a * sqrt(1-es);
158
333
 
159
 
    if( Set_Geocentric_Parameters( a, b ) != 0 )
 
334
    if( pj_Set_Geocentric_Parameters( a, b ) != 0 )
160
335
    {
161
336
        pj_errno = PJD_ERR_GEOCENTRIC;
162
337
        return pj_errno;
166
341
    {
167
342
        long io = i * point_offset;
168
343
 
169
 
        if( Convert_Geodetic_To_Geocentric( y[io], x[io], z[io], 
 
344
        if( pj_Convert_Geodetic_To_Geocentric( y[io], x[io], z[io], 
170
345
                                            x+io, y+io, z+io ) != 0 )
171
346
        {
172
347
            pj_errno = PJD_ERR_GEOCENTRIC;
194
369
    else
195
370
        b = a * sqrt(1-es);
196
371
 
197
 
    if( Set_Geocentric_Parameters( a, b ) != 0 )
 
372
    if( pj_Set_Geocentric_Parameters( a, b ) != 0 )
198
373
    {
199
374
        pj_errno = PJD_ERR_GEOCENTRIC;
200
375
        return pj_errno;
204
379
    {
205
380
        long io = i * point_offset;
206
381
 
207
 
        Convert_Geocentric_To_Geodetic( x[io], y[io], z[io], 
 
382
        if( x[io] == HUGE_VAL )
 
383
            continue;
 
384
 
 
385
        pj_Convert_Geocentric_To_Geodetic( x[io], y[io], z[io], 
208
386
                                        y+io, x+io, z+io );
209
387
    }
210
388
 
276
454
        {
277
455
            long io = i * point_offset;
278
456
            
 
457
            if( x[io] == HUGE_VAL )
 
458
                continue;
 
459
 
279
460
            x[io] = x[io] + defn->datum_params[0];
280
461
            y[io] = y[io] + defn->datum_params[1];
281
462
            z[io] = z[io] + defn->datum_params[2];
288
469
            long io = i * point_offset;
289
470
            double x_out, y_out, z_out;
290
471
 
 
472
            if( x[io] == HUGE_VAL )
 
473
                continue;
 
474
 
291
475
            x_out = M_BF*(       x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
292
476
            y_out = M_BF*( Rz_BF*x[io] +       y[io] - Rx_BF*z[io]) + Dy_BF;
293
477
            z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] +       z[io]) + Dz_BF;
319
503
        for( i = 0; i < point_count; i++ )
320
504
        {
321
505
            long io = i * point_offset;
 
506
 
 
507
            if( x[io] == HUGE_VAL )
 
508
                continue;
322
509
            
323
510
            x[io] = x[io] - defn->datum_params[0];
324
511
            y[io] = y[io] - defn->datum_params[1];
330
517
        for( i = 0; i < point_count; i++ )
331
518
        {
332
519
            long io = i * point_offset;
333
 
            double x_out, y_out, z_out;
334
 
 
335
 
            x_out = M_BF*(       x[io] + Rz_BF*y[io] - Ry_BF*z[io]) - Dx_BF;
336
 
            y_out = M_BF*(-Rz_BF*x[io] +       y[io] + Rx_BF*z[io]) - Dy_BF;
337
 
            z_out = M_BF*( Ry_BF*x[io] - Rx_BF*y[io] +       z[io]) - Dz_BF;
338
 
 
339
 
            x[io] = x_out;
340
 
            y[io] = y_out;
341
 
            z[io] = z_out;
 
520
            double x_tmp, y_tmp, z_tmp;
 
521
 
 
522
            if( x[io] == HUGE_VAL )
 
523
                continue;
 
524
 
 
525
            x_tmp = (x[io] - Dx_BF) / M_BF;
 
526
            y_tmp = (y[io] - Dy_BF) / M_BF;
 
527
            z_tmp = (z[io] - Dz_BF) / M_BF;
 
528
 
 
529
            x[io] =        x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp;
 
530
            y[io] = -Rz_BF*x_tmp +       y_tmp + Rx_BF*z_tmp;
 
531
            z[io] =  Ry_BF*x_tmp - Rx_BF*y_tmp +       z_tmp;
342
532
        }
343
533
    }
344
534
 
355
545
 
356
546
{
357
547
    double      src_a, src_es, dst_a, dst_es;
 
548
    int         z_is_temp = FALSE;
358
549
 
359
550
    pj_errno = 0;
360
551
 
371
562
    dst_es = dstdefn->es;
372
563
 
373
564
/* -------------------------------------------------------------------- */
 
565
/*      Create a temporary Z array if one is not provided.              */
 
566
/* -------------------------------------------------------------------- */
 
567
    if( z == NULL )
 
568
    {
 
569
        int     bytes = sizeof(double) * point_count * point_offset;
 
570
        z = (double *) pj_malloc(bytes);
 
571
        memset( z, 0, bytes );
 
572
        z_is_temp = TRUE;
 
573
    }
 
574
 
 
575
#define CHECK_RETURN {if( pj_errno != 0 ) { if( z_is_temp ) pj_dalloc(z); return pj_errno; }}
 
576
 
 
577
/* -------------------------------------------------------------------- */
374
578
/*      If this datum requires grid shifts, then apply it to geodetic   */
375
579
/*      coordinates.                                                    */
376
580
/* -------------------------------------------------------------------- */
378
582
    {
379
583
        pj_apply_gridshift( pj_param(srcdefn->params,"snadgrids").s, 0, 
380
584
                            point_count, point_offset, x, y, z );
381
 
 
382
 
        if( pj_errno != 0 )
383
 
            return pj_errno;
 
585
        CHECK_RETURN;
384
586
 
385
587
        src_a = SRS_WGS84_SEMIMAJOR;
386
588
        src_es = 0.006694379990;
395
597
/* ==================================================================== */
396
598
/*      Do we need to go through geocentric coordinates?                */
397
599
/* ==================================================================== */
398
 
    if( srcdefn->datum_type == PJD_3PARAM 
 
600
    if( src_es != dst_es || src_a != dst_a
 
601
        || srcdefn->datum_type == PJD_3PARAM 
399
602
        || srcdefn->datum_type == PJD_7PARAM
400
603
        || dstdefn->datum_type == PJD_3PARAM 
401
604
        || dstdefn->datum_type == PJD_7PARAM)
405
608
/* -------------------------------------------------------------------- */
406
609
        pj_geodetic_to_geocentric( src_a, src_es,
407
610
                                   point_count, point_offset, x, y, z );
408
 
 
409
 
        if( pj_errno )
410
 
            return pj_errno;
 
611
        CHECK_RETURN;
411
612
 
412
613
/* -------------------------------------------------------------------- */
413
614
/*      Convert between datums.                                         */
414
615
/* -------------------------------------------------------------------- */
415
 
        if( srcdefn->datum_type != PJD_UNKNOWN
416
 
            && dstdefn->datum_type != PJD_UNKNOWN )
 
616
        if( srcdefn->datum_type == PJD_3PARAM 
 
617
            || srcdefn->datum_type == PJD_7PARAM )
417
618
        {
418
619
            pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
419
 
            if( pj_errno != 0 )
420
 
                return pj_errno;
421
 
            
 
620
            CHECK_RETURN;
 
621
        }
 
622
 
 
623
        if( dstdefn->datum_type == PJD_3PARAM 
 
624
            || dstdefn->datum_type == PJD_7PARAM )
 
625
        {
422
626
            pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
423
 
            if( pj_errno != 0 )
424
 
                return pj_errno;
 
627
            CHECK_RETURN;
425
628
        }
426
629
 
427
630
/* -------------------------------------------------------------------- */
429
632
/* -------------------------------------------------------------------- */
430
633
        pj_geocentric_to_geodetic( dst_a, dst_es,
431
634
                                   point_count, point_offset, x, y, z );
432
 
        
433
 
        if( pj_errno )
434
 
            return pj_errno;
 
635
        CHECK_RETURN;
435
636
    }
436
637
 
437
638
/* -------------------------------------------------------------------- */
441
642
    {
442
643
        pj_apply_gridshift( pj_param(dstdefn->params,"snadgrids").s, 1,
443
644
                            point_count, point_offset, x, y, z );
444
 
 
445
 
        if( pj_errno != 0 )
446
 
            return pj_errno;
 
645
        CHECK_RETURN;
447
646
    }
448
647
 
 
648
    if( z_is_temp )
 
649
        pj_dalloc( z );
 
650
 
449
651
    return 0;
450
652
}
451
653