← Back to branch summary

~ubuntu-branches/ubuntu/wily/proj/wily

« back to all changes in this revision

Viewing changes to src/pj_transform.c

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine, Jerome Villeneuve Larouche, Francesco Paolo Lovergine
  • Date: 2013-11-25 15:11:25 UTC
  • mfrom: (1.2.7)
  • Revision ID: package-import@ubuntu.com-20131125151125-mvcw144wvgep1hev
Tags: 4.8.0-1
[ Jerome Villeneuve Larouche ]
* New upstream release
* Modified libproj-dev.install to remove nad_list.h and projects.h
* Modified proj-bin.install to remove nad2nad
* Modified proj-bin.manpages to remove nad2nad man
* Added symbols for libproj0

[ Francesco Paolo Lovergine ]
* Properly merged with current git master and sid version 4.7.0-2.
* Removed proj transitional package and obsolete conflicts/replaces against
  series 4.6 and older.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/pj_transform.c
1
1
/******************************************************************************
2
 
 * $Id: pj_transform.c 1504 2009-01-06 02:11:57Z warmerdam $
 
2
 * $Id: pj_transform.c 2000 2011-05-10 17:06:33Z warmerdam $
3
3
 *
4
4
 * Project:  PROJ.4
5
5
 * Purpose:  Perform overall coordinate system to coordinate system 
34
34
#include <math.h>
35
35
#include "geocent.h"
36
36
 
37
 
PJ_CVSID("$Id: pj_transform.c 1504 2009-01-06 02:11:57Z warmerdam $");
 
37
PJ_CVSID("$Id: pj_transform.c 2000 2011-05-10 17:06:33Z warmerdam $");
 
38
 
 
39
static int pj_adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
 
40
                           long point_count, int point_offset, 
 
41
                           double *x, double *y, double *z );
38
42
 
39
43
#ifndef SRS_WGS84_SEMIMAJOR
40
44
#define SRS_WGS84_SEMIMAJOR 6378137.0
63
67
** list or something, but while experimenting with it this should be fine. 
64
68
*/
65
69
 
66
 
static const int transient_error[45] = {
 
70
static const int transient_error[50] = {
67
71
    /*             0  1  2  3  4  5  6  7  8  9   */
68
72
    /* 0 to 9 */   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
69
73
    /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,  
70
74
    /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
71
 
    /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
72
 
    /* 40 to 44 */ 0, 0, 0, 0, 0 };
 
75
    /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 
76
    /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 };
73
77
 
74
78
/************************************************************************/
75
79
/*                            pj_transform()                            */
85
89
 
86
90
{
87
91
    long      i;
88
 
    int       need_datum_shift;
 
92
    int       err;
89
93
 
90
 
    pj_errno = 0;
 
94
    srcdefn->ctx->last_errno = 0;
 
95
    dstdefn->ctx->last_errno = 0;
91
96
 
92
97
    if( point_offset == 0 )
93
98
        point_offset = 1;
94
99
 
95
100
/* -------------------------------------------------------------------- */
 
101
/*      Transform unusual input coordinate axis orientation to          */
 
102
/*      standard form if needed.                                        */
 
103
/* -------------------------------------------------------------------- */
 
104
    if( strcmp(srcdefn->axis,"enu") != 0 )
 
105
    {
 
106
        int err;
 
107
 
 
108
        err = pj_adjust_axis( srcdefn->ctx, srcdefn->axis, 
 
109
                              0, point_count, point_offset, x, y, z );
 
110
        if( err != 0 )
 
111
            return err;
 
112
    }
 
113
 
 
114
/* -------------------------------------------------------------------- */
 
115
/*      Transform Z to meters if it isn't already.                      */
 
116
/* -------------------------------------------------------------------- */
 
117
    if( srcdefn->vto_meter != 1.0 && z != NULL )
 
118
    {
 
119
        for( i = 0; i < point_count; i++ )
 
120
            z[point_offset*i] *= srcdefn->vto_meter;
 
121
    }
 
122
 
 
123
/* -------------------------------------------------------------------- */
96
124
/*      Transform geocentric source coordinates to lat/long.            */
97
125
/* -------------------------------------------------------------------- */
98
126
    if( srcdefn->is_geocent )
99
127
    {
100
128
        if( z == NULL )
101
129
        {
102
 
            pj_errno = PJD_ERR_GEOCENTRIC;
 
130
            pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
103
131
            return PJD_ERR_GEOCENTRIC;
104
132
        }
105
133
 
115
143
            }
116
144
        }
117
145
 
118
 
        if( pj_geocentric_to_geodetic( srcdefn->a_orig, srcdefn->es_orig,
119
 
                                       point_count, point_offset, 
120
 
                                       x, y, z ) != 0) 
121
 
            return pj_errno;
 
146
        err = pj_geocentric_to_geodetic( srcdefn->a_orig, srcdefn->es_orig,
 
147
                                         point_count, point_offset, 
 
148
                                         x, y, z );
 
149
        if( err != 0 )
 
150
            return err;
122
151
    }
123
152
 
124
153
/* -------------------------------------------------------------------- */
129
158
    {
130
159
        if( srcdefn->inv == NULL )
131
160
        {
132
 
            pj_errno = -17; /* this isn't correct, we need a no inverse err */
133
 
            if( getenv( "PROJ_DEBUG" ) != NULL )
134
 
            {
135
 
                fprintf( stderr, 
136
 
                       "pj_transform(): source projection not invertable\n" );
137
 
            }
138
 
            return pj_errno;
 
161
            pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 );
 
162
            pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR, 
 
163
                    "pj_transform(): source projection not invertable" );
 
164
            return -17;
139
165
        }
140
166
 
141
167
        for( i = 0; i < point_count; i++ )
150
176
                continue;
151
177
 
152
178
            geodetic_loc = pj_inv( projected_loc, srcdefn );
153
 
            if( pj_errno != 0 )
 
179
            if( srcdefn->ctx->last_errno != 0 )
154
180
            {
155
 
                if( (pj_errno != 33 /*EDOM*/ && pj_errno != 34 /*ERANGE*/ )
156
 
                    && (pj_errno > 0 || pj_errno < -44 || point_count == 1
157
 
                        || transient_error[-pj_errno] == 0 ) )
158
 
                    return pj_errno;
 
181
                if( (srcdefn->ctx->last_errno != 33 /*EDOM*/ 
 
182
                     && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
 
183
                    && (srcdefn->ctx->last_errno > 0 
 
184
                        || srcdefn->ctx->last_errno < -44 || point_count == 1
 
185
                        || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
 
186
                    return srcdefn->ctx->last_errno;
159
187
                else
160
188
                {
161
189
                    geodetic_loc.u = HUGE_VAL;
181
209
    }
182
210
 
183
211
/* -------------------------------------------------------------------- */
 
212
/*      Do we need to translate from geoid to ellipsoidal vertical      */
 
213
/*      datum?                                                          */
 
214
/* -------------------------------------------------------------------- */
 
215
    if( srcdefn->has_geoid_vgrids )
 
216
    {
 
217
        if( pj_apply_vgridshift( srcdefn, "sgeoidgrids", 
 
218
                                 &(srcdefn->vgridlist_geoid), 
 
219
                                 &(srcdefn->vgridlist_geoid_count),
 
220
                                 0, point_count, point_offset, x, y, z ) != 0 )
 
221
            return pj_ctx_get_errno(srcdefn->ctx);
 
222
    }
 
223
        
 
224
/* -------------------------------------------------------------------- */
184
225
/*      Convert datums if needed, and possible.                         */
185
226
/* -------------------------------------------------------------------- */
186
227
    if( pj_datum_transform( srcdefn, dstdefn, point_count, point_offset, 
187
228
                            x, y, z ) != 0 )
188
 
        return pj_errno;
 
229
    {
 
230
        if( srcdefn->ctx->last_errno != 0 )
 
231
            return srcdefn->ctx->last_errno;
 
232
        else
 
233
            return dstdefn->ctx->last_errno;
 
234
    }
189
235
 
190
236
/* -------------------------------------------------------------------- */
 
237
/*      Do we need to translate from geoid to ellipsoidal vertical      */
 
238
/*      datum?                                                          */
 
239
/* -------------------------------------------------------------------- */
 
240
    if( dstdefn->has_geoid_vgrids )
 
241
    {
 
242
        if( pj_apply_vgridshift( dstdefn, "sgeoidgrids", 
 
243
                                 &(dstdefn->vgridlist_geoid), 
 
244
                                 &(dstdefn->vgridlist_geoid_count),
 
245
                                 1, point_count, point_offset, x, y, z ) != 0 )
 
246
            return dstdefn->ctx->last_errno;
 
247
    }
 
248
        
 
249
/* -------------------------------------------------------------------- */
191
250
/*      But if they are staying lat long, adjust for the prime          */
192
251
/*      meridian if there is one in effect.                             */
193
252
/* -------------------------------------------------------------------- */
208
267
    {
209
268
        if( z == NULL )
210
269
        {
211
 
            pj_errno = PJD_ERR_GEOCENTRIC;
 
270
            pj_ctx_set_errno( dstdefn->ctx, PJD_ERR_GEOCENTRIC );
212
271
            return PJD_ERR_GEOCENTRIC;
213
272
        }
214
273
 
246
305
                continue;
247
306
 
248
307
            projected_loc = pj_fwd( geodetic_loc, dstdefn );
249
 
            if( pj_errno != 0 )
 
308
            if( dstdefn->ctx->last_errno != 0 )
250
309
            {
251
 
                if( (pj_errno != 33 /*EDOM*/ && pj_errno != 34 /*ERANGE*/ )
252
 
                    && (pj_errno > 0 || pj_errno < -44 || point_count == 1
253
 
                        || transient_error[-pj_errno] == 0 ) )
254
 
                    return pj_errno;
 
310
                if( (dstdefn->ctx->last_errno != 33 /*EDOM*/ 
 
311
                     && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
 
312
                    && (dstdefn->ctx->last_errno > 0 
 
313
                        || dstdefn->ctx->last_errno < -44 || point_count == 1
 
314
                        || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
 
315
                    return dstdefn->ctx->last_errno;
255
316
                else
256
317
                {
257
318
                    projected_loc.u = HUGE_VAL;
268
329
/*      If a wrapping center other than 0 is provided, rewrap around    */
269
330
/*      the suggested center (for latlong coordinate systems only).     */
270
331
/* -------------------------------------------------------------------- */
271
 
    else if( dstdefn->is_latlong && dstdefn->long_wrap_center != 0 )
 
332
    else if( dstdefn->is_latlong && dstdefn->is_long_wrap_set )
272
333
    {
273
334
        for( i = 0; i < point_count; i++ )
274
335
        {
275
336
            if( x[point_offset*i] == HUGE_VAL )
276
337
                continue;
277
338
 
278
 
            while( x[point_offset*i] < dstdefn->long_wrap_center - HALFPI )
279
 
                x[point_offset*i] += PI;
280
 
            while( x[point_offset*i] > dstdefn->long_wrap_center + HALFPI )
281
 
                x[point_offset*i] -= PI;
 
339
            while( x[point_offset*i] < dstdefn->long_wrap_center - PI )
 
340
                x[point_offset*i] += TWOPI;
 
341
            while( x[point_offset*i] > dstdefn->long_wrap_center + PI )
 
342
                x[point_offset*i] -= TWOPI;
282
343
        }
283
344
    }
284
345
 
 
346
/* -------------------------------------------------------------------- */
 
347
/*      Transform Z from meters if needed.                              */
 
348
/* -------------------------------------------------------------------- */
 
349
    if( dstdefn->vto_meter != 1.0 && z != NULL )
 
350
    {
 
351
        for( i = 0; i < point_count; i++ )
 
352
            z[point_offset*i] *= dstdefn->vfr_meter;
 
353
    }
 
354
 
 
355
/* -------------------------------------------------------------------- */
 
356
/*      Transform normalized axes into unusual output coordinate axis   */
 
357
/*      orientation if needed.                                          */
 
358
/* -------------------------------------------------------------------- */
 
359
    if( strcmp(dstdefn->axis,"enu") != 0 )
 
360
    {
 
361
        int err;
 
362
 
 
363
        err = pj_adjust_axis( dstdefn->ctx, dstdefn->axis, 
 
364
                              1, point_count, point_offset, x, y, z );
 
365
        if( err != 0 )
 
366
            return err;
 
367
    }
 
368
 
285
369
    return 0;
286
370
}
287
371
 
297
381
    double b;
298
382
    int    i;
299
383
    GeocentricInfo gi;
300
 
 
301
 
    pj_errno = 0;
 
384
    int ret_errno = 0;
302
385
 
303
386
    if( es == 0.0 )
304
387
        b = a;
307
390
 
308
391
    if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
309
392
    {
310
 
        pj_errno = PJD_ERR_GEOCENTRIC;
311
 
        return pj_errno;
 
393
        return PJD_ERR_GEOCENTRIC;
312
394
    }
313
395
 
314
396
    for( i = 0; i < point_count; i++ )
321
403
        if( pj_Convert_Geodetic_To_Geocentric( &gi, y[io], x[io], z[io], 
322
404
                                               x+io, y+io, z+io ) != 0 )
323
405
        {
324
 
            pj_errno = -14;
 
406
            ret_errno = -14;
325
407
            x[io] = y[io] = HUGE_VAL;
326
408
            /* but keep processing points! */
327
409
        }
328
410
    }
329
411
 
330
 
    return pj_errno;
 
412
    return ret_errno;
331
413
}
332
414
 
333
415
/************************************************************************/
350
432
 
351
433
    if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
352
434
    {
353
 
        pj_errno = PJD_ERR_GEOCENTRIC;
354
 
        return pj_errno;
 
435
        return PJD_ERR_GEOCENTRIC;
355
436
    }
356
437
 
357
438
    for( i = 0; i < point_count; i++ )
407
488
    }
408
489
    else if( srcdefn->datum_type == PJD_GRIDSHIFT )
409
490
    {
410
 
        return strcmp( pj_param(srcdefn->params,"snadgrids").s,
411
 
                       pj_param(dstdefn->params,"snadgrids").s ) == 0;
 
491
        return strcmp( pj_param(srcdefn->ctx, srcdefn->params,"snadgrids").s,
 
492
                       pj_param(dstdefn->ctx, dstdefn->params,"snadgrids").s ) == 0;
412
493
    }
413
494
    else
414
495
        return 1;
425
506
{
426
507
    int       i;
427
508
 
428
 
    pj_errno = 0;
429
 
 
430
509
    if( defn->datum_type == PJD_3PARAM )
431
510
    {
432
511
        for( i = 0; i < point_count; i++ )
475
554
{
476
555
    int       i;
477
556
 
478
 
    pj_errno = 0;
479
 
 
480
557
    if( defn->datum_type == PJD_3PARAM )
481
558
    {
482
559
        for( i = 0; i < point_count; i++ )
530
607
    double      src_a, src_es, dst_a, dst_es;
531
608
    int         z_is_temp = FALSE;
532
609
 
533
 
    pj_errno = 0;
534
 
 
535
610
/* -------------------------------------------------------------------- */
536
611
/*      We cannot do any meaningful datum transformation if either      */
537
612
/*      the source or destination are of an unknown datum type          */
565
640
        z_is_temp = TRUE;
566
641
    }
567
642
 
568
 
#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; }}
 
643
#define CHECK_RETURN(defn) {if( defn->ctx->last_errno != 0 && (defn->ctx->last_errno > 0 || transient_error[-defn->ctx->last_errno] == 0) ) { if( z_is_temp ) pj_dalloc(z); return defn->ctx->last_errno; }}
569
644
 
570
645
/* -------------------------------------------------------------------- */
571
646
/*      If this datum requires grid shifts, then apply it to geodetic   */
573
648
/* -------------------------------------------------------------------- */
574
649
    if( srcdefn->datum_type == PJD_GRIDSHIFT )
575
650
    {
576
 
        pj_apply_gridshift( pj_param(srcdefn->params,"snadgrids").s, 0, 
577
 
                            point_count, point_offset, x, y, z );
578
 
        CHECK_RETURN;
 
651
        pj_apply_gridshift_2( srcdefn, 0, point_count, point_offset, x, y, z );
 
652
        CHECK_RETURN(srcdefn);
579
653
 
580
654
        src_a = SRS_WGS84_SEMIMAJOR;
581
655
        src_es = SRS_WGS84_ESQUARED;
586
660
        dst_a = SRS_WGS84_SEMIMAJOR;
587
661
        dst_es = SRS_WGS84_ESQUARED;
588
662
    }
589
 
        
 
663
 
590
664
/* ==================================================================== */
591
665
/*      Do we need to go through geocentric coordinates?                */
592
666
/* ==================================================================== */
599
673
/* -------------------------------------------------------------------- */
600
674
/*      Convert to geocentric coordinates.                              */
601
675
/* -------------------------------------------------------------------- */
602
 
        pj_geodetic_to_geocentric( src_a, src_es,
603
 
                                   point_count, point_offset, x, y, z );
604
 
        CHECK_RETURN;
 
676
        srcdefn->ctx->last_errno = 
 
677
            pj_geodetic_to_geocentric( src_a, src_es,
 
678
                                       point_count, point_offset, x, y, z );
 
679
        CHECK_RETURN(srcdefn);
605
680
 
606
681
/* -------------------------------------------------------------------- */
607
682
/*      Convert between datums.                                         */
610
685
            || srcdefn->datum_type == PJD_7PARAM )
611
686
        {
612
687
            pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
613
 
            CHECK_RETURN;
 
688
            CHECK_RETURN(srcdefn);
614
689
        }
615
690
 
616
691
        if( dstdefn->datum_type == PJD_3PARAM 
617
692
            || dstdefn->datum_type == PJD_7PARAM )
618
693
        {
619
694
            pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
620
 
            CHECK_RETURN;
 
695
            CHECK_RETURN(dstdefn);
621
696
        }
622
697
 
623
698
/* -------------------------------------------------------------------- */
624
699
/*      Convert back to geodetic coordinates.                           */
625
700
/* -------------------------------------------------------------------- */
626
 
        pj_geocentric_to_geodetic( dst_a, dst_es,
627
 
                                   point_count, point_offset, x, y, z );
628
 
        CHECK_RETURN;
 
701
        dstdefn->ctx->last_errno = 
 
702
            pj_geocentric_to_geodetic( dst_a, dst_es,
 
703
                                       point_count, point_offset, x, y, z );
 
704
        CHECK_RETURN(dstdefn);
629
705
    }
630
706
 
631
707
/* -------------------------------------------------------------------- */
633
709
/* -------------------------------------------------------------------- */
634
710
    if( dstdefn->datum_type == PJD_GRIDSHIFT )
635
711
    {
636
 
        pj_apply_gridshift( pj_param(dstdefn->params,"snadgrids").s, 1,
637
 
                            point_count, point_offset, x, y, z );
638
 
        CHECK_RETURN;
 
712
        pj_apply_gridshift_2( dstdefn, 1, point_count, point_offset, x, y, z );
 
713
        CHECK_RETURN(dstdefn);
639
714
    }
640
715
 
641
716
    if( z_is_temp )
644
719
    return 0;
645
720
}
646
721
 
 
722
/************************************************************************/
 
723
/*                           pj_adjust_axis()                           */
 
724
/*                                                                      */
 
725
/*      Normalize or de-normalized the x/y/z axes.  The normal form     */
 
726
/*      is "enu" (easting, northing, up).                               */
 
727
/************************************************************************/
 
728
static int pj_adjust_axis( projCtx ctx, 
 
729
                           const char *axis, int denormalize_flag,
 
730
                           long point_count, int point_offset, 
 
731
                           double *x, double *y, double *z )
 
732
 
 
733
{
 
734
    double x_in, y_in, z_in = 0.0;
 
735
    int i, i_axis;
 
736
 
 
737
    if( !denormalize_flag )
 
738
    {
 
739
        for( i = 0; i < point_count; i++ )
 
740
        {
 
741
            x_in = x[point_offset*i];
 
742
            y_in = y[point_offset*i];
 
743
            if( z )
 
744
                z_in = z[point_offset*i];
 
745
     
 
746
            for( i_axis = 0; i_axis < 3; i_axis++ )
 
747
            {
 
748
                double value;
 
749
 
 
750
                if( i_axis == 0 )
 
751
                    value = x_in;
 
752
                else if( i_axis == 1 )
 
753
                    value = y_in;
 
754
                else
 
755
                    value = z_in;
 
756
                
 
757
                switch( axis[i_axis] )
 
758
                {
 
759
                  case 'e':
 
760
                    x[point_offset*i] = value; break;
 
761
                  case 'w':
 
762
                    x[point_offset*i] = -value; break;
 
763
                  case 'n':
 
764
                    y[point_offset*i] = value; break;
 
765
                  case 's':
 
766
                    y[point_offset*i] = -value; break;
 
767
                  case 'u':
 
768
                    if( z ) z[point_offset*i] = value; break;
 
769
                  case 'd':
 
770
                    if( z ) z[point_offset*i] = -value; break;
 
771
                  default:
 
772
                    pj_ctx_set_errno( ctx, PJD_ERR_AXIS );
 
773
                    return PJD_ERR_AXIS;
 
774
                }
 
775
            } /* i_axis */
 
776
        } /* i (point) */
 
777
    }
 
778
 
 
779
    else /* denormalize */
 
780
    {
 
781
        for( i = 0; i < point_count; i++ )
 
782
        {
 
783
            x_in = x[point_offset*i];
 
784
            y_in = y[point_offset*i];
 
785
            if( z )
 
786
                z_in = z[point_offset*i];
 
787
     
 
788
            for( i_axis = 0; i_axis < 3; i_axis++ )
 
789
            {
 
790
                double *target;
 
791
 
 
792
                if( i_axis == 2 && z == NULL )
 
793
                    continue;
 
794
 
 
795
                if( i_axis == 0 )
 
796
                    target = x;
 
797
                else if( i_axis == 1 )
 
798
                    target = y;
 
799
                else
 
800
                    target = z;
 
801
                
 
802
                switch( axis[i_axis] )
 
803
                {
 
804
                  case 'e':
 
805
                    target[point_offset*i] = x_in; break;
 
806
                  case 'w':
 
807
                    target[point_offset*i] = -x_in; break;
 
808
                  case 'n':
 
809
                    target[point_offset*i] = y_in; break;
 
810
                  case 's':
 
811
                    target[point_offset*i] = -y_in; break;
 
812
                  case 'u':
 
813
                    target[point_offset*i] = z_in; break;
 
814
                  case 'd':
 
815
                    target[point_offset*i] = -z_in; break;
 
816
                  default:
 
817
                    pj_ctx_set_errno( ctx, PJD_ERR_AXIS );
 
818
                    return PJD_ERR_AXIS;
 
819
                }
 
820
            } /* i_axis */
 
821
        } /* i (point) */
 
822
    }
 
823
    
 
824
    return 0;
 
825
}
 
826