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 <projects.h>
31 #include <string.h>
32 #include <math.h>
33 #include "geocent.h"
34
35 static int pj_adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
36 long point_count, int point_offset,
37 double *x, double *y, double *z );
38
39 #ifndef SRS_WGS84_SEMIMAJOR
40 #define SRS_WGS84_SEMIMAJOR 6378137.0
41 #endif
42
43 #ifndef SRS_WGS84_ESQUARED
44 #define SRS_WGS84_ESQUARED 0.0066943799901413165
45 #endif
46
47 #define Dx_BF (defn->datum_params[0])
48 #define Dy_BF (defn->datum_params[1])
49 #define Dz_BF (defn->datum_params[2])
50 #define Rx_BF (defn->datum_params[3])
51 #define Ry_BF (defn->datum_params[4])
52 #define Rz_BF (defn->datum_params[5])
53 #define M_BF (defn->datum_params[6])
54
55 /*
56 ** This table is intended to indicate for any given error code in
57 ** the range 0 to -44, whether that error will occur for all locations (ie.
58 ** it is a problem with the coordinate system as a whole) in which case the
59 ** value would be 0, or if the problem is with the point being transformed
60 ** in which case the value is 1.
61 **
62 ** At some point we might want to move this array in with the error message
63 ** list or something, but while experimenting with it this should be fine.
64 */
65
66 static const int transient_error[50] = {
67 /* 0 1 2 3 4 5 6 7 8 9 */
68 /* 0 to 9 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
69 /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,
70 /* 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, 0, 0,
72 /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 };
73
74 /************************************************************************/
75 /* pj_transform() */
76 /* */
77 /* Currently this function doesn't recognise if two projections */
78 /* are identical (to short circuit reprojection) because it is */
79 /* difficult to compare PJ structures (since there are some */
80 /* projection specific components). */
81 /************************************************************************/
82
pj_transform(PJ * srcdefn,PJ * dstdefn,long point_count,int point_offset,double * x,double * y,double * z)83 int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
84 double *x, double *y, double *z )
85
86 {
87 long i;
88 int err;
89
90 srcdefn->ctx->last_errno = 0;
91 dstdefn->ctx->last_errno = 0;
92
93 if( point_offset == 0 )
94 point_offset = 1;
95
96 /* -------------------------------------------------------------------- */
97 /* Transform unusual input coordinate axis orientation to */
98 /* standard form if needed. */
99 /* -------------------------------------------------------------------- */
100 if( strcmp(srcdefn->axis,"enu") != 0 )
101 {
102 int err;
103
104 err = pj_adjust_axis( srcdefn->ctx, srcdefn->axis,
105 0, point_count, point_offset, x, y, z );
106 if( err != 0 )
107 return err;
108 }
109
110 /* -------------------------------------------------------------------- */
111 /* Transform Z to meters if it isn't already. */
112 /* -------------------------------------------------------------------- */
113 if( srcdefn->vto_meter != 1.0 && z != NULL )
114 {
115 for( i = 0; i < point_count; i++ )
116 z[point_offset*i] *= srcdefn->vto_meter;
117 }
118
119 /* -------------------------------------------------------------------- */
120 /* Transform geocentric source coordinates to lat/long. */
121 /* -------------------------------------------------------------------- */
122 if( srcdefn->is_geocent )
123 {
124 if( z == NULL )
125 {
126 pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
127 return PJD_ERR_GEOCENTRIC;
128 }
129
130 if( srcdefn->to_meter != 1.0 )
131 {
132 for( i = 0; i < point_count; i++ )
133 {
134 if( x[point_offset*i] != HUGE_VAL )
135 {
136 x[point_offset*i] *= srcdefn->to_meter;
137 y[point_offset*i] *= srcdefn->to_meter;
138 }
139 }
140 }
141
142 err = pj_geocentric_to_geodetic( srcdefn->a_orig, srcdefn->es_orig,
143 point_count, point_offset,
144 x, y, z );
145 if( err != 0 )
146 return err;
147 }
148
149 /* -------------------------------------------------------------------- */
150 /* Transform source points to lat/long, if they aren't */
151 /* already. */
152 /* -------------------------------------------------------------------- */
153 else if( !srcdefn->is_latlong )
154 {
155
156 /* Check first if projection is invertible. */
157 if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL))
158 {
159 pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 );
160 pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR,
161 "pj_transform(): source projection not invertable" );
162 return -17;
163 }
164
165 /* If invertible - First try inv3d if defined */
166 if (srcdefn->inv3d != NULL)
167 {
168 /* Three dimensions must be defined */
169 if ( z == NULL)
170 {
171 pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
172 return PJD_ERR_GEOCENTRIC;
173 }
174
175 for (i=0; i < point_count; i++)
176 {
177 XYZ projected_loc;
178 XYZ geodetic_loc;
179
180 projected_loc.u = x[point_offset*i];
181 projected_loc.v = y[point_offset*i];
182 projected_loc.w = z[point_offset*i];
183
184 if (projected_loc.u == HUGE_VAL)
185 continue;
186
187 geodetic_loc = pj_inv3d(projected_loc, srcdefn);
188 if( srcdefn->ctx->last_errno != 0 )
189 {
190 if( (srcdefn->ctx->last_errno != 33 /*EDOM*/
191 && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
192 && (srcdefn->ctx->last_errno > 0
193 || srcdefn->ctx->last_errno < -44 || point_count == 1
194 || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
195 return srcdefn->ctx->last_errno;
196 else
197 {
198 geodetic_loc.u = HUGE_VAL;
199 geodetic_loc.v = HUGE_VAL;
200 geodetic_loc.w = HUGE_VAL;
201 }
202 }
203
204 x[point_offset*i] = geodetic_loc.u;
205 y[point_offset*i] = geodetic_loc.v;
206 z[point_offset*i] = geodetic_loc.w;
207
208 }
209
210 }
211 else
212 {
213 /* Fallback to the original PROJ.4 API 2d inversion - inv */
214 for( i = 0; i < point_count; i++ )
215 {
216 XY projected_loc;
217 LP geodetic_loc;
218
219 projected_loc.u = x[point_offset*i];
220 projected_loc.v = y[point_offset*i];
221
222 if( projected_loc.u == HUGE_VAL )
223 continue;
224
225 geodetic_loc = pj_inv( projected_loc, srcdefn );
226 if( srcdefn->ctx->last_errno != 0 )
227 {
228 if( (srcdefn->ctx->last_errno != 33 /*EDOM*/
229 && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
230 && (srcdefn->ctx->last_errno > 0
231 || srcdefn->ctx->last_errno < -44 || point_count == 1
232 || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
233 return srcdefn->ctx->last_errno;
234 else
235 {
236 geodetic_loc.u = HUGE_VAL;
237 geodetic_loc.v = HUGE_VAL;
238 }
239 }
240
241 x[point_offset*i] = geodetic_loc.u;
242 y[point_offset*i] = geodetic_loc.v;
243 }
244 }
245 }
246 /* -------------------------------------------------------------------- */
247 /* But if they are already lat long, adjust for the prime */
248 /* meridian if there is one in effect. */
249 /* -------------------------------------------------------------------- */
250 if( srcdefn->from_greenwich != 0.0 )
251 {
252 for( i = 0; i < point_count; i++ )
253 {
254 if( x[point_offset*i] != HUGE_VAL )
255 x[point_offset*i] += srcdefn->from_greenwich;
256 }
257 }
258
259 /* -------------------------------------------------------------------- */
260 /* Do we need to translate from geoid to ellipsoidal vertical */
261 /* datum? */
262 /* -------------------------------------------------------------------- */
263 if( srcdefn->has_geoid_vgrids && z != NULL )
264 {
265 if( pj_apply_vgridshift( srcdefn, "sgeoidgrids",
266 &(srcdefn->vgridlist_geoid),
267 &(srcdefn->vgridlist_geoid_count),
268 0, point_count, point_offset, x, y, z ) != 0 )
269 return pj_ctx_get_errno(srcdefn->ctx);
270 }
271
272 /* -------------------------------------------------------------------- */
273 /* Convert datums if needed, and possible. */
274 /* -------------------------------------------------------------------- */
275 if( pj_datum_transform( srcdefn, dstdefn, point_count, point_offset,
276 x, y, z ) != 0 )
277 {
278 if( srcdefn->ctx->last_errno != 0 )
279 return srcdefn->ctx->last_errno;
280 else
281 return dstdefn->ctx->last_errno;
282 }
283
284 /* -------------------------------------------------------------------- */
285 /* Do we need to translate from geoid to ellipsoidal vertical */
286 /* datum? */
287 /* -------------------------------------------------------------------- */
288 if( dstdefn->has_geoid_vgrids && z != NULL )
289 {
290 if( pj_apply_vgridshift( dstdefn, "sgeoidgrids",
291 &(dstdefn->vgridlist_geoid),
292 &(dstdefn->vgridlist_geoid_count),
293 1, point_count, point_offset, x, y, z ) != 0 )
294 return dstdefn->ctx->last_errno;
295 }
296
297 /* -------------------------------------------------------------------- */
298 /* But if they are staying lat long, adjust for the prime */
299 /* meridian if there is one in effect. */
300 /* -------------------------------------------------------------------- */
301 if( dstdefn->from_greenwich != 0.0 )
302 {
303 for( i = 0; i < point_count; i++ )
304 {
305 if( x[point_offset*i] != HUGE_VAL )
306 x[point_offset*i] -= dstdefn->from_greenwich;
307 }
308 }
309
310
311 /* -------------------------------------------------------------------- */
312 /* Transform destination latlong to geocentric if required. */
313 /* -------------------------------------------------------------------- */
314 if( dstdefn->is_geocent )
315 {
316 if( z == NULL )
317 {
318 pj_ctx_set_errno( dstdefn->ctx, PJD_ERR_GEOCENTRIC );
319 return PJD_ERR_GEOCENTRIC;
320 }
321
322 pj_geodetic_to_geocentric( dstdefn->a_orig, dstdefn->es_orig,
323 point_count, point_offset, x, y, z );
324
325 if( dstdefn->fr_meter != 1.0 )
326 {
327 for( i = 0; i < point_count; i++ )
328 {
329 if( x[point_offset*i] != HUGE_VAL )
330 {
331 x[point_offset*i] *= dstdefn->fr_meter;
332 y[point_offset*i] *= dstdefn->fr_meter;
333 }
334 }
335 }
336 }
337
338 /* -------------------------------------------------------------------- */
339 /* Transform destination points to projection coordinates, if */
340 /* desired. */
341 /* -------------------------------------------------------------------- */
342 else if( !dstdefn->is_latlong )
343 {
344
345 if( dstdefn->fwd3d != NULL)
346 {
347 for( i = 0; i < point_count; i++ )
348 {
349 XYZ projected_loc;
350 LPZ geodetic_loc;
351
352 geodetic_loc.u = x[point_offset*i];
353 geodetic_loc.v = y[point_offset*i];
354 geodetic_loc.w = z[point_offset*i];
355
356 if (geodetic_loc.u == HUGE_VAL)
357 continue;
358
359 projected_loc = pj_fwd3d( geodetic_loc, dstdefn);
360 if( dstdefn->ctx->last_errno != 0 )
361 {
362 if( (dstdefn->ctx->last_errno != 33 /*EDOM*/
363 && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
364 && (dstdefn->ctx->last_errno > 0
365 || dstdefn->ctx->last_errno < -44 || point_count == 1
366 || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
367 return dstdefn->ctx->last_errno;
368 else
369 {
370 projected_loc.u = HUGE_VAL;
371 projected_loc.v = HUGE_VAL;
372 projected_loc.w = HUGE_VAL;
373 }
374 }
375
376 x[point_offset*i] = projected_loc.u;
377 y[point_offset*i] = projected_loc.v;
378 z[point_offset*i] = projected_loc.w;
379 }
380
381 }
382 else
383 {
384 for( i = 0; i < point_count; i++ )
385 {
386 XY projected_loc;
387 LP geodetic_loc;
388
389 geodetic_loc.u = x[point_offset*i];
390 geodetic_loc.v = y[point_offset*i];
391
392 if( geodetic_loc.u == HUGE_VAL )
393 continue;
394
395 projected_loc = pj_fwd( geodetic_loc, dstdefn );
396 if( dstdefn->ctx->last_errno != 0 )
397 {
398 if( (dstdefn->ctx->last_errno != 33 /*EDOM*/
399 && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
400 && (dstdefn->ctx->last_errno > 0
401 || dstdefn->ctx->last_errno < -44 || point_count == 1
402 || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
403 return dstdefn->ctx->last_errno;
404 else
405 {
406 projected_loc.u = HUGE_VAL;
407 projected_loc.v = HUGE_VAL;
408 }
409 }
410
411 x[point_offset*i] = projected_loc.u;
412 y[point_offset*i] = projected_loc.v;
413 }
414 }
415 }
416
417 /* -------------------------------------------------------------------- */
418 /* If a wrapping center other than 0 is provided, rewrap around */
419 /* the suggested center (for latlong coordinate systems only). */
420 /* -------------------------------------------------------------------- */
421 else if( dstdefn->is_latlong && dstdefn->is_long_wrap_set )
422 {
423 for( i = 0; i < point_count; i++ )
424 {
425 if( x[point_offset*i] == HUGE_VAL )
426 continue;
427
428 while( x[point_offset*i] < dstdefn->long_wrap_center - M_PI )
429 x[point_offset*i] += M_TWOPI;
430 while( x[point_offset*i] > dstdefn->long_wrap_center + M_PI )
431 x[point_offset*i] -= M_TWOPI;
432 }
433 }
434
435 /* -------------------------------------------------------------------- */
436 /* Transform Z from meters if needed. */
437 /* -------------------------------------------------------------------- */
438 if( dstdefn->vto_meter != 1.0 && z != NULL )
439 {
440 for( i = 0; i < point_count; i++ )
441 z[point_offset*i] *= dstdefn->vfr_meter;
442 }
443
444 /* -------------------------------------------------------------------- */
445 /* Transform normalized axes into unusual output coordinate axis */
446 /* orientation if needed. */
447 /* -------------------------------------------------------------------- */
448 if( strcmp(dstdefn->axis,"enu") != 0 )
449 {
450 int err;
451
452 err = pj_adjust_axis( dstdefn->ctx, dstdefn->axis,
453 1, point_count, point_offset, x, y, z );
454 if( err != 0 )
455 return err;
456 }
457
458 return 0;
459 }
460
461 /************************************************************************/
462 /* pj_geodetic_to_geocentric() */
463 /************************************************************************/
464
pj_geodetic_to_geocentric(double a,double es,long point_count,int point_offset,double * x,double * y,double * z)465 int pj_geodetic_to_geocentric( double a, double es,
466 long point_count, int point_offset,
467 double *x, double *y, double *z )
468
469 {
470 double b;
471 int i;
472 GeocentricInfo gi;
473 int ret_errno = 0;
474
475 if( es == 0.0 )
476 b = a;
477 else
478 b = a * sqrt(1-es);
479
480 if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
481 {
482 return PJD_ERR_GEOCENTRIC;
483 }
484
485 for( i = 0; i < point_count; i++ )
486 {
487 long io = i * point_offset;
488
489 if( x[io] == HUGE_VAL )
490 continue;
491
492 if( pj_Convert_Geodetic_To_Geocentric( &gi, y[io], x[io], z[io],
493 x+io, y+io, z+io ) != 0 )
494 {
495 ret_errno = -14;
496 x[io] = y[io] = HUGE_VAL;
497 /* but keep processing points! */
498 }
499 }
500
501 return ret_errno;
502 }
503
504 /************************************************************************/
505 /* pj_geodetic_to_geocentric() */
506 /************************************************************************/
507
pj_geocentric_to_geodetic(double a,double es,long point_count,int point_offset,double * x,double * y,double * z)508 int pj_geocentric_to_geodetic( double a, double es,
509 long point_count, int point_offset,
510 double *x, double *y, double *z )
511
512 {
513 double b;
514 int i;
515 GeocentricInfo gi;
516
517 if( es == 0.0 )
518 b = a;
519 else
520 b = a * sqrt(1-es);
521
522 if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 )
523 {
524 return PJD_ERR_GEOCENTRIC;
525 }
526
527 for( i = 0; i < point_count; i++ )
528 {
529 long io = i * point_offset;
530
531 if( x[io] == HUGE_VAL )
532 continue;
533
534 pj_Convert_Geocentric_To_Geodetic( &gi, x[io], y[io], z[io],
535 y+io, x+io, z+io );
536 }
537
538 return 0;
539 }
540
541 /************************************************************************/
542 /* pj_compare_datums() */
543 /* */
544 /* Returns TRUE if the two datums are identical, otherwise */
545 /* FALSE. */
546 /************************************************************************/
547
pj_compare_datums(PJ * srcdefn,PJ * dstdefn)548 int pj_compare_datums( PJ *srcdefn, PJ *dstdefn )
549
550 {
551 if( srcdefn->datum_type != dstdefn->datum_type )
552 {
553 return 0;
554 }
555 else if( srcdefn->a_orig != dstdefn->a_orig
556 || ABS(srcdefn->es_orig - dstdefn->es_orig) > 0.000000000050 )
557 {
558 /* the tolerence for es is to ensure that GRS80 and WGS84 are
559 considered identical */
560 return 0;
561 }
562 else if( srcdefn->datum_type == PJD_3PARAM )
563 {
564 return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
565 && srcdefn->datum_params[1] == dstdefn->datum_params[1]
566 && srcdefn->datum_params[2] == dstdefn->datum_params[2]);
567 }
568 else if( srcdefn->datum_type == PJD_7PARAM )
569 {
570 return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
571 && srcdefn->datum_params[1] == dstdefn->datum_params[1]
572 && srcdefn->datum_params[2] == dstdefn->datum_params[2]
573 && srcdefn->datum_params[3] == dstdefn->datum_params[3]
574 && srcdefn->datum_params[4] == dstdefn->datum_params[4]
575 && srcdefn->datum_params[5] == dstdefn->datum_params[5]
576 && srcdefn->datum_params[6] == dstdefn->datum_params[6]);
577 }
578 else if( srcdefn->datum_type == PJD_GRIDSHIFT )
579 {
580 return strcmp( pj_param(srcdefn->ctx, srcdefn->params,"snadgrids").s,
581 pj_param(dstdefn->ctx, dstdefn->params,"snadgrids").s ) == 0;
582 }
583 else
584 return 1;
585 }
586
587 /************************************************************************/
588 /* pj_geocentic_to_wgs84() */
589 /************************************************************************/
590
pj_geocentric_to_wgs84(PJ * defn,long point_count,int point_offset,double * x,double * y,double * z)591 int pj_geocentric_to_wgs84( PJ *defn,
592 long point_count, int point_offset,
593 double *x, double *y, double *z )
594
595 {
596 int i;
597
598 if( defn->datum_type == PJD_3PARAM )
599 {
600 for( i = 0; i < point_count; i++ )
601 {
602 long io = i * point_offset;
603
604 if( x[io] == HUGE_VAL )
605 continue;
606
607 x[io] = x[io] + Dx_BF;
608 y[io] = y[io] + Dy_BF;
609 z[io] = z[io] + Dz_BF;
610 }
611 }
612 else if( defn->datum_type == PJD_7PARAM )
613 {
614 for( i = 0; i < point_count; i++ )
615 {
616 long io = i * point_offset;
617 double x_out, y_out, z_out;
618
619 if( x[io] == HUGE_VAL )
620 continue;
621
622 x_out = M_BF*( x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
623 y_out = M_BF*( Rz_BF*x[io] + y[io] - Rx_BF*z[io]) + Dy_BF;
624 z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] + z[io]) + Dz_BF;
625
626 x[io] = x_out;
627 y[io] = y_out;
628 z[io] = z_out;
629 }
630 }
631
632 return 0;
633 }
634
635 /************************************************************************/
636 /* pj_geocentic_from_wgs84() */
637 /************************************************************************/
638
pj_geocentric_from_wgs84(PJ * defn,long point_count,int point_offset,double * x,double * y,double * z)639 int pj_geocentric_from_wgs84( PJ *defn,
640 long point_count, int point_offset,
641 double *x, double *y, double *z )
642
643 {
644 int i;
645
646 if( defn->datum_type == PJD_3PARAM )
647 {
648 for( i = 0; i < point_count; i++ )
649 {
650 long io = i * point_offset;
651
652 if( x[io] == HUGE_VAL )
653 continue;
654
655 x[io] = x[io] - Dx_BF;
656 y[io] = y[io] - Dy_BF;
657 z[io] = z[io] - Dz_BF;
658 }
659 }
660 else if( defn->datum_type == PJD_7PARAM )
661 {
662 for( i = 0; i < point_count; i++ )
663 {
664 long io = i * point_offset;
665 double x_tmp, y_tmp, z_tmp;
666
667 if( x[io] == HUGE_VAL )
668 continue;
669
670 x_tmp = (x[io] - Dx_BF) / M_BF;
671 y_tmp = (y[io] - Dy_BF) / M_BF;
672 z_tmp = (z[io] - Dz_BF) / M_BF;
673
674 x[io] = x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp;
675 y[io] = -Rz_BF*x_tmp + y_tmp + Rx_BF*z_tmp;
676 z[io] = Ry_BF*x_tmp - Rx_BF*y_tmp + z_tmp;
677 }
678 }
679
680 return 0;
681 }
682
683 /************************************************************************/
684 /* pj_datum_transform() */
685 /* */
686 /* The input should be long/lat/z coordinates in radians in the */
687 /* source datum, and the output should be long/lat/z */
688 /* coordinates in radians in the destination datum. */
689 /************************************************************************/
690
pj_datum_transform(PJ * srcdefn,PJ * dstdefn,long point_count,int point_offset,double * x,double * y,double * z)691 int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
692 long point_count, int point_offset,
693 double *x, double *y, double *z )
694
695 {
696 double src_a, src_es, dst_a, dst_es;
697 int z_is_temp = FALSE;
698
699 /* -------------------------------------------------------------------- */
700 /* We cannot do any meaningful datum transformation if either */
701 /* the source or destination are of an unknown datum type */
702 /* (ie. only a +ellps declaration, no +datum). This is new */
703 /* behavior for PROJ 4.6.0. */
704 /* -------------------------------------------------------------------- */
705 if( srcdefn->datum_type == PJD_UNKNOWN
706 || dstdefn->datum_type == PJD_UNKNOWN )
707 return 0;
708
709 /* -------------------------------------------------------------------- */
710 /* Short cut if the datums are identical. */
711 /* -------------------------------------------------------------------- */
712 if( pj_compare_datums( srcdefn, dstdefn ) )
713 return 0;
714
715 src_a = srcdefn->a_orig;
716 src_es = srcdefn->es_orig;
717
718 dst_a = dstdefn->a_orig;
719 dst_es = dstdefn->es_orig;
720
721 /* -------------------------------------------------------------------- */
722 /* Create a temporary Z array if one is not provided. */
723 /* -------------------------------------------------------------------- */
724 if( z == NULL )
725 {
726 int bytes = sizeof(double) * point_count * point_offset;
727 z = (double *) pj_malloc(bytes);
728 memset( z, 0, bytes );
729 z_is_temp = TRUE;
730 }
731
732 #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; }}
733
734 /* -------------------------------------------------------------------- */
735 /* If this datum requires grid shifts, then apply it to geodetic */
736 /* coordinates. */
737 /* -------------------------------------------------------------------- */
738 if( srcdefn->datum_type == PJD_GRIDSHIFT )
739 {
740 pj_apply_gridshift_2( srcdefn, 0, point_count, point_offset, x, y, z );
741 CHECK_RETURN(srcdefn);
742
743 src_a = SRS_WGS84_SEMIMAJOR;
744 src_es = SRS_WGS84_ESQUARED;
745 }
746
747 if( dstdefn->datum_type == PJD_GRIDSHIFT )
748 {
749 dst_a = SRS_WGS84_SEMIMAJOR;
750 dst_es = SRS_WGS84_ESQUARED;
751 }
752
753 /* ==================================================================== */
754 /* Do we need to go through geocentric coordinates? */
755 /* ==================================================================== */
756 if( src_es != dst_es || src_a != dst_a
757 || srcdefn->datum_type == PJD_3PARAM
758 || srcdefn->datum_type == PJD_7PARAM
759 || dstdefn->datum_type == PJD_3PARAM
760 || dstdefn->datum_type == PJD_7PARAM)
761 {
762 /* -------------------------------------------------------------------- */
763 /* Convert to geocentric coordinates. */
764 /* -------------------------------------------------------------------- */
765 srcdefn->ctx->last_errno =
766 pj_geodetic_to_geocentric( src_a, src_es,
767 point_count, point_offset, x, y, z );
768 CHECK_RETURN(srcdefn);
769
770 /* -------------------------------------------------------------------- */
771 /* Convert between datums. */
772 /* -------------------------------------------------------------------- */
773 if( srcdefn->datum_type == PJD_3PARAM
774 || srcdefn->datum_type == PJD_7PARAM )
775 {
776 pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
777 CHECK_RETURN(srcdefn);
778 }
779
780 if( dstdefn->datum_type == PJD_3PARAM
781 || dstdefn->datum_type == PJD_7PARAM )
782 {
783 pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
784 CHECK_RETURN(dstdefn);
785 }
786
787 /* -------------------------------------------------------------------- */
788 /* Convert back to geodetic coordinates. */
789 /* -------------------------------------------------------------------- */
790 dstdefn->ctx->last_errno =
791 pj_geocentric_to_geodetic( dst_a, dst_es,
792 point_count, point_offset, x, y, z );
793 CHECK_RETURN(dstdefn);
794 }
795
796 /* -------------------------------------------------------------------- */
797 /* Apply grid shift to destination if required. */
798 /* -------------------------------------------------------------------- */
799 if( dstdefn->datum_type == PJD_GRIDSHIFT )
800 {
801 pj_apply_gridshift_2( dstdefn, 1, point_count, point_offset, x, y, z );
802 CHECK_RETURN(dstdefn);
803 }
804
805 if( z_is_temp )
806 pj_dalloc( z );
807
808 return 0;
809 }
810
811 /************************************************************************/
812 /* pj_adjust_axis() */
813 /* */
814 /* Normalize or de-normalized the x/y/z axes. The normal form */
815 /* is "enu" (easting, northing, up). */
816 /************************************************************************/
pj_adjust_axis(projCtx ctx,const char * axis,int denormalize_flag,long point_count,int point_offset,double * x,double * y,double * z)817 static int pj_adjust_axis( projCtx ctx,
818 const char *axis, int denormalize_flag,
819 long point_count, int point_offset,
820 double *x, double *y, double *z )
821
822 {
823 double x_in, y_in, z_in = 0.0;
824 int i, i_axis;
825
826 if( !denormalize_flag )
827 {
828 for( i = 0; i < point_count; i++ )
829 {
830 x_in = x[point_offset*i];
831 y_in = y[point_offset*i];
832 if( z )
833 z_in = z[point_offset*i];
834
835 for( i_axis = 0; i_axis < 3; i_axis++ )
836 {
837 double value;
838
839 if( i_axis == 0 )
840 value = x_in;
841 else if( i_axis == 1 )
842 value = y_in;
843 else
844 value = z_in;
845
846 switch( axis[i_axis] )
847 {
848 case 'e':
849 x[point_offset*i] = value; break;
850 case 'w':
851 x[point_offset*i] = -value; break;
852 case 'n':
853 y[point_offset*i] = value; break;
854 case 's':
855 y[point_offset*i] = -value; break;
856 case 'u':
857 if( z ) z[point_offset*i] = value; break;
858 case 'd':
859 if( z ) z[point_offset*i] = -value; break;
860 default:
861 pj_ctx_set_errno( ctx, PJD_ERR_AXIS );
862 return PJD_ERR_AXIS;
863 }
864 } /* i_axis */
865 } /* i (point) */
866 }
867
868 else /* denormalize */
869 {
870 for( i = 0; i < point_count; i++ )
871 {
872 x_in = x[point_offset*i];
873 y_in = y[point_offset*i];
874 if( z )
875 z_in = z[point_offset*i];
876
877 for( i_axis = 0; i_axis < 3; i_axis++ )
878 {
879 double *target;
880
881 if( i_axis == 2 && z == NULL )
882 continue;
883
884 if( i_axis == 0 )
885 target = x;
886 else if( i_axis == 1 )
887 target = y;
888 else
889 target = z;
890
891 switch( axis[i_axis] )
892 {
893 case 'e':
894 target[point_offset*i] = x_in; break;
895 case 'w':
896 target[point_offset*i] = -x_in; break;
897 case 'n':
898 target[point_offset*i] = y_in; break;
899 case 's':
900 target[point_offset*i] = -y_in; break;
901 case 'u':
902 target[point_offset*i] = z_in; break;
903 case 'd':
904 target[point_offset*i] = -z_in; break;
905 default:
906 pj_ctx_set_errno( ctx, PJD_ERR_AXIS );
907 return PJD_ERR_AXIS;
908 }
909 } /* i_axis */
910 } /* i (point) */
911 }
912
913 return 0;
914 }
915