1 /******************************************************************************
2 * Project: PROJ.4
3 * Purpose: Implement a (currently minimalistic) proj API based primarily
4 * on the PJ_COORD 4D geodetic spatiotemporal data type.
5 *
6 * Author: Thomas Knudsen, thokn@sdfe.dk, 2016-06-09/2016-11-06
7 *
8 ******************************************************************************
9 * Copyright (c) 2016, 2017 Thomas Knudsen/SDFE
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 #define FROM_PROJ_CPP
31
32 #include <assert.h>
33 #include <errno.h>
34 #include <stddef.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #ifndef _MSC_VER
39 #include <strings.h>
40 #endif
41
42 #include <algorithm>
43 #include <limits>
44
45 #include "proj.h"
46 #include "proj_experimental.h"
47 #include "proj_internal.h"
48 #include <math.h>
49 #include "geodesic.h"
50 #include "grids.hpp"
51 #include "filemanager.hpp"
52
53 #include "proj/common.hpp"
54 #include "proj/coordinateoperation.hpp"
55 #include "proj/internal/internal.hpp"
56 #include "proj/internal/io_internal.hpp"
57
58 using namespace NS_PROJ::internal;
59
60 /* Initialize PJ_COORD struct */
proj_coord(double x,double y,double z,double t)61 PJ_COORD proj_coord (double x, double y, double z, double t) {
62 PJ_COORD res;
63 res.v[0] = x;
64 res.v[1] = y;
65 res.v[2] = z;
66 res.v[3] = t;
67 return res;
68 }
69
opposite_direction(PJ_DIRECTION dir)70 static PJ_DIRECTION opposite_direction(PJ_DIRECTION dir) {
71 return static_cast<PJ_DIRECTION>(-dir);
72 }
73
74 /*****************************************************************************/
proj_angular_input(PJ * P,enum PJ_DIRECTION dir)75 int proj_angular_input (PJ *P, enum PJ_DIRECTION dir) {
76 /******************************************************************************
77 Returns 1 if the operator P expects angular input coordinates when
78 operating in direction dir, 0 otherwise.
79 dir: {PJ_FWD, PJ_INV}
80 ******************************************************************************/
81 if (PJ_FWD==dir)
82 return pj_left (P)==PJ_IO_UNITS_RADIANS;
83 return pj_right (P)==PJ_IO_UNITS_RADIANS;
84 }
85
86 /*****************************************************************************/
proj_angular_output(PJ * P,enum PJ_DIRECTION dir)87 int proj_angular_output (PJ *P, enum PJ_DIRECTION dir) {
88 /******************************************************************************
89 Returns 1 if the operator P provides angular output coordinates when
90 operating in direction dir, 0 otherwise.
91 dir: {PJ_FWD, PJ_INV}
92 ******************************************************************************/
93 return proj_angular_input (P, opposite_direction(dir));
94 }
95
96 /*****************************************************************************/
proj_degree_input(PJ * P,enum PJ_DIRECTION dir)97 int proj_degree_input (PJ *P, enum PJ_DIRECTION dir) {
98 /******************************************************************************
99 Returns 1 if the operator P expects degree input coordinates when
100 operating in direction dir, 0 otherwise.
101 dir: {PJ_FWD, PJ_INV}
102 ******************************************************************************/
103 if (PJ_FWD==dir)
104 return pj_left (P)==PJ_IO_UNITS_DEGREES;
105 return pj_right (P)==PJ_IO_UNITS_DEGREES;
106 }
107
108 /*****************************************************************************/
proj_degree_output(PJ * P,enum PJ_DIRECTION dir)109 int proj_degree_output (PJ *P, enum PJ_DIRECTION dir) {
110 /******************************************************************************
111 Returns 1 if the operator P provides degree output coordinates when
112 operating in direction dir, 0 otherwise.
113 dir: {PJ_FWD, PJ_INV}
114 ******************************************************************************/
115 return proj_degree_input (P, opposite_direction(dir));
116 }
117
118 /* Geodesic distance (in meter) + fwd and rev azimuth between two points on the ellipsoid */
proj_geod(const PJ * P,PJ_COORD a,PJ_COORD b)119 PJ_COORD proj_geod (const PJ *P, PJ_COORD a, PJ_COORD b) {
120 PJ_COORD c;
121 if( !P->geod ) {
122 return proj_coord_error();
123 }
124 /* Note: the geodesic code takes arguments in degrees */
125 geod_inverse (P->geod,
126 PJ_TODEG(a.lpz.phi), PJ_TODEG(a.lpz.lam),
127 PJ_TODEG(b.lpz.phi), PJ_TODEG(b.lpz.lam),
128 c.v, c.v+1, c.v+2
129 );
130
131 // cppcheck-suppress uninitvar
132 return c;
133 }
134
135
136 /* Geodesic distance (in meter) between two points with angular 2D coordinates */
proj_lp_dist(const PJ * P,PJ_COORD a,PJ_COORD b)137 double proj_lp_dist (const PJ *P, PJ_COORD a, PJ_COORD b) {
138 double s12, azi1, azi2;
139 /* Note: the geodesic code takes arguments in degrees */
140 if( !P->geod ) {
141 return HUGE_VAL;
142 }
143 geod_inverse (P->geod,
144 PJ_TODEG(a.lpz.phi), PJ_TODEG(a.lpz.lam),
145 PJ_TODEG(b.lpz.phi), PJ_TODEG(b.lpz.lam),
146 &s12, &azi1, &azi2
147 );
148 return s12;
149 }
150
151 /* The geodesic distance AND the vertical offset */
proj_lpz_dist(const PJ * P,PJ_COORD a,PJ_COORD b)152 double proj_lpz_dist (const PJ *P, PJ_COORD a, PJ_COORD b) {
153 if (HUGE_VAL==a.lpz.lam || HUGE_VAL==b.lpz.lam)
154 return HUGE_VAL;
155 return hypot (proj_lp_dist (P, a, b), a.lpz.z - b.lpz.z);
156 }
157
158 /* Euclidean distance between two points with linear 2D coordinates */
proj_xy_dist(PJ_COORD a,PJ_COORD b)159 double proj_xy_dist (PJ_COORD a, PJ_COORD b) {
160 return hypot (a.xy.x - b.xy.x, a.xy.y - b.xy.y);
161 }
162
163 /* Euclidean distance between two points with linear 3D coordinates */
proj_xyz_dist(PJ_COORD a,PJ_COORD b)164 double proj_xyz_dist (PJ_COORD a, PJ_COORD b) {
165 return hypot (proj_xy_dist (a, b), a.xyz.z - b.xyz.z);
166 }
167
168
169
170 /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
proj_roundtrip(PJ * P,PJ_DIRECTION direction,int n,PJ_COORD * coord)171 double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coord) {
172 int i;
173 PJ_COORD t, org;
174
175 if (nullptr==P)
176 return HUGE_VAL;
177
178 if (n < 1) {
179 proj_errno_set (P, EINVAL);
180 return HUGE_VAL;
181 }
182
183 /* in the first half-step, we generate the output value */
184 org = *coord;
185 *coord = proj_trans (P, direction, org);
186 t = *coord;
187
188 /* now we take n-1 full steps in inverse direction: We are */
189 /* out of phase due to the half step already taken */
190 for (i = 0; i < n - 1; i++)
191 t = proj_trans (P, direction, proj_trans (P, opposite_direction(direction), t) );
192
193 /* finally, we take the last half-step */
194 t = proj_trans (P, opposite_direction(direction), t);
195
196 /* checking for angular *input* since we do a roundtrip, and end where we begin */
197 if (proj_angular_input (P, direction))
198 return proj_lpz_dist (P, org, t);
199
200 return proj_xyz_dist (org, t);
201 }
202
203 /**************************************************************************************/
pj_get_suggested_operation(PJ_CONTEXT *,const std::vector<CoordOperation> & opList,const int iExcluded[2],PJ_DIRECTION direction,PJ_COORD coord)204 int pj_get_suggested_operation(PJ_CONTEXT*,
205 const std::vector<CoordOperation>& opList,
206 const int iExcluded[2],
207 PJ_DIRECTION direction,
208 PJ_COORD coord)
209 /**************************************************************************************/
210 {
211 // Select the operations that match the area of use
212 // and has the best accuracy.
213 int iBest = -1;
214 double bestAccuracy = std::numeric_limits<double>::max();
215 const int nOperations = static_cast<int>(opList.size());
216 for( int i = 0; i < nOperations; i++ ) {
217 if( i == iExcluded[0] || i == iExcluded[1] ) {
218 continue;
219 }
220 const auto &alt = opList[i];
221 bool spatialCriterionOK = false;
222 if( direction == PJ_FWD ) {
223 if( coord.xyzt.x >= alt.minxSrc &&
224 coord.xyzt.y >= alt.minySrc &&
225 coord.xyzt.x <= alt.maxxSrc &&
226 coord.xyzt.y <= alt.maxySrc) {
227 spatialCriterionOK = true;
228 }
229 } else {
230 if( coord.xyzt.x >= alt.minxDst &&
231 coord.xyzt.y >= alt.minyDst &&
232 coord.xyzt.x <= alt.maxxDst &&
233 coord.xyzt.y <= alt.maxyDst ) {
234 spatialCriterionOK = true;
235 }
236 }
237
238 if( spatialCriterionOK ) {
239 // The offshore test is for the "Test bug 245 (use +datum=carthage)"
240 // of testvarious. The long=10 lat=34 point belongs both to the
241 // onshore and offshore Tunisia area of uses, but is slightly
242 // onshore. So in a general way, prefer a onshore area to a
243 // offshore one.
244 if( iBest < 0 ||
245 (alt.accuracy >= 0 && alt.accuracy < bestAccuracy &&
246 !alt.isOffshore) ) {
247 iBest = i;
248 bestAccuracy = alt.accuracy;
249 }
250 }
251 }
252
253 return iBest;
254 }
255
256 /**************************************************************************************/
proj_trans(PJ * P,PJ_DIRECTION direction,PJ_COORD coord)257 PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coord) {
258 /***************************************************************************************
259 Apply the transformation P to the coordinate coord, preferring the 4D interfaces if
260 available.
261
262 See also pj_approx_2D_trans and pj_approx_3D_trans in pj_internal.c, which work
263 similarly, but prefers the 2D resp. 3D interfaces if available.
264 ***************************************************************************************/
265 if (nullptr==P || direction == PJ_IDENT)
266 return coord;
267 if (P->inverted)
268 direction = opposite_direction(direction);
269
270 if( !P->alternativeCoordinateOperations.empty() ) {
271 constexpr int N_MAX_RETRY = 2;
272 int iExcluded[N_MAX_RETRY] = {-1, -1};
273
274 const int nOperations = static_cast<int>(
275 P->alternativeCoordinateOperations.size());
276
277 // We may need several attempts. For example the point at
278 // lon=-111.5 lat=45.26 falls into the bounding box of the Canadian
279 // ntv2_0.gsb grid, except that it is not in any of the subgrids, being
280 // in the US. We thus need another retry that will select the conus
281 // grid.
282 for( int iRetry = 0; iRetry <= N_MAX_RETRY; iRetry++ )
283 {
284 // Do a first pass and select the operations that match the area of use
285 // and has the best accuracy.
286 int iBest = pj_get_suggested_operation(P->ctx,
287 P->alternativeCoordinateOperations,
288 iExcluded,
289 direction,
290 coord);
291 if( iBest < 0 ) {
292 break;
293 }
294 if( iRetry > 0 ) {
295 const int oldErrno = proj_errno_reset(P);
296 if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
297 pj_log(P->ctx, PJ_LOG_DEBUG, proj_errno_string(oldErrno));
298 }
299 pj_log(P->ctx, PJ_LOG_DEBUG,
300 "Did not result in valid result. "
301 "Attempting a retry with another operation.");
302 }
303
304 const auto& alt = P->alternativeCoordinateOperations[iBest];
305 if( P->iCurCoordOp != iBest ) {
306 if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
307 std::string msg("Using coordinate operation ");
308 msg += alt.name;
309 pj_log(P->ctx, PJ_LOG_DEBUG, msg.c_str());
310 }
311 P->iCurCoordOp = iBest;
312 }
313 PJ_COORD res = direction == PJ_FWD ?
314 pj_fwd4d( coord, alt.pj ) : pj_inv4d( coord, alt.pj );
315 if( proj_errno(alt.pj) == PJD_ERR_NETWORK_ERROR ) {
316 return proj_coord_error ();
317 }
318 if( res.xyzt.x != HUGE_VAL ) {
319 return res;
320 }
321 if( iRetry == N_MAX_RETRY ) {
322 break;
323 }
324 iExcluded[iRetry] = iBest;
325 }
326
327 // In case we did not find an operation whose area of use is compatible
328 // with the input coordinate, then goes through again the list, and
329 // use the first operation that does not require grids.
330 NS_PROJ::io::DatabaseContextPtr dbContext;
331 try
332 {
333 if( P->ctx->cpp_context ) {
334 dbContext = P->ctx->cpp_context->getDatabaseContext().as_nullable();
335 }
336 }
337 catch( const std::exception& ) {}
338 for( int i = 0; i < nOperations; i++ ) {
339 const auto &alt = P->alternativeCoordinateOperations[i];
340 auto coordOperation = dynamic_cast<
341 NS_PROJ::operation::CoordinateOperation*>(alt.pj->iso_obj.get());
342 if( coordOperation ) {
343 if( coordOperation->gridsNeeded(dbContext, true).empty() ) {
344 if( P->iCurCoordOp != i ) {
345 if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
346 std::string msg("Using coordinate operation ");
347 msg += alt.name;
348 pj_log(P->ctx, PJ_LOG_DEBUG, msg.c_str());
349 }
350 P->iCurCoordOp = i;
351 }
352 if( direction == PJ_FWD ) {
353 return pj_fwd4d( coord, alt.pj );
354 }
355 else {
356 return pj_inv4d( coord, alt.pj );
357 }
358 }
359 }
360 }
361
362 proj_errno_set (P, EINVAL);
363 return proj_coord_error ();
364 }
365
366 switch (direction) {
367 case PJ_FWD:
368 return pj_fwd4d (coord, P);
369 case PJ_INV:
370 return pj_inv4d (coord, P);
371 default:
372 break;
373 }
374
375 proj_errno_set (P, EINVAL);
376 return proj_coord_error ();
377 }
378
379
380
381 /*****************************************************************************/
proj_trans_array(PJ * P,PJ_DIRECTION direction,size_t n,PJ_COORD * coord)382 int proj_trans_array (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord) {
383 /******************************************************************************
384 Batch transform an array of PJ_COORD.
385
386 Returns 0 if all coordinates are transformed without error, otherwise
387 returns error number.
388 ******************************************************************************/
389 size_t i;
390
391 for (i = 0; i < n; i++) {
392 coord[i] = proj_trans (P, direction, coord[i]);
393 if (proj_errno(P))
394 return proj_errno (P);
395 }
396
397 return 0;
398 }
399
400
401
402 /*************************************************************************************/
proj_trans_generic(PJ * P,PJ_DIRECTION direction,double * x,size_t sx,size_t nx,double * y,size_t sy,size_t ny,double * z,size_t sz,size_t nz,double * t,size_t st,size_t nt)403 size_t proj_trans_generic (
404 PJ *P,
405 PJ_DIRECTION direction,
406 double *x, size_t sx, size_t nx,
407 double *y, size_t sy, size_t ny,
408 double *z, size_t sz, size_t nz,
409 double *t, size_t st, size_t nt
410 ) {
411 /**************************************************************************************
412
413 Transform a series of coordinates, where the individual coordinate dimension
414 may be represented by an array that is either
415
416 1. fully populated
417 2. a null pointer and/or a length of zero, which will be treated as a
418 fully populated array of zeroes
419 3. of length one, i.e. a constant, which will be treated as a fully
420 populated array of that constant value
421
422 The strides, sx, sy, sz, st, represent the step length, in bytes, between
423 consecutive elements of the corresponding array. This makes it possible for
424 proj_transform to handle transformation of a large class of application
425 specific data structures, without necessarily understanding the data structure
426 format, as in:
427
428 typedef struct {double x, y; int quality_level; char surveyor_name[134];} XYQS;
429 XYQS survey[345];
430 double height = 23.45;
431 PJ *P = {...};
432 size_t stride = sizeof (XYQS);
433 ...
434 proj_transform (
435 P, PJ_INV, sizeof(XYQS),
436 &(survey[0].x), stride, 345, (* We have 345 eastings *)
437 &(survey[0].y), stride, 345, (* ...and 345 northings. *)
438 &height, 1, (* The height is the constant 23.45 m *)
439 0, 0 (* and the time is the constant 0.00 s *)
440 );
441
442 This is similar to the inner workings of the pj_transform function, but the
443 stride functionality has been generalized to work for any size of basic unit,
444 not just a fixed number of doubles.
445
446 In most cases, the stride will be identical for x, y,z, and t, since they will
447 typically be either individual arrays (stride = sizeof(double)), or strided
448 views into an array of application specific data structures (stride = sizeof (...)).
449
450 But in order to support cases where x, y, z, and t come from heterogeneous
451 sources, individual strides, sx, sy, sz, st, are used.
452
453 Caveat: Since proj_transform does its work *in place*, this means that even the
454 supposedly constants (i.e. length 1 arrays) will return from the call in altered
455 state. Hence, remember to reinitialize between repeated calls.
456
457 Return value: Number of transformations completed.
458
459 **************************************************************************************/
460 PJ_COORD coord = {{0,0,0,0}};
461 size_t i, nmin;
462 double null_broadcast = 0;
463 double invalid_time = HUGE_VAL;
464
465 if (nullptr==P)
466 return 0;
467
468 if (P->inverted)
469 direction = opposite_direction(direction);
470
471 /* ignore lengths of null arrays */
472 if (nullptr==x) nx = 0;
473 if (nullptr==y) ny = 0;
474 if (nullptr==z) nz = 0;
475 if (nullptr==t) nt = 0;
476
477 /* and make the nullities point to some real world memory for broadcasting nulls */
478 if (0==nx) x = &null_broadcast;
479 if (0==ny) y = &null_broadcast;
480 if (0==nz) z = &null_broadcast;
481 if (0==nt) t = &invalid_time;
482
483 /* nothing to do? */
484 if (0==nx+ny+nz+nt)
485 return 0;
486
487 /* arrays of length 1 are constants, which we broadcast along the longer arrays */
488 /* so we need to find the length of the shortest non-unity array to figure out */
489 /* how many coordinate pairs we must transform */
490 nmin = (nx > 1)? nx: (ny > 1)? ny: (nz > 1)? nz: (nt > 1)? nt: 1;
491 if ((nx > 1) && (nx < nmin)) nmin = nx;
492 if ((ny > 1) && (ny < nmin)) nmin = ny;
493 if ((nz > 1) && (nz < nmin)) nmin = nz;
494 if ((nt > 1) && (nt < nmin)) nmin = nt;
495
496 /* Check validity of direction flag */
497 switch (direction) {
498 case PJ_FWD:
499 case PJ_INV:
500 break;
501 case PJ_IDENT:
502 return nmin;
503 default:
504 proj_errno_set (P, EINVAL);
505 return 0;
506 }
507
508 /* Arrays of length==0 are broadcast as the constant 0 */
509 /* Arrays of length==1 are broadcast as their single value */
510 /* Arrays of length >1 are iterated over (for the first nmin values) */
511 /* The slightly convolved incremental indexing is used due */
512 /* to the stride, which may be any size supported by the platform */
513 for (i = 0; i < nmin; i++) {
514 coord.xyzt.x = *x;
515 coord.xyzt.y = *y;
516 coord.xyzt.z = *z;
517 coord.xyzt.t = *t;
518
519 coord = proj_trans(P, direction, coord);
520
521 /* in all full length cases, we overwrite the input with the output, */
522 /* and step on to the next element. */
523 /* The casts are somewhat funky, but they compile down to no-ops and */
524 /* they tell compilers and static analyzers that we know what we do */
525 if (nx > 1) {
526 *x = coord.xyzt.x;
527 x = (double *) ((void *) ( ((char *) x) + sx));
528 }
529 if (ny > 1) {
530 *y = coord.xyzt.y;
531 y = (double *) ((void *) ( ((char *) y) + sy));
532 }
533 if (nz > 1) {
534 *z = coord.xyzt.z;
535 z = (double *) ((void *) ( ((char *) z) + sz));
536 }
537 if (nt > 1) {
538 *t = coord.xyzt.t;
539 t = (double *) ((void *) ( ((char *) t) + st));
540 }
541 }
542
543 /* Last time around, we update the length 1 cases with their transformed alter egos */
544 if (nx==1)
545 *x = coord.xyzt.x;
546 if (ny==1)
547 *y = coord.xyzt.y;
548 if (nz==1)
549 *z = coord.xyzt.z;
550 if (nt==1)
551 *t = coord.xyzt.t;
552
553 return i;
554 }
555
556
557 /*************************************************************************************/
pj_geocentric_latitude(const PJ * P,PJ_DIRECTION direction,PJ_COORD coord)558 PJ_COORD pj_geocentric_latitude (const PJ *P, PJ_DIRECTION direction, PJ_COORD coord) {
559 /**************************************************************************************
560 Convert geographical latitude to geocentric (or the other way round if
561 direction = PJ_INV)
562
563 The conversion involves a call to the tangent function, which goes through the
564 roof at the poles, so very close (the last centimeter) to the poles no
565 conversion takes place and the input latitude is copied directly to the output.
566
567 Fortunately, the geocentric latitude converges to the geographical at the
568 poles, so the difference is negligible.
569
570 For the spherical case, the geographical latitude equals the geocentric, and
571 consequently, the input is copied directly to the output.
572 **************************************************************************************/
573 const double limit = M_HALFPI - 1e-9;
574 PJ_COORD res = coord;
575 if ((coord.lp.phi > limit) || (coord.lp.phi < -limit) || (P->es==0))
576 return res;
577 if (direction==PJ_FWD)
578 res.lp.phi = atan (P->one_es * tan (coord.lp.phi) );
579 else
580 res.lp.phi = atan (P->rone_es * tan (coord.lp.phi) );
581
582 return res;
583 }
584
proj_torad(double angle_in_degrees)585 double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);}
proj_todeg(double angle_in_radians)586 double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);}
587
proj_dmstor(const char * is,char ** rs)588 double proj_dmstor(const char *is, char **rs) {
589 return dmstor(is, rs);
590 }
591
proj_rtodms(char * s,double r,int pos,int neg)592 char* proj_rtodms(char *s, double r, int pos, int neg) {
593 return rtodms(s, r, pos, neg);
594 }
595
596 /*************************************************************************************/
skip_prep_fin(PJ * P)597 static PJ* skip_prep_fin(PJ *P) {
598 /**************************************************************************************
599 Skip prepare and finalize function for the various "helper operations" added to P when
600 in cs2cs compatibility mode.
601 **************************************************************************************/
602 P->skip_fwd_prepare = 1;
603 P->skip_fwd_finalize = 1;
604 P->skip_inv_prepare = 1;
605 P->skip_inv_finalize = 1;
606 return P;
607 }
608
609 /*************************************************************************************/
cs2cs_emulation_setup(PJ * P)610 static int cs2cs_emulation_setup (PJ *P) {
611 /**************************************************************************************
612 If any cs2cs style modifiers are given (axis=..., towgs84=..., ) create the 4D API
613 equivalent operations, so the preparation and finalization steps in the pj_inv/pj_fwd
614 invocators can emulate the behavior of pj_transform and the cs2cs app.
615
616 Returns 1 on success, 0 on failure
617 **************************************************************************************/
618 PJ *Q;
619 paralist *p;
620 int do_cart = 0;
621 if (nullptr==P)
622 return 0;
623
624 /* Don't recurse when calling proj_create (which calls us back) */
625 if (pj_param_exists (P->params, "break_cs2cs_recursion"))
626 return 1;
627
628 /* Swap axes? */
629 p = pj_param_exists (P->params, "axis");
630
631 const bool disable_grid_presence_check = pj_param_exists (
632 P->params, "disable_grid_presence_check") != nullptr;
633
634 /* Don't axisswap if data are already in "enu" order */
635 if (p && (0!=strcmp ("enu", p->param))) {
636 char *def = static_cast<char*>(malloc (100+strlen(P->axis)));
637 if (nullptr==def)
638 return 0;
639 sprintf (def, "break_cs2cs_recursion proj=axisswap axis=%s", P->axis);
640 Q = pj_create_internal (P->ctx, def);
641 free (def);
642 if (nullptr==Q)
643 return 0;
644 P->axisswap = skip_prep_fin(Q);
645 }
646
647 /* Geoid grid(s) given? */
648 p = pj_param_exists (P->params, "geoidgrids");
649 if (!disable_grid_presence_check && p && strlen (p->param) > strlen ("geoidgrids=")) {
650 char *gridnames = p->param + strlen ("geoidgrids=");
651 char *def = static_cast<char*>(malloc (100+2*strlen(gridnames)));
652 if (nullptr==def)
653 return 0;
654 sprintf (def, "break_cs2cs_recursion proj=vgridshift grids=%s",
655 pj_double_quote_string_param_if_needed(gridnames).c_str());
656 Q = pj_create_internal (P->ctx, def);
657 free (def);
658 if (nullptr==Q)
659 return 0;
660 P->vgridshift = skip_prep_fin(Q);
661 }
662
663 /* Datum shift grid(s) given? */
664 p = pj_param_exists (P->params, "nadgrids");
665 if (!disable_grid_presence_check && p && strlen (p->param) > strlen ("nadgrids=")) {
666 char *gridnames = p->param + strlen ("nadgrids=");
667 char *def = static_cast<char*>(malloc (100+2*strlen(gridnames)));
668 if (nullptr==def)
669 return 0;
670 sprintf (def, "break_cs2cs_recursion proj=hgridshift grids=%s",
671 pj_double_quote_string_param_if_needed(gridnames).c_str());
672 Q = pj_create_internal (P->ctx, def);
673 free (def);
674 if (nullptr==Q)
675 return 0;
676 P->hgridshift = skip_prep_fin(Q);
677 }
678
679 /* We ignore helmert if we have grid shift */
680 p = P->hgridshift ? nullptr : pj_param_exists (P->params, "towgs84");
681 while (p) {
682 char *def;
683 char *s = p->param;
684 double *d = P->datum_params;
685 size_t n = strlen (s);
686
687 /* We ignore null helmert shifts (common in auto-translated resource files, e.g. epsg) */
688 if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6]) {
689 /* If the current ellipsoid is not WGS84, then make sure the */
690 /* change in ellipsoid is still done. */
691 if (!(fabs(P->a_orig - 6378137.0) < 1e-8 && fabs(P->es_orig - 0.0066943799901413) < 1e-15)) {
692 do_cart = 1;
693 }
694 break;
695 }
696
697 if (n <= 8) /* 8==strlen ("towgs84=") */
698 return 0;
699
700 def = static_cast<char*>(malloc (100+n));
701 if (nullptr==def)
702 return 0;
703 sprintf (def, "break_cs2cs_recursion proj=helmert exact %s convention=position_vector", s);
704 Q = pj_create_internal (P->ctx, def);
705 free(def);
706 if (nullptr==Q)
707 return 0;
708 pj_inherit_ellipsoid_def (P, Q);
709 P->helmert = skip_prep_fin (Q);
710
711 break;
712 }
713
714 /* We also need cartesian/geographical transformations if we are working in */
715 /* geocentric/cartesian space or we need to do a Helmert transform. */
716 if (P->is_geocent || P->helmert || do_cart) {
717 char def[150];
718 sprintf (def, "break_cs2cs_recursion proj=cart a=%40.20g es=%40.20g", P->a_orig, P->es_orig);
719 {
720 /* In case the current locale does not use dot but comma as decimal */
721 /* separator, replace it with dot, so that proj_atof() behaves */
722 /* correctly. */
723 /* TODO later: use C++ ostringstream with imbue(std::locale::classic()) */
724 /* to be locale unaware */
725 char* next_pos;
726 for (next_pos = def; (next_pos = strchr (next_pos, ',')) != nullptr; next_pos++) {
727 *next_pos = '.';
728 }
729 }
730 Q = pj_create_internal (P->ctx, def);
731 if (nullptr==Q)
732 return 0;
733 P->cart = skip_prep_fin (Q);
734
735 if (!P->is_geocent) {
736 sprintf (def, "break_cs2cs_recursion proj=cart ellps=WGS84");
737 Q = pj_create_internal (P->ctx, def);
738 if (nullptr==Q)
739 return 0;
740 P->cart_wgs84 = skip_prep_fin (Q);
741 }
742 }
743
744 return 1;
745 }
746
747
748 /*************************************************************************************/
pj_create_internal(PJ_CONTEXT * ctx,const char * definition)749 PJ *pj_create_internal (PJ_CONTEXT *ctx, const char *definition) {
750 /*************************************************************************************/
751
752 /**************************************************************************************
753 Create a new PJ object in the context ctx, using the given definition. If ctx==0,
754 the default context is used, if definition==0, or invalid, a null-pointer is
755 returned. The definition may use '+' as argument start indicator, as in
756 "+proj=utm +zone=32", or leave it out, as in "proj=utm zone=32".
757
758 It may even use free formatting "proj = utm; zone =32 ellps= GRS80".
759 Note that the semicolon separator is allowed, but not required.
760 **************************************************************************************/
761 PJ *P;
762 char *args, **argv;
763 size_t argc, n;
764 int ret;
765 int allow_init_epsg;
766
767 if (nullptr==ctx)
768 ctx = pj_get_default_ctx ();
769
770 /* Make a copy that we can manipulate */
771 n = strlen (definition);
772 args = (char *) malloc (n + 1);
773 if (nullptr==args) {
774 proj_context_errno_set(ctx, ENOMEM);
775 return nullptr;
776 }
777 strcpy (args, definition);
778
779 argc = pj_trim_argc (args);
780 if (argc==0) {
781 pj_dealloc (args);
782 proj_context_errno_set(ctx, PJD_ERR_NO_ARGS);
783 return nullptr;
784 }
785
786 argv = pj_trim_argv (argc, args);
787 if (!argv) {
788 pj_dealloc(args);
789 proj_context_errno_set(ctx, ENOMEM);
790 return nullptr;
791 }
792
793 /* ...and let pj_init_ctx do the hard work */
794 /* New interface: forbid init=epsg:XXXX syntax by default */
795 allow_init_epsg = proj_context_get_use_proj4_init_rules(ctx, FALSE);
796 P = pj_init_ctx_with_allow_init_epsg (ctx, (int) argc, argv, allow_init_epsg);
797
798 pj_dealloc (argv);
799 pj_dealloc (args);
800
801 /* Support cs2cs-style modifiers */
802 ret = cs2cs_emulation_setup (P);
803 if (0==ret)
804 return proj_destroy (P);
805
806 return P;
807 }
808
809 /*************************************************************************************/
proj_create_argv(PJ_CONTEXT * ctx,int argc,char ** argv)810 PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv) {
811 /**************************************************************************************
812 Create a new PJ object in the context ctx, using the given definition argument
813 array argv. If ctx==0, the default context is used, if definition==0, or invalid,
814 a null-pointer is returned. The definition arguments may use '+' as argument start
815 indicator, as in {"+proj=utm", "+zone=32"}, or leave it out, as in {"proj=utm",
816 "zone=32"}.
817 **************************************************************************************/
818 PJ *P;
819 const char *c;
820
821 if (nullptr==ctx)
822 ctx = pj_get_default_ctx ();
823 if (nullptr==argv) {
824 proj_context_errno_set(ctx, PJD_ERR_NO_ARGS);
825 return nullptr;
826 }
827
828 /* We assume that free format is used, and build a full proj_create compatible string */
829 c = pj_make_args (argc, argv);
830 if (nullptr==c) {
831 proj_context_errno_set(ctx, ENOMEM);
832 return nullptr;
833 }
834
835 P = proj_create (ctx, c);
836
837 pj_dealloc ((char *) c);
838 return P;
839 }
840
841 /*************************************************************************************/
pj_create_argv_internal(PJ_CONTEXT * ctx,int argc,char ** argv)842 PJ *pj_create_argv_internal (PJ_CONTEXT *ctx, int argc, char **argv) {
843 /**************************************************************************************
844 Same as proj_create_argv() but calls pj_create_internal() instead of proj_create() internally
845 **************************************************************************************/
846 PJ *P;
847 const char *c;
848
849 if (nullptr==ctx)
850 ctx = pj_get_default_ctx ();
851 if (nullptr==argv) {
852 proj_context_errno_set(ctx, PJD_ERR_NO_ARGS);
853 return nullptr;
854 }
855
856 /* We assume that free format is used, and build a full proj_create compatible string */
857 c = pj_make_args (argc, argv);
858 if (nullptr==c) {
859 proj_context_errno_set(ctx, ENOMEM);
860 return nullptr;
861 }
862
863 P = pj_create_internal (ctx, c);
864
865 pj_dealloc ((char *) c);
866 return P;
867 }
868
869 /** Create an area of use */
proj_area_create(void)870 PJ_AREA * proj_area_create(void) {
871 return static_cast<PJ_AREA*>(pj_calloc(1, sizeof(PJ_AREA)));
872 }
873
874 /** Assign a bounding box to an area of use. */
proj_area_set_bbox(PJ_AREA * area,double west_lon_degree,double south_lat_degree,double east_lon_degree,double north_lat_degree)875 void proj_area_set_bbox(PJ_AREA *area,
876 double west_lon_degree,
877 double south_lat_degree,
878 double east_lon_degree,
879 double north_lat_degree) {
880 area->bbox_set = TRUE;
881 area->west_lon_degree = west_lon_degree;
882 area->south_lat_degree = south_lat_degree;
883 area->east_lon_degree = east_lon_degree;
884 area->north_lat_degree = north_lat_degree;
885 }
886
887 /** Free an area of use */
proj_area_destroy(PJ_AREA * area)888 void proj_area_destroy(PJ_AREA* area) {
889 pj_dealloc(area);
890 }
891
892 /************************************************************************/
893 /* proj_context_use_proj4_init_rules() */
894 /************************************************************************/
895
proj_context_use_proj4_init_rules(PJ_CONTEXT * ctx,int enable)896 void proj_context_use_proj4_init_rules(PJ_CONTEXT *ctx, int enable) {
897 if( ctx == nullptr ) {
898 ctx = pj_get_default_ctx();
899 }
900 ctx->use_proj4_init_rules = enable;
901 }
902
903 /************************************************************************/
904 /* EQUAL() */
905 /************************************************************************/
906
EQUAL(const char * a,const char * b)907 static int EQUAL(const char* a, const char* b) {
908 #ifdef _MSC_VER
909 return _stricmp(a, b) == 0;
910 #else
911 return strcasecmp(a, b) == 0;
912 #endif
913 }
914
915 /************************************************************************/
916 /* proj_context_get_use_proj4_init_rules() */
917 /************************************************************************/
918
proj_context_get_use_proj4_init_rules(PJ_CONTEXT * ctx,int from_legacy_code_path)919 int proj_context_get_use_proj4_init_rules(PJ_CONTEXT *ctx, int from_legacy_code_path) {
920 const char* val = getenv("PROJ_USE_PROJ4_INIT_RULES");
921
922 if( ctx == nullptr ) {
923 ctx = pj_get_default_ctx();
924 }
925
926 if( val ) {
927 if( EQUAL(val, "yes") || EQUAL(val, "on") || EQUAL(val, "true") ) {
928 return TRUE;
929 }
930 if( EQUAL(val, "no") || EQUAL(val, "off") || EQUAL(val, "false") ) {
931 return FALSE;
932 }
933 pj_log(ctx, PJ_LOG_ERROR, "Invalid value for PROJ_USE_PROJ4_INIT_RULES");
934 }
935
936 if( ctx->use_proj4_init_rules >= 0 ) {
937 return ctx->use_proj4_init_rules;
938 }
939 return from_legacy_code_path;
940 }
941
942 /** Adds a " +type=crs" suffix to a PROJ string (if it is a PROJ string) */
pj_add_type_crs_if_needed(const std::string & str)943 std::string pj_add_type_crs_if_needed(const std::string& str)
944 {
945 std::string ret(str);
946 if( (starts_with(str, "proj=") ||
947 starts_with(str, "+proj=") ||
948 starts_with(str, "+init=") ||
949 starts_with(str, "+title=")) &&
950 str.find("type=crs") == std::string::npos )
951 {
952 ret += " +type=crs";
953 }
954 return ret;
955 }
956
957 /*****************************************************************************/
reproject_bbox(PJ * pjGeogToCrs,double west_lon,double south_lat,double east_lon,double north_lat,double & minx,double & miny,double & maxx,double & maxy)958 static void reproject_bbox(PJ* pjGeogToCrs,
959 double west_lon, double south_lat,
960 double east_lon, double north_lat,
961 double& minx,
962 double& miny,
963 double& maxx,
964 double& maxy) {
965 /*****************************************************************************/
966
967 minx = -std::numeric_limits<double>::max();
968 miny = -std::numeric_limits<double>::max();
969 maxx = std::numeric_limits<double>::max();
970 maxy = std::numeric_limits<double>::max();
971
972 if( !(west_lon == -180.0 && east_lon == 180.0 &&
973 south_lat == -90.0 && north_lat == 90.0) )
974 {
975 minx = -minx;
976 miny = -miny;
977 maxx = -maxx;
978 maxy = -maxy;
979
980 std::vector<double> x(21 * 4), y(21 * 4);
981 for( int j = 0; j <= 20; j++ )
982 {
983 x[j] = west_lon + j * (east_lon - west_lon) / 20;
984 y[j] = south_lat;
985 x[21+j] = west_lon + j * (east_lon - west_lon) / 20;
986 y[21+j] = north_lat;
987 x[21*2+j] = west_lon;
988 y[21*2+j] = south_lat + j * (north_lat - south_lat) / 20;
989 x[21*3+j] = east_lon;
990 y[21*3+j] = south_lat + j * (north_lat - south_lat) / 20;
991 }
992 proj_trans_generic (
993 pjGeogToCrs, PJ_FWD,
994 &x[0], sizeof(double), 21 * 4,
995 &y[0], sizeof(double), 21 * 4,
996 nullptr, 0, 0,
997 nullptr, 0, 0);
998 for( int j = 0; j < 21 * 4; j++ )
999 {
1000 if( x[j] != HUGE_VAL && y[j] != HUGE_VAL )
1001 {
1002 minx = std::min(minx, x[j]);
1003 miny = std::min(miny, y[j]);
1004 maxx = std::max(maxx, x[j]);
1005 maxy = std::max(maxy, y[j]);
1006 }
1007 }
1008 }
1009 }
1010
1011
1012 /*****************************************************************************/
add_coord_op_to_list(int idxInOriginalList,PJ * op,double west_lon,double south_lat,double east_lon,double north_lat,PJ * pjGeogToSrc,PJ * pjGeogToDst,bool isOffshore,std::vector<CoordOperation> & altCoordOps)1013 static PJ* add_coord_op_to_list(
1014 int idxInOriginalList,
1015 PJ* op,
1016 double west_lon, double south_lat,
1017 double east_lon, double north_lat,
1018 PJ* pjGeogToSrc,
1019 PJ* pjGeogToDst,
1020 bool isOffshore,
1021 std::vector<CoordOperation>& altCoordOps) {
1022 /*****************************************************************************/
1023
1024 double minxSrc;
1025 double minySrc;
1026 double maxxSrc;
1027 double maxySrc;
1028 double minxDst;
1029 double minyDst;
1030 double maxxDst;
1031 double maxyDst;
1032
1033 reproject_bbox(pjGeogToSrc, west_lon, south_lat, east_lon, north_lat,
1034 minxSrc, minySrc, maxxSrc, maxySrc);
1035 reproject_bbox(pjGeogToDst, west_lon, south_lat, east_lon, north_lat,
1036 minxDst, minyDst, maxxDst, maxyDst);
1037
1038 if( minxSrc <= maxxSrc && minxDst <= maxxDst )
1039 {
1040 const char* c_name = proj_get_name(op);
1041 std::string name(c_name ? c_name : "");
1042
1043 const double accuracy = proj_coordoperation_get_accuracy(op->ctx, op);
1044 altCoordOps.emplace_back(idxInOriginalList,
1045 minxSrc, minySrc, maxxSrc, maxySrc,
1046 minxDst, minyDst, maxxDst, maxyDst,
1047 op, name, accuracy, isOffshore);
1048 op = nullptr;
1049 }
1050 return op;
1051 }
1052
1053 /*****************************************************************************/
create_operation_to_geog_crs(PJ_CONTEXT * ctx,const PJ * crs)1054 static PJ* create_operation_to_geog_crs(PJ_CONTEXT* ctx, const PJ* crs) {
1055 /*****************************************************************************/
1056 // Create a geographic 2D long-lat degrees CRS that is related to the
1057 // CRS
1058 auto geodetic_crs = proj_crs_get_geodetic_crs(ctx, crs);
1059 if( !geodetic_crs ) {
1060 proj_context_log_debug(ctx, "Cannot find geodetic CRS matching CRS");
1061 return nullptr;
1062 }
1063
1064 auto geodetic_crs_type = proj_get_type(geodetic_crs);
1065 if( geodetic_crs_type == PJ_TYPE_GEOCENTRIC_CRS ||
1066 geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
1067 geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_3D_CRS )
1068 {
1069 auto datum = proj_crs_get_datum(ctx, geodetic_crs);
1070 auto datum_ensemble = proj_crs_get_datum_ensemble(ctx, geodetic_crs);
1071 auto cs = proj_create_ellipsoidal_2D_cs(
1072 ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, nullptr, 0);
1073 auto temp = proj_create_geographic_crs_from_datum(
1074 ctx, "unnamed crs", datum ? datum : datum_ensemble,
1075 cs);
1076 proj_destroy(datum);
1077 proj_destroy(datum_ensemble);
1078 proj_destroy(cs);
1079 proj_destroy(geodetic_crs);
1080 geodetic_crs = temp;
1081 geodetic_crs_type = proj_get_type(geodetic_crs);
1082 }
1083 if( geodetic_crs_type != PJ_TYPE_GEOGRAPHIC_2D_CRS )
1084 {
1085 // Shouldn't happen
1086 proj_context_log_debug(ctx, "Cannot find geographic CRS matching CRS");
1087 proj_destroy(geodetic_crs);
1088 return nullptr;
1089 }
1090
1091 // Create the transformation from this geographic 2D CRS to the source CRS
1092 auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
1093 proj_operation_factory_context_set_spatial_criterion(
1094 ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
1095 proj_operation_factory_context_set_grid_availability_use(
1096 ctx, operation_ctx, PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
1097 auto target_crs_2D = proj_crs_demote_to_2D(ctx, nullptr, crs);
1098 auto op_list_to_geodetic = proj_create_operations(
1099 ctx, geodetic_crs, target_crs_2D, operation_ctx);
1100 proj_destroy(target_crs_2D);
1101 proj_operation_factory_context_destroy(operation_ctx);
1102 proj_destroy(geodetic_crs);
1103
1104 const int nOpCount = op_list_to_geodetic == nullptr ? 0 :
1105 proj_list_get_count(op_list_to_geodetic);
1106 if( nOpCount == 0 )
1107 {
1108 proj_context_log_debug(ctx, "Cannot compute transformation from geographic CRS to CRS");
1109 proj_list_destroy(op_list_to_geodetic);
1110 return nullptr;
1111 }
1112 PJ* opGeogToCrs = nullptr;
1113 // Use in priority operations *without* grids
1114 for(int i = 0; i < nOpCount; i++ )
1115 {
1116 auto op = proj_list_get(ctx, op_list_to_geodetic, i);
1117 assert(op);
1118 if( proj_coordoperation_get_grid_used_count(ctx, op) == 0 )
1119 {
1120 opGeogToCrs = op;
1121 break;
1122 }
1123 proj_destroy(op);
1124 }
1125 if( opGeogToCrs == nullptr )
1126 {
1127 opGeogToCrs = proj_list_get(ctx, op_list_to_geodetic, 0);
1128 assert(opGeogToCrs);
1129 }
1130 proj_list_destroy(op_list_to_geodetic);
1131 return opGeogToCrs;
1132 }
1133
1134 /*****************************************************************************/
proj_create_crs_to_crs(PJ_CONTEXT * ctx,const char * source_crs,const char * target_crs,PJ_AREA * area)1135 PJ *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *source_crs, const char *target_crs, PJ_AREA *area) {
1136 /******************************************************************************
1137 Create a transformation pipeline between two known coordinate reference
1138 systems.
1139
1140 See docs/source/development/reference/functions.rst
1141
1142 ******************************************************************************/
1143 if( !ctx ) {
1144 ctx = pj_get_default_ctx();
1145 }
1146
1147 PJ* src;
1148 PJ* dst;
1149 try
1150 {
1151 std::string source_crs_modified(pj_add_type_crs_if_needed(source_crs));
1152 std::string target_crs_modified(pj_add_type_crs_if_needed(target_crs));
1153
1154 src = proj_create(ctx, source_crs_modified.c_str());
1155 if( !src ) {
1156 proj_context_log_debug(ctx, "Cannot instantiate source_crs");
1157 return nullptr;
1158 }
1159
1160 dst = proj_create(ctx, target_crs_modified.c_str());
1161 if( !dst ) {
1162 proj_context_log_debug(ctx, "Cannot instantiate target_crs");
1163 proj_destroy(src);
1164 return nullptr;
1165 }
1166 }
1167 catch( const std::exception& )
1168 {
1169 return nullptr;
1170 }
1171
1172 auto ret = proj_create_crs_to_crs_from_pj(ctx, src, dst, area, nullptr);
1173 proj_destroy(src);
1174 proj_destroy(dst);
1175 return ret;
1176 }
1177
1178
1179 /*****************************************************************************/
pj_create_prepared_operations(PJ_CONTEXT * ctx,const PJ * source_crs,const PJ * target_crs,PJ_OBJ_LIST * op_list)1180 std::vector<CoordOperation> pj_create_prepared_operations(PJ_CONTEXT *ctx,
1181 const PJ *source_crs,
1182 const PJ *target_crs,
1183 PJ_OBJ_LIST* op_list)
1184 /*****************************************************************************/
1185 {
1186 auto pjGeogToSrc = create_operation_to_geog_crs(ctx, source_crs);
1187 if( !pjGeogToSrc )
1188 {
1189 proj_context_log_debug(ctx,
1190 "Cannot create transformation from geographic CRS of source CRS to source CRS");
1191 return {};
1192 }
1193
1194 auto pjGeogToDst = create_operation_to_geog_crs(ctx, target_crs);
1195 if( !pjGeogToDst )
1196 {
1197 proj_context_log_debug(ctx,
1198 "Cannot create transformation from geographic CRS of target CRS to target CRS");
1199 proj_destroy(pjGeogToSrc);
1200 return {};
1201 }
1202
1203 try
1204 {
1205 std::vector<CoordOperation> preparedOpList;
1206
1207 // Iterate over source->target candidate transformations and reproject
1208 // their long-lat bounding box into the source CRS.
1209 const auto op_count = proj_list_get_count(op_list);
1210 for( int i = 0; i < op_count; i++ )
1211 {
1212 auto op = proj_list_get(ctx, op_list, i);
1213 assert(op);
1214 double west_lon = 0.0;
1215 double south_lat = 0.0;
1216 double east_lon = 0.0;
1217 double north_lat = 0.0;
1218
1219 const char* areaName = nullptr;
1220 if( proj_get_area_of_use(ctx, op,
1221 &west_lon, &south_lat, &east_lon, &north_lat, &areaName) )
1222 {
1223 const bool isOffshore =
1224 areaName && strstr(areaName, "- offshore");
1225 if( west_lon <= east_lon )
1226 {
1227 op = add_coord_op_to_list(i, op,
1228 west_lon, south_lat, east_lon, north_lat,
1229 pjGeogToSrc, pjGeogToDst, isOffshore,
1230 preparedOpList);
1231 }
1232 else
1233 {
1234 auto op_clone = proj_clone(ctx, op);
1235
1236 op = add_coord_op_to_list(i, op,
1237 west_lon, south_lat, 180, north_lat,
1238 pjGeogToSrc, pjGeogToDst, isOffshore,
1239 preparedOpList);
1240 op_clone = add_coord_op_to_list(i, op_clone,
1241 -180, south_lat, east_lon, north_lat,
1242 pjGeogToSrc, pjGeogToDst, isOffshore,
1243 preparedOpList);
1244 proj_destroy(op_clone);
1245 }
1246 }
1247
1248 proj_destroy(op);
1249 }
1250
1251 proj_destroy(pjGeogToSrc);
1252 proj_destroy(pjGeogToDst);
1253 return preparedOpList;
1254 }
1255 catch( const std::exception& )
1256 {
1257 proj_destroy(pjGeogToSrc);
1258 proj_destroy(pjGeogToDst);
1259 return {};
1260 }
1261 }
1262
1263 /*****************************************************************************/
proj_create_crs_to_crs_from_pj(PJ_CONTEXT * ctx,const PJ * source_crs,const PJ * target_crs,PJ_AREA * area,const char * const *)1264 PJ *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, const PJ *target_crs, PJ_AREA *area, const char* const *) {
1265 /******************************************************************************
1266 Create a transformation pipeline between two known coordinate reference
1267 systems.
1268
1269 See docs/source/development/reference/functions.rst
1270
1271 ******************************************************************************/
1272 if( !ctx ) {
1273 ctx = pj_get_default_ctx();
1274 }
1275
1276 auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
1277 if( !operation_ctx ) {
1278 return nullptr;
1279 }
1280
1281 if( area && area->bbox_set ) {
1282 proj_operation_factory_context_set_area_of_interest(
1283 ctx,
1284 operation_ctx,
1285 area->west_lon_degree,
1286 area->south_lat_degree,
1287 area->east_lon_degree,
1288 area->north_lat_degree);
1289 }
1290
1291 proj_operation_factory_context_set_spatial_criterion(
1292 ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
1293 proj_operation_factory_context_set_grid_availability_use(
1294 ctx, operation_ctx,
1295 proj_context_is_network_enabled(ctx) ?
1296 PROJ_GRID_AVAILABILITY_KNOWN_AVAILABLE:
1297 PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
1298
1299 auto op_list = proj_create_operations(ctx, source_crs, target_crs, operation_ctx);
1300 proj_operation_factory_context_destroy(operation_ctx);
1301
1302 if( !op_list ) {
1303 return nullptr;
1304 }
1305
1306 auto op_count = proj_list_get_count(op_list);
1307 if( op_count == 0 ) {
1308 proj_list_destroy(op_list);
1309
1310 proj_context_log_debug(ctx, "No operation found matching criteria");
1311 return nullptr;
1312 }
1313
1314 PJ* P = proj_list_get(ctx, op_list, 0);
1315 assert(P);
1316
1317 if( P == nullptr || op_count == 1 || (area && area->bbox_set) ||
1318 proj_get_type(source_crs) == PJ_TYPE_GEOCENTRIC_CRS ||
1319 proj_get_type(target_crs) == PJ_TYPE_GEOCENTRIC_CRS ) {
1320 proj_list_destroy(op_list);
1321 return P;
1322 }
1323
1324 auto preparedOpList = pj_create_prepared_operations(ctx, source_crs, target_crs,
1325 op_list);
1326 proj_list_destroy(op_list);
1327
1328 if( preparedOpList.empty() )
1329 {
1330 proj_destroy(P);
1331 return nullptr;
1332 }
1333
1334 // If there's finally juste a single result, return it directly
1335 if( preparedOpList.size() == 1 )
1336 {
1337 auto retP = preparedOpList[0].pj;
1338 preparedOpList[0].pj = nullptr;
1339 proj_destroy(P);
1340 return retP;
1341 }
1342
1343 P->alternativeCoordinateOperations = std::move(preparedOpList);
1344 // The returned P is rather dummy
1345 P->iso_obj = nullptr;
1346 P->fwd = nullptr;
1347 P->inv = nullptr;
1348 P->fwd3d = nullptr;
1349 P->inv3d = nullptr;
1350 P->fwd4d = nullptr;
1351 P->inv4d = nullptr;
1352
1353 return P;
1354 }
1355
proj_destroy(PJ * P)1356 PJ *proj_destroy (PJ *P) {
1357 pj_free (P);
1358 return nullptr;
1359 }
1360
1361 /*****************************************************************************/
proj_errno(const PJ * P)1362 int proj_errno (const PJ *P) {
1363 /******************************************************************************
1364 Read an error level from the context of a PJ.
1365 ******************************************************************************/
1366 return pj_ctx_get_errno (pj_get_ctx ((PJ *) P));
1367 }
1368
1369 /*****************************************************************************/
proj_context_errno(PJ_CONTEXT * ctx)1370 int proj_context_errno (PJ_CONTEXT *ctx) {
1371 /******************************************************************************
1372 Read an error directly from a context, without going through a PJ
1373 belonging to that context.
1374 ******************************************************************************/
1375 if (nullptr==ctx)
1376 ctx = pj_get_default_ctx();
1377 return pj_ctx_get_errno (ctx);
1378 }
1379
1380 /*****************************************************************************/
proj_errno_set(const PJ * P,int err)1381 int proj_errno_set (const PJ *P, int err) {
1382 /******************************************************************************
1383 Set context-errno, bubble it up to the thread local errno, return err
1384 ******************************************************************************/
1385 /* Use proj_errno_reset to explicitly clear the error status */
1386 if (0==err)
1387 return 0;
1388
1389 /* For P==0 err goes to the default context */
1390 proj_context_errno_set (pj_get_ctx ((PJ *) P), err);
1391 errno = err;
1392 return err;
1393 }
1394
1395 /*****************************************************************************/
proj_errno_restore(const PJ * P,int err)1396 int proj_errno_restore (const PJ *P, int err) {
1397 /******************************************************************************
1398 Use proj_errno_restore when the current function succeeds, but the
1399 error flag was set on entry, and stored/reset using proj_errno_reset
1400 in order to monitor for new errors.
1401
1402 See usage example under proj_errno_reset ()
1403 ******************************************************************************/
1404 if (0==err)
1405 return 0;
1406 proj_errno_set (P, err);
1407 return 0;
1408 }
1409
1410 /*****************************************************************************/
proj_errno_reset(const PJ * P)1411 int proj_errno_reset (const PJ *P) {
1412 /******************************************************************************
1413 Clears errno in the context and thread local levels
1414 through the low level pj_ctx interface.
1415
1416 Returns the previous value of the errno, for convenient reset/restore
1417 operations:
1418
1419 int foo (PJ *P) {
1420 // errno may be set on entry, but we need to reset it to be able to
1421 // check for errors from "do_something_with_P(P)"
1422 int last_errno = proj_errno_reset (P);
1423
1424 // local failure
1425 if (0==P)
1426 return proj_errno_set (P, 42);
1427
1428 // call to function that may fail
1429 do_something_with_P (P);
1430
1431 // failure in do_something_with_P? - keep latest error status
1432 if (proj_errno(P))
1433 return proj_errno (P);
1434
1435 // success - restore previous error status, return 0
1436 return proj_errno_restore (P, last_errno);
1437 }
1438 ******************************************************************************/
1439 int last_errno;
1440 last_errno = proj_errno (P);
1441
1442 pj_ctx_set_errno (pj_get_ctx ((PJ *) P), 0);
1443 errno = 0;
1444 pj_errno = 0;
1445 return last_errno;
1446 }
1447
1448
1449 /* Create a new context based on the default context */
proj_context_create(void)1450 PJ_CONTEXT *proj_context_create (void) {
1451 return pj_ctx_alloc ();
1452 }
1453
1454
proj_context_destroy(PJ_CONTEXT * ctx)1455 PJ_CONTEXT *proj_context_destroy (PJ_CONTEXT *ctx) {
1456 if (nullptr==ctx)
1457 return nullptr;
1458
1459 /* Trying to free the default context is a no-op (since it is statically allocated) */
1460 if (pj_get_default_ctx ()==ctx)
1461 return nullptr;
1462
1463 pj_ctx_free (ctx);
1464 return nullptr;
1465 }
1466
1467
1468
1469
1470
1471
1472 /*****************************************************************************/
path_append(char * buf,const char * app,size_t * buf_size)1473 static char *path_append (char *buf, const char *app, size_t *buf_size) {
1474 /******************************************************************************
1475 Helper for proj_info() below. Append app to buf, separated by a
1476 semicolon. Also handle allocation of longer buffer if needed.
1477
1478 Returns buffer and adjusts *buf_size through provided pointer arg.
1479 ******************************************************************************/
1480 char *p;
1481 size_t len, applen = 0, buflen = 0;
1482 #ifdef _WIN32
1483 const char *delim = ";";
1484 #else
1485 const char *delim = ":";
1486 #endif
1487
1488 /* Nothing to do? */
1489 if (nullptr == app)
1490 return buf;
1491 applen = strlen (app);
1492 if (0 == applen)
1493 return buf;
1494
1495 /* Start checking whether buf is long enough */
1496 if (nullptr != buf)
1497 buflen = strlen (buf);
1498 len = buflen+applen+strlen (delim) + 1;
1499
1500 /* "pj_realloc", so to speak */
1501 if (*buf_size < len) {
1502 p = static_cast<char*>(pj_calloc (2 * len, sizeof (char)));
1503 if (nullptr==p) {
1504 pj_dealloc (buf);
1505 return nullptr;
1506 }
1507 *buf_size = 2 * len;
1508 if (buf != nullptr)
1509 strcpy (p, buf);
1510 pj_dealloc (buf);
1511 buf = p;
1512 }
1513 assert(buf);
1514
1515 /* Only append a delimiter if something's already there */
1516 if (0 != buflen)
1517 strcat (buf, delim);
1518 strcat (buf, app);
1519 return buf;
1520 }
1521
1522 static const char *empty = {""};
1523 static char version[64] = {""};
1524 static PJ_INFO info = {0, 0, 0, nullptr, nullptr, nullptr, nullptr, 0};
1525
1526 /*****************************************************************************/
proj_info(void)1527 PJ_INFO proj_info (void) {
1528 /******************************************************************************
1529 Basic info about the current instance of the PROJ.4 library.
1530
1531 Returns PJ_INFO struct.
1532 ******************************************************************************/
1533 size_t buf_size = 0;
1534 char *buf = nullptr;
1535
1536 pj_acquire_lock ();
1537
1538 info.major = PROJ_VERSION_MAJOR;
1539 info.minor = PROJ_VERSION_MINOR;
1540 info.patch = PROJ_VERSION_PATCH;
1541
1542 /* This is a controlled environment, so no risk of sprintf buffer
1543 overflow. A normal version string is xx.yy.zz which is 8 characters
1544 long and there is room for 64 bytes in the version string. */
1545 sprintf (version, "%d.%d.%d", info.major, info.minor, info.patch);
1546
1547 info.version = version;
1548 info.release = pj_get_release ();
1549
1550 /* build search path string */
1551 auto ctx = pj_get_default_ctx();
1552 if (!ctx || ctx->search_paths.empty()) {
1553 const auto searchpaths = pj_get_default_searchpaths(ctx);
1554 for( const auto& path: searchpaths ) {
1555 buf = path_append(buf, path.c_str(), &buf_size);
1556 }
1557 } else {
1558 for (const auto &path : ctx->search_paths) {
1559 buf = path_append(buf, path.c_str(), &buf_size);
1560 }
1561 }
1562
1563 pj_dalloc(const_cast<char*>(info.searchpath));
1564 info.searchpath = buf ? buf : empty;
1565
1566 info.paths = ctx ? ctx->c_compat_paths : nullptr;
1567 info.path_count = ctx ? static_cast<int>(ctx->search_paths.size()) : 0;
1568
1569 pj_release_lock ();
1570 return info;
1571 }
1572
1573
1574 /*****************************************************************************/
proj_pj_info(PJ * P)1575 PJ_PROJ_INFO proj_pj_info(PJ *P) {
1576 /******************************************************************************
1577 Basic info about a particular instance of a projection object.
1578
1579 Returns PJ_PROJ_INFO struct.
1580 ******************************************************************************/
1581 PJ_PROJ_INFO pjinfo;
1582 char *def;
1583
1584 memset(&pjinfo, 0, sizeof(PJ_PROJ_INFO));
1585
1586 pjinfo.accuracy = -1.0;
1587
1588 if (nullptr==P)
1589 return pjinfo;
1590
1591 /* coordinate operation description */
1592 if( P->iCurCoordOp >= 0 ) {
1593 P = P->alternativeCoordinateOperations[P->iCurCoordOp].pj;
1594 } else if( !P->alternativeCoordinateOperations.empty() ) {
1595 pjinfo.id = "unknown";
1596 pjinfo.description = "unavailable until proj_trans is called";
1597 pjinfo.definition = "unavailable until proj_trans is called";
1598 return pjinfo;
1599 }
1600
1601 /* projection id */
1602 if (pj_param(P->ctx, P->params, "tproj").i)
1603 pjinfo.id = pj_param(P->ctx, P->params, "sproj").s;
1604
1605 if( P->iso_obj ) {
1606 pjinfo.description = P->iso_obj->nameStr().c_str();
1607 } else {
1608 pjinfo.description = P->descr;
1609 }
1610
1611 // accuracy
1612 if( P->iso_obj ) {
1613 auto conv = dynamic_cast<const NS_PROJ::operation::Conversion*>(P->iso_obj.get());
1614 if( conv ) {
1615 pjinfo.accuracy = 0.0;
1616 } else {
1617 auto op = dynamic_cast<const NS_PROJ::operation::CoordinateOperation*>(P->iso_obj.get());
1618 if( op ) {
1619 const auto& accuracies = op->coordinateOperationAccuracies();
1620 if( !accuracies.empty() ) {
1621 try {
1622 pjinfo.accuracy = std::stod(accuracies[0]->value());
1623 } catch ( const std::exception& ) {}
1624 }
1625 }
1626 }
1627 }
1628
1629 /* projection definition */
1630 if (P->def_full)
1631 def = P->def_full;
1632 else
1633 def = pj_get_def(P, 0); /* pj_get_def takes a non-const PJ pointer */
1634 if (nullptr==def)
1635 pjinfo.definition = empty;
1636 else
1637 pjinfo.definition = pj_shrink (def);
1638 /* Make pj_free clean this up eventually */
1639 P->def_full = def;
1640
1641 pjinfo.has_inverse = pj_has_inverse(P);
1642 return pjinfo;
1643 }
1644
1645
1646 /*****************************************************************************/
proj_grid_info(const char * gridname)1647 PJ_GRID_INFO proj_grid_info(const char *gridname) {
1648 /******************************************************************************
1649 Information about a named datum grid.
1650
1651 Returns PJ_GRID_INFO struct.
1652 ******************************************************************************/
1653 PJ_GRID_INFO grinfo;
1654
1655 /*PJ_CONTEXT *ctx = proj_context_create(); */
1656 PJ_CONTEXT *ctx = pj_get_default_ctx();
1657 memset(&grinfo, 0, sizeof(PJ_GRID_INFO));
1658
1659 const auto fillGridInfo = [&grinfo, ctx, gridname]
1660 (const NS_PROJ::Grid& grid, const std::string& format)
1661 {
1662 const auto& extent = grid.extentAndRes();
1663
1664 /* name of grid */
1665 strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1);
1666
1667 /* full path of grid */
1668 pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1);
1669
1670 /* grid format */
1671 strncpy (grinfo.format, format.c_str(), sizeof(grinfo.format) - 1);
1672
1673 /* grid size */
1674 grinfo.n_lon = grid.width();
1675 grinfo.n_lat = grid.height();
1676
1677 /* cell size */
1678 grinfo.cs_lon = extent.resX;
1679 grinfo.cs_lat = extent.resY;
1680
1681 /* bounds of grid */
1682 grinfo.lowerleft.lam = extent.west;
1683 grinfo.lowerleft.phi = extent.south;
1684 grinfo.upperright.lam = extent.east;
1685 grinfo.upperright.phi = extent.north;
1686 };
1687
1688 {
1689 const auto gridSet = NS_PROJ::VerticalShiftGridSet::open(ctx, gridname);
1690 if( gridSet )
1691 {
1692 const auto& grids = gridSet->grids();
1693 if( !grids.empty() )
1694 {
1695 const auto& grid = grids.front();
1696 fillGridInfo(*grid, gridSet->format());
1697 return grinfo;
1698 }
1699 }
1700 }
1701
1702 {
1703 const auto gridSet = NS_PROJ::HorizontalShiftGridSet::open(ctx, gridname);
1704 if( gridSet )
1705 {
1706 const auto& grids = gridSet->grids();
1707 if( !grids.empty() )
1708 {
1709 const auto& grid = grids.front();
1710 fillGridInfo(*grid, gridSet->format());
1711 return grinfo;
1712 }
1713 }
1714 }
1715 strcpy(grinfo.format, "missing");
1716 return grinfo;
1717 }
1718
1719
1720
1721 /*****************************************************************************/
proj_init_info(const char * initname)1722 PJ_INIT_INFO proj_init_info(const char *initname){
1723 /******************************************************************************
1724 Information about a named init file.
1725
1726 Maximum length of initname is 64.
1727
1728 Returns PJ_INIT_INFO struct.
1729
1730 If the init file is not found all members of the return struct are set
1731 to the empty string.
1732
1733 If the init file is found, but the metadata is missing, the value is
1734 set to "Unknown".
1735 ******************************************************************************/
1736 int file_found;
1737 char param[80], key[74];
1738 paralist *start, *next;
1739 PJ_INIT_INFO ininfo;
1740 PJ_CONTEXT *ctx = pj_get_default_ctx();
1741
1742 memset(&ininfo, 0, sizeof(PJ_INIT_INFO));
1743
1744 file_found = pj_find_file(ctx, initname, ininfo.filename, sizeof(ininfo.filename));
1745 if (!file_found || strlen(initname) > 64) {
1746 if( strcmp(initname, "epsg") == 0 || strcmp(initname, "EPSG") == 0 ) {
1747 const char* val;
1748
1749 pj_ctx_set_errno( ctx, 0 );
1750
1751 strncpy (ininfo.name, initname, sizeof(ininfo.name) - 1);
1752 strcpy(ininfo.origin, "EPSG");
1753 val = proj_context_get_database_metadata(ctx, "EPSG.VERSION");
1754 if( val ) {
1755 strncpy(ininfo.version, val, sizeof(ininfo.version) - 1);
1756 }
1757 val = proj_context_get_database_metadata(ctx, "EPSG.DATE");
1758 if( val ) {
1759 strncpy(ininfo.lastupdate, val, sizeof(ininfo.lastupdate) - 1);
1760 }
1761 return ininfo;
1762 }
1763
1764 if( strcmp(initname, "IGNF") == 0 ) {
1765 const char* val;
1766
1767 pj_ctx_set_errno( ctx, 0 );
1768
1769 strncpy (ininfo.name, initname, sizeof(ininfo.name) - 1);
1770 strcpy(ininfo.origin, "IGNF");
1771 val = proj_context_get_database_metadata(ctx, "IGNF.VERSION");
1772 if( val ) {
1773 strncpy(ininfo.version, val, sizeof(ininfo.version) - 1);
1774 }
1775 val = proj_context_get_database_metadata(ctx, "IGNF.DATE");
1776 if( val ) {
1777 strncpy(ininfo.lastupdate, val, sizeof(ininfo.lastupdate) - 1);
1778 }
1779 return ininfo;
1780 }
1781
1782 return ininfo;
1783 }
1784
1785 /* The initial memset (0) makes strncpy safe here */
1786 strncpy (ininfo.name, initname, sizeof(ininfo.name) - 1);
1787 strcpy(ininfo.origin, "Unknown");
1788 strcpy(ininfo.version, "Unknown");
1789 strcpy(ininfo.lastupdate, "Unknown");
1790
1791 strncpy (key, initname, 64); /* make room for ":metadata\0" at the end */
1792 key[64] = 0;
1793 memcpy(key + strlen(key), ":metadata", 9 + 1);
1794 strcpy(param, "+init=");
1795 /* The +strlen(param) avoids a cppcheck false positive warning */
1796 strncat(param + strlen(param), key, sizeof(param)-1-strlen(param));
1797
1798 start = pj_mkparam(param);
1799 pj_expand_init(ctx, start);
1800
1801 if (pj_param(ctx, start, "tversion").i)
1802 strncpy(ininfo.version, pj_param(ctx, start, "sversion").s, sizeof(ininfo.version) - 1);
1803
1804 if (pj_param(ctx, start, "torigin").i)
1805 strncpy(ininfo.origin, pj_param(ctx, start, "sorigin").s, sizeof(ininfo.origin) - 1);
1806
1807 if (pj_param(ctx, start, "tlastupdate").i)
1808 strncpy(ininfo.lastupdate, pj_param(ctx, start, "slastupdate").s, sizeof(ininfo.lastupdate) - 1);
1809
1810 for ( ; start; start = next) {
1811 next = start->next;
1812 pj_dalloc(start);
1813 }
1814
1815 return ininfo;
1816 }
1817
1818
1819
1820 /*****************************************************************************/
proj_factors(PJ * P,PJ_COORD lp)1821 PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
1822 /******************************************************************************
1823 Cartographic characteristics at point lp.
1824
1825 Characteristics include meridian, parallel and areal scales, angular
1826 distortion, meridian/parallel, meridian convergence and scale error.
1827
1828 returns PJ_FACTORS. If unsuccessful, error number is set and the
1829 struct returned contains NULL data.
1830 ******************************************************************************/
1831 PJ_FACTORS factors = {0,0,0, 0,0,0, 0,0, 0,0,0,0};
1832 struct FACTORS f;
1833
1834 if (nullptr==P)
1835 return factors;
1836
1837 if (pj_factors(lp.lp, P, 0.0, &f))
1838 return factors;
1839
1840 factors.meridional_scale = f.h;
1841 factors.parallel_scale = f.k;
1842 factors.areal_scale = f.s;
1843
1844 factors.angular_distortion = f.omega;
1845 factors.meridian_parallel_angle = f.thetap;
1846 factors.meridian_convergence = f.conv;
1847
1848 factors.tissot_semimajor = f.a;
1849 factors.tissot_semiminor = f.b;
1850
1851 /* Raw derivatives, for completeness's sake */
1852 factors.dx_dlam = f.der.x_l;
1853 factors.dx_dphi = f.der.x_p;
1854 factors.dy_dlam = f.der.y_l;
1855 factors.dy_dphi = f.der.y_p;
1856
1857 return factors;
1858 }
1859
1860