1 : /******************************************************************************
2 : * $Id: ogrct.cpp 25007 2012-09-30 08:49:33Z rouault $
3 : *
4 : * Project: OpenGIS Simple Features Reference Implementation
5 : * Purpose: The OGRSCoordinateTransformation class.
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 "ogr_spatialref.h"
31 : #include "cpl_port.h"
32 : #include "cpl_error.h"
33 : #include "cpl_conv.h"
34 : #include "cpl_string.h"
35 : #include "cpl_multiproc.h"
36 :
37 : #ifdef PROJ_STATIC
38 : #include "proj_api.h"
39 : #endif
40 :
41 : CPL_CVSID("$Id: ogrct.cpp 25007 2012-09-30 08:49:33Z rouault $");
42 :
43 : /* ==================================================================== */
44 : /* PROJ.4 interface stuff. */
45 : /* ==================================================================== */
46 : #ifndef PROJ_STATIC
47 : typedef struct { double u, v; } projUV;
48 :
49 : #define projPJ void *
50 : #define projCtx void *
51 : #define RAD_TO_DEG 57.29577951308232
52 : #define DEG_TO_RAD .0174532925199432958
53 :
54 : #else
55 :
56 : #if PJ_VERSION < 480
57 : #define projCtx void *
58 : #endif
59 :
60 : #endif
61 :
62 : static void *hPROJMutex = NULL;
63 :
64 : static projPJ (*pfn_pj_init_plus)(const char *) = NULL;
65 : static projPJ (*pfn_pj_init)(int, char**) = NULL;
66 : static void (*pfn_pj_free)(projPJ) = NULL;
67 : static int (*pfn_pj_transform)(projPJ, projPJ, long, int,
68 : double *, double *, double * ) = NULL;
69 : static int *(*pfn_pj_get_errno_ref)(void) = NULL;
70 : static char *(*pfn_pj_strerrno)(int) = NULL;
71 : static char *(*pfn_pj_get_def)(projPJ,int) = NULL;
72 : static void (*pfn_pj_dalloc)(void *) = NULL;
73 :
74 : static projPJ (*pfn_pj_init_plus_ctx)( projCtx, const char * ) = NULL;
75 : static int (*pfn_pj_ctx_get_errno)( projCtx ) = NULL;
76 : static projCtx (*pfn_pj_ctx_alloc)(void) = NULL;
77 : static void (*pfn_pj_ctx_free)( projCtx ) = NULL;
78 :
79 : #if (defined(WIN32) || defined(WIN32CE)) && !defined(__MINGW32__)
80 : # define LIBNAME "proj.dll"
81 : #elif defined(__MINGW32__)
82 : // XXX: If PROJ.4 library was properly built using libtool in Cygwin or MinGW
83 : // environments it has the interface version number embedded in the file name
84 : // (it is CURRENT-AGE number). If DLL came somewhere else (e.g. from MSVC
85 : // build) it can be named either way, so use PROJSO environment variable to
86 : // specify the right library name. By default assume that in Cygwin/MinGW all
87 : // components were buit in the same way.
88 : # define LIBNAME "libproj-0.dll"
89 : #elif defined(__CYGWIN__)
90 : # define LIBNAME "cygproj-0.dll"
91 : #elif defined(__APPLE__)
92 : # define LIBNAME "libproj.dylib"
93 : #else
94 : # define LIBNAME "libproj.so"
95 : #endif
96 :
97 : /************************************************************************/
98 : /* OGRProj4CT */
99 : /************************************************************************/
100 :
101 : class OGRProj4CT : public OGRCoordinateTransformation
102 : {
103 : OGRSpatialReference *poSRSSource;
104 : void *psPJSource;
105 : int bSourceLatLong;
106 : double dfSourceToRadians;
107 : double dfSourceFromRadians;
108 : int bSourceWrap;
109 : double dfSourceWrapLong;
110 :
111 :
112 : OGRSpatialReference *poSRSTarget;
113 : void *psPJTarget;
114 : int bTargetLatLong;
115 : double dfTargetToRadians;
116 : double dfTargetFromRadians;
117 : int bTargetWrap;
118 : double dfTargetWrapLong;
119 :
120 : int nErrorCount;
121 :
122 : int bCheckWithInvertProj;
123 : double dfThreshold;
124 :
125 : projCtx pjctx;
126 :
127 : int InitializeNoLock( OGRSpatialReference *poSource,
128 : OGRSpatialReference *poTarget );
129 :
130 : int nMaxCount;
131 : double *padfOriX;
132 : double *padfOriY;
133 : double *padfOriZ;
134 : double *padfTargetX;
135 : double *padfTargetY;
136 : double *padfTargetZ;
137 :
138 : public:
139 : OGRProj4CT();
140 : virtual ~OGRProj4CT();
141 :
142 : int Initialize( OGRSpatialReference *poSource,
143 : OGRSpatialReference *poTarget );
144 :
145 : virtual OGRSpatialReference *GetSourceCS();
146 : virtual OGRSpatialReference *GetTargetCS();
147 : virtual int Transform( int nCount,
148 : double *x, double *y, double *z = NULL );
149 : virtual int TransformEx( int nCount,
150 : double *x, double *y, double *z = NULL,
151 : int *panSuccess = NULL );
152 :
153 : };
154 :
155 : /************************************************************************/
156 : /* GetProjLibraryName() */
157 : /************************************************************************/
158 :
159 39 : static const char* GetProjLibraryName()
160 : {
161 39 : const char *pszLibName = LIBNAME;
162 : #if !defined(WIN32CE)
163 39 : if( CPLGetConfigOption("PROJSO",NULL) != NULL )
164 0 : pszLibName = CPLGetConfigOption("PROJSO",NULL);
165 : #endif
166 39 : return pszLibName;
167 : }
168 :
169 : /************************************************************************/
170 : /* LoadProjLibrary() */
171 : /************************************************************************/
172 :
173 598 : static int LoadProjLibrary()
174 :
175 : {
176 598 : CPLMutexHolderD( &hPROJMutex );
177 : static int bTriedToLoad = FALSE;
178 : const char *pszLibName;
179 :
180 598 : if( bTriedToLoad )
181 559 : return( pfn_pj_transform != NULL );
182 :
183 39 : bTriedToLoad = TRUE;
184 :
185 39 : pszLibName = GetProjLibraryName();
186 :
187 : #ifdef PROJ_STATIC
188 : pfn_pj_init = pj_init;
189 : pfn_pj_init_plus = pj_init_plus;
190 : pfn_pj_free = pj_free;
191 : pfn_pj_transform = pj_transform;
192 : pfn_pj_get_errno_ref = (int *(*)(void)) pj_get_errno_ref;
193 : pfn_pj_strerrno = pj_strerrno;
194 : pfn_pj_dalloc = pj_dalloc;
195 : #if PJ_VERSION >= 446
196 : pfn_pj_get_def = pj_get_def;
197 : #endif
198 : #if PJ_VERSION >= 480
199 : pfn_pj_ctx_alloc = pj_ctx_alloc;
200 : pfn_pj_ctx_free = pj_ctx_free;
201 : pfn_pj_init_plus_ctx = pj_init_plus_ctx;
202 : pfn_pj_ctx_get_errno = pj_ctx_get_errno;
203 : #endif
204 : #else
205 39 : CPLPushErrorHandler( CPLQuietErrorHandler );
206 :
207 : pfn_pj_init = (projPJ (*)(int, char**)) CPLGetSymbol( pszLibName,
208 39 : "pj_init" );
209 39 : CPLPopErrorHandler();
210 :
211 39 : if( pfn_pj_init == NULL )
212 0 : return( FALSE );
213 :
214 : pfn_pj_init_plus = (projPJ (*)(const char *))
215 39 : CPLGetSymbol( pszLibName, "pj_init_plus" );
216 : pfn_pj_free = (void (*)(projPJ))
217 39 : CPLGetSymbol( pszLibName, "pj_free" );
218 : pfn_pj_transform = (int (*)(projPJ,projPJ,long,int,double*,
219 : double*,double*))
220 39 : CPLGetSymbol( pszLibName, "pj_transform" );
221 : pfn_pj_get_errno_ref = (int *(*)(void))
222 39 : CPLGetSymbol( pszLibName, "pj_get_errno_ref" );
223 : pfn_pj_strerrno = (char *(*)(int))
224 39 : CPLGetSymbol( pszLibName, "pj_strerrno" );
225 :
226 39 : CPLPushErrorHandler( CPLQuietErrorHandler );
227 : pfn_pj_get_def = (char *(*)(projPJ,int))
228 39 : CPLGetSymbol( pszLibName, "pj_get_def" );
229 : pfn_pj_dalloc = (void (*)(void*))
230 39 : CPLGetSymbol( pszLibName, "pj_dalloc" );
231 :
232 : /* PROJ 4.8.0 symbols */
233 : pfn_pj_ctx_alloc = (projCtx (*)( void ))
234 39 : CPLGetSymbol( pszLibName, "pj_ctx_alloc" );
235 : pfn_pj_ctx_free = (void (*)( projCtx ))
236 39 : CPLGetSymbol( pszLibName, "pj_ctx_free" );
237 : pfn_pj_init_plus_ctx = (projPJ (*)( projCtx, const char * ))
238 39 : CPLGetSymbol( pszLibName, "pj_init_plus_ctx" );
239 : pfn_pj_ctx_get_errno = (int (*)( projCtx ))
240 39 : CPLGetSymbol( pszLibName, "pj_ctx_get_errno" );
241 :
242 39 : CPLPopErrorHandler();
243 39 : CPLErrorReset();
244 : #endif
245 :
246 39 : if (pfn_pj_ctx_alloc != NULL &&
247 : pfn_pj_ctx_free != NULL &&
248 : pfn_pj_init_plus_ctx != NULL &&
249 : pfn_pj_ctx_get_errno != NULL &&
250 : CSLTestBoolean(CPLGetConfigOption("USE_PROJ_480_FEATURES", "YES")))
251 : {
252 39 : CPLDebug("OGRCT", "PROJ >= 4.8.0 features enabled");
253 : }
254 : else
255 : {
256 0 : pfn_pj_ctx_alloc = NULL;
257 0 : pfn_pj_ctx_free = NULL;
258 0 : pfn_pj_init_plus_ctx = NULL;
259 0 : pfn_pj_ctx_get_errno = NULL;
260 : }
261 :
262 39 : if( pfn_pj_transform == NULL )
263 : {
264 : CPLError( CE_Failure, CPLE_AppDefined,
265 : "Attempt to load %s, but couldn't find pj_transform.\n"
266 : "Please upgrade to PROJ 4.1.2 or later.",
267 0 : pszLibName );
268 :
269 0 : return FALSE;
270 : }
271 :
272 39 : return( TRUE );
273 : }
274 :
275 : /************************************************************************/
276 : /* OCTProj4Normalize() */
277 : /* */
278 : /* This function is really just here since we already have all */
279 : /* the code to load libproj.so. It is intended to "normalize" */
280 : /* a proj.4 definition, expanding +init= definitions and so */
281 : /* forth as possible. */
282 : /************************************************************************/
283 :
284 567 : char *OCTProj4Normalize( const char *pszProj4Src )
285 :
286 : {
287 : char *pszNewProj4Def, *pszCopy;
288 567 : projPJ psPJSource = NULL;
289 :
290 567 : CPLMutexHolderD( &hPROJMutex );
291 :
292 567 : if( !LoadProjLibrary() || pfn_pj_dalloc == NULL || pfn_pj_get_def == NULL )
293 0 : return CPLStrdup( pszProj4Src );
294 :
295 567 : psPJSource = pfn_pj_init_plus( pszProj4Src );
296 :
297 567 : if( psPJSource == NULL )
298 438 : return CPLStrdup( pszProj4Src );
299 :
300 129 : pszNewProj4Def = pfn_pj_get_def( psPJSource, 0 );
301 :
302 129 : pfn_pj_free( psPJSource );
303 :
304 129 : if( pszNewProj4Def == NULL )
305 0 : return CPLStrdup( pszProj4Src );
306 :
307 129 : pszCopy = CPLStrdup( pszNewProj4Def );
308 129 : pfn_pj_dalloc( pszNewProj4Def );
309 :
310 129 : return pszCopy;
311 : }
312 :
313 : /************************************************************************/
314 : /* OCTDestroyCoordinateTransformation() */
315 : /************************************************************************/
316 :
317 : /**
318 : * \brief OGRCoordinateTransformation destructor.
319 : *
320 : * This function is the same as OGRCoordinateTransformation::DestroyCT()
321 : *
322 : * @param hCT the object to delete
323 : */
324 :
325 : void CPL_STDCALL
326 104 : OCTDestroyCoordinateTransformation( OGRCoordinateTransformationH hCT )
327 :
328 : {
329 104 : delete (OGRCoordinateTransformation *) hCT;
330 104 : }
331 :
332 : /************************************************************************/
333 : /* DestroyCT() */
334 : /************************************************************************/
335 :
336 : /**
337 : * \brief OGRCoordinateTransformation destructor.
338 : *
339 : * This function is the same as OGRCoordinateTransformation::~OGRCoordinateTransformation()
340 : * and OCTDestroyCoordinateTransformation()
341 : *
342 : * This static method will destroy a OGRCoordinateTransformation. It is
343 : * equivalent to calling delete on the object, but it ensures that the
344 : * deallocation is properly executed within the OGR libraries heap on
345 : * platforms where this can matter (win32).
346 : *
347 : * @param poCT the object to delete
348 : *
349 : * @since GDAL 1.7.0
350 : */
351 :
352 19 : void OGRCoordinateTransformation::DestroyCT(OGRCoordinateTransformation* poCT)
353 : {
354 19 : delete poCT;
355 19 : }
356 :
357 : /************************************************************************/
358 : /* OGRCreateCoordinateTransformation() */
359 : /************************************************************************/
360 :
361 : /**
362 : * Create transformation object.
363 : *
364 : * This is the same as the C function OCTNewCoordinateTransformation().
365 : *
366 : * Input spatial reference system objects are assigned
367 : * by copy (calling clone() method) and no ownership transfer occurs.
368 : *
369 : * The delete operator, or OCTDestroyCoordinateTransformation() should
370 : * be used to destroy transformation objects.
371 : *
372 : * The PROJ.4 library must be available at run-time.
373 : *
374 : * @param poSource source spatial reference system.
375 : * @param poTarget target spatial reference system.
376 : * @return NULL on failure or a ready to use transformation object.
377 : */
378 :
379 : OGRCoordinateTransformation*
380 478 : OGRCreateCoordinateTransformation( OGRSpatialReference *poSource,
381 : OGRSpatialReference *poTarget )
382 :
383 : {
384 : OGRProj4CT *poCT;
385 :
386 478 : if( pfn_pj_init == NULL && !LoadProjLibrary() )
387 : {
388 : CPLError( CE_Failure, CPLE_NotSupported,
389 : "Unable to load PROJ.4 library (%s), creation of\n"
390 : "OGRCoordinateTransformation failed.",
391 0 : GetProjLibraryName() );
392 0 : return NULL;
393 : }
394 :
395 478 : poCT = new OGRProj4CT();
396 :
397 478 : if( !poCT->Initialize( poSource, poTarget ) )
398 : {
399 1 : delete poCT;
400 1 : return NULL;
401 : }
402 : else
403 : {
404 477 : return poCT;
405 : }
406 : }
407 :
408 : /************************************************************************/
409 : /* OCTNewCoordinateTransformation() */
410 : /************************************************************************/
411 :
412 : /**
413 : * Create transformation object.
414 : *
415 : * This is the same as the C++ function OGRCreateCoordinateTransformation().
416 : *
417 : * Input spatial reference system objects are assigned
418 : * by copy (calling clone() method) and no ownership transfer occurs.
419 : *
420 : * OCTDestroyCoordinateTransformation() should
421 : * be used to destroy transformation objects.
422 : *
423 : * The PROJ.4 library must be available at run-time.
424 : *
425 : * @param hSourceSRS source spatial reference system.
426 : * @param hTargetSRS target spatial reference system.
427 : * @return NULL on failure or a ready to use transformation object.
428 : */
429 :
430 : OGRCoordinateTransformationH CPL_STDCALL
431 102 : OCTNewCoordinateTransformation(
432 : OGRSpatialReferenceH hSourceSRS, OGRSpatialReferenceH hTargetSRS )
433 :
434 : {
435 : return (OGRCoordinateTransformationH)
436 : OGRCreateCoordinateTransformation(
437 : (OGRSpatialReference *) hSourceSRS,
438 102 : (OGRSpatialReference *) hTargetSRS );
439 : }
440 :
441 : /************************************************************************/
442 : /* OGRProj4CT() */
443 : /************************************************************************/
444 :
445 478 : OGRProj4CT::OGRProj4CT()
446 :
447 : {
448 478 : poSRSSource = NULL;
449 478 : poSRSTarget = NULL;
450 478 : psPJSource = NULL;
451 478 : psPJTarget = NULL;
452 :
453 478 : nErrorCount = 0;
454 :
455 478 : bCheckWithInvertProj = FALSE;
456 478 : dfThreshold = 0;
457 :
458 478 : nMaxCount = 0;
459 478 : padfOriX = NULL;
460 478 : padfOriY = NULL;
461 478 : padfOriZ = NULL;
462 478 : padfTargetX = NULL;
463 478 : padfTargetY = NULL;
464 478 : padfTargetZ = NULL;
465 :
466 478 : if (pfn_pj_ctx_alloc != NULL)
467 478 : pjctx = pfn_pj_ctx_alloc();
468 : else
469 0 : pjctx = NULL;
470 478 : }
471 :
472 : /************************************************************************/
473 : /* ~OGRProj4CT() */
474 : /************************************************************************/
475 :
476 478 : OGRProj4CT::~OGRProj4CT()
477 :
478 : {
479 478 : if( poSRSSource != NULL )
480 : {
481 477 : if( poSRSSource->Dereference() <= 0 )
482 477 : delete poSRSSource;
483 : }
484 :
485 478 : if( poSRSTarget != NULL )
486 : {
487 477 : if( poSRSTarget->Dereference() <= 0 )
488 474 : delete poSRSTarget;
489 : }
490 :
491 478 : if (pjctx != NULL)
492 : {
493 478 : pfn_pj_ctx_free(pjctx);
494 :
495 478 : if( psPJSource != NULL )
496 477 : pfn_pj_free( psPJSource );
497 :
498 478 : if( psPJTarget != NULL )
499 477 : pfn_pj_free( psPJTarget );
500 : }
501 : else
502 : {
503 0 : CPLMutexHolderD( &hPROJMutex );
504 :
505 0 : if( psPJSource != NULL )
506 0 : pfn_pj_free( psPJSource );
507 :
508 0 : if( psPJTarget != NULL )
509 0 : pfn_pj_free( psPJTarget );
510 : }
511 :
512 478 : CPLFree(padfOriX);
513 478 : CPLFree(padfOriY);
514 478 : CPLFree(padfOriZ);
515 478 : CPLFree(padfTargetX);
516 478 : CPLFree(padfTargetY);
517 478 : CPLFree(padfTargetZ);
518 478 : }
519 :
520 : /************************************************************************/
521 : /* Initialize() */
522 : /************************************************************************/
523 :
524 478 : int OGRProj4CT::Initialize( OGRSpatialReference * poSourceIn,
525 : OGRSpatialReference * poTargetIn )
526 :
527 : {
528 478 : if (pjctx != NULL)
529 : {
530 478 : return InitializeNoLock(poSourceIn, poTargetIn);
531 : }
532 :
533 0 : CPLMutexHolderD( &hPROJMutex );
534 0 : return InitializeNoLock(poSourceIn, poTargetIn);
535 : }
536 :
537 : /************************************************************************/
538 : /* InitializeNoLock() */
539 : /************************************************************************/
540 :
541 478 : int OGRProj4CT::InitializeNoLock( OGRSpatialReference * poSourceIn,
542 : OGRSpatialReference * poTargetIn )
543 :
544 : {
545 478 : if( poSourceIn == NULL || poTargetIn == NULL )
546 1 : return FALSE;
547 :
548 477 : poSRSSource = poSourceIn->Clone();
549 477 : poSRSTarget = poTargetIn->Clone();
550 :
551 477 : bSourceLatLong = poSRSSource->IsGeographic();
552 477 : bTargetLatLong = poSRSTarget->IsGeographic();
553 :
554 : /* -------------------------------------------------------------------- */
555 : /* Setup source and target translations to radians for lat/long */
556 : /* systems. */
557 : /* -------------------------------------------------------------------- */
558 477 : dfSourceToRadians = DEG_TO_RAD;
559 477 : dfSourceFromRadians = RAD_TO_DEG;
560 477 : bSourceWrap = FALSE;
561 477 : dfSourceWrapLong = 0.0;
562 :
563 477 : if( bSourceLatLong )
564 : {
565 277 : OGR_SRSNode *poUNITS = poSRSSource->GetAttrNode( "GEOGCS|UNIT" );
566 277 : if( poUNITS && poUNITS->GetChildCount() >= 2 )
567 : {
568 277 : dfSourceToRadians = atof(poUNITS->GetChild(1)->GetValue());
569 277 : if( dfSourceToRadians == 0.0 )
570 0 : dfSourceToRadians = DEG_TO_RAD;
571 : else
572 277 : dfSourceFromRadians = 1 / dfSourceToRadians;
573 : }
574 : }
575 :
576 477 : dfTargetToRadians = DEG_TO_RAD;
577 477 : dfTargetFromRadians = RAD_TO_DEG;
578 477 : bTargetWrap = FALSE;
579 477 : dfTargetWrapLong = 0.0;
580 :
581 477 : if( bTargetLatLong )
582 : {
583 202 : OGR_SRSNode *poUNITS = poSRSTarget->GetAttrNode( "GEOGCS|UNIT" );
584 202 : if( poUNITS && poUNITS->GetChildCount() >= 2 )
585 : {
586 202 : dfTargetToRadians = atof(poUNITS->GetChild(1)->GetValue());
587 202 : if( dfTargetToRadians == 0.0 )
588 0 : dfTargetToRadians = DEG_TO_RAD;
589 : else
590 202 : dfTargetFromRadians = 1 / dfTargetToRadians;
591 : }
592 : }
593 :
594 : /* -------------------------------------------------------------------- */
595 : /* Preliminary logic to setup wrapping. */
596 : /* -------------------------------------------------------------------- */
597 : const char *pszCENTER_LONG;
598 :
599 477 : if( CPLGetConfigOption( "CENTER_LONG", NULL ) != NULL )
600 : {
601 0 : bSourceWrap = bTargetWrap = TRUE;
602 : dfSourceWrapLong = dfTargetWrapLong =
603 0 : atof(CPLGetConfigOption( "CENTER_LONG", "" ));
604 0 : CPLDebug( "OGRCT", "Wrap at %g.", dfSourceWrapLong );
605 : }
606 :
607 477 : pszCENTER_LONG = poSRSSource->GetExtension( "GEOGCS", "CENTER_LONG" );
608 477 : if( pszCENTER_LONG != NULL )
609 : {
610 36 : dfSourceWrapLong = atof(pszCENTER_LONG);
611 36 : bSourceWrap = TRUE;
612 36 : CPLDebug( "OGRCT", "Wrap source at %g.", dfSourceWrapLong );
613 : }
614 :
615 477 : pszCENTER_LONG = poSRSTarget->GetExtension( "GEOGCS", "CENTER_LONG" );
616 477 : if( pszCENTER_LONG != NULL )
617 : {
618 36 : dfTargetWrapLong = atof(pszCENTER_LONG);
619 36 : bTargetWrap = TRUE;
620 36 : CPLDebug( "OGRCT", "Wrap target at %g.", dfTargetWrapLong );
621 : }
622 :
623 477 : bCheckWithInvertProj = CSLTestBoolean(CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", "NO" ));
624 :
625 : /* The threshold is rather experimental... Works well with the cases of ticket #2305 */
626 477 : if (bSourceLatLong)
627 277 : dfThreshold = atof(CPLGetConfigOption( "THRESHOLD", ".1" ));
628 : else
629 : /* 1 works well for most projections, except for +proj=aeqd that requires */
630 : /* a tolerance of 10000 */
631 200 : dfThreshold = atof(CPLGetConfigOption( "THRESHOLD", "10000" ));
632 :
633 : /* -------------------------------------------------------------------- */
634 : /* Establish PROJ.4 handle for source if projection. */
635 : /* -------------------------------------------------------------------- */
636 : // OGRThreadSafety: The following variable is not a thread safety issue
637 : // since the only issue is incrementing while accessing which at worse
638 : // means debug output could be one "increment" late.
639 : static int nDebugReportCount = 0;
640 :
641 477 : char *pszProj4Defn = NULL;
642 :
643 477 : if( poSRSSource->exportToProj4( &pszProj4Defn ) != OGRERR_NONE )
644 : {
645 0 : CPLFree( pszProj4Defn );
646 0 : return FALSE;
647 : }
648 :
649 477 : if( strlen(pszProj4Defn) == 0 )
650 : {
651 0 : CPLFree( pszProj4Defn );
652 : CPLError( CE_Failure, CPLE_AppDefined,
653 : "No PROJ.4 translation for source SRS, coordinate\n"
654 0 : "transformation initialization has failed." );
655 0 : return FALSE;
656 : }
657 :
658 477 : if (pjctx)
659 477 : psPJSource = pfn_pj_init_plus_ctx( pjctx, pszProj4Defn );
660 : else
661 0 : psPJSource = pfn_pj_init_plus( pszProj4Defn );
662 :
663 477 : if( psPJSource == NULL )
664 : {
665 0 : if( pjctx != NULL)
666 : {
667 0 : int pj_errno = pfn_pj_ctx_get_errno(pjctx);
668 :
669 : /* pfn_pj_strerrno not yet thread-safe in PROJ 4.8.0 */
670 0 : CPLMutexHolderD(&hPROJMutex);
671 : CPLError( CE_Failure, CPLE_NotSupported,
672 : "Failed to initialize PROJ.4 with `%s'.\n%s",
673 0 : pszProj4Defn, pfn_pj_strerrno(pj_errno) );
674 : }
675 0 : else if( pfn_pj_get_errno_ref != NULL
676 : && pfn_pj_strerrno != NULL )
677 : {
678 0 : int *p_pj_errno = pfn_pj_get_errno_ref();
679 :
680 : CPLError( CE_Failure, CPLE_NotSupported,
681 : "Failed to initialize PROJ.4 with `%s'.\n%s",
682 0 : pszProj4Defn, pfn_pj_strerrno(*p_pj_errno) );
683 : }
684 : else
685 : {
686 : CPLError( CE_Failure, CPLE_NotSupported,
687 : "Failed to initialize PROJ.4 with `%s'.\n",
688 0 : pszProj4Defn );
689 : }
690 : }
691 :
692 477 : if( nDebugReportCount < 10 )
693 95 : CPLDebug( "OGRCT", "Source: %s", pszProj4Defn );
694 :
695 477 : CPLFree( pszProj4Defn );
696 :
697 477 : if( psPJSource == NULL )
698 0 : return FALSE;
699 :
700 : /* -------------------------------------------------------------------- */
701 : /* Establish PROJ.4 handle for target if projection. */
702 : /* -------------------------------------------------------------------- */
703 477 : pszProj4Defn = NULL;
704 :
705 477 : if( poSRSTarget->exportToProj4( &pszProj4Defn ) != OGRERR_NONE )
706 : {
707 0 : CPLFree( pszProj4Defn );
708 0 : return FALSE;
709 : }
710 :
711 477 : if( strlen(pszProj4Defn) == 0 )
712 : {
713 0 : CPLFree( pszProj4Defn );
714 : CPLError( CE_Failure, CPLE_AppDefined,
715 : "No PROJ.4 translation for destination SRS, coordinate\n"
716 0 : "transformation initialization has failed." );
717 0 : return FALSE;
718 : }
719 :
720 477 : if (pjctx)
721 477 : psPJTarget = pfn_pj_init_plus_ctx( pjctx, pszProj4Defn );
722 : else
723 0 : psPJTarget = pfn_pj_init_plus( pszProj4Defn );
724 :
725 477 : if( psPJTarget == NULL )
726 : CPLError( CE_Failure, CPLE_NotSupported,
727 : "Failed to initialize PROJ.4 with `%s'.",
728 0 : pszProj4Defn );
729 :
730 477 : if( nDebugReportCount < 10 )
731 : {
732 95 : CPLDebug( "OGRCT", "Target: %s", pszProj4Defn );
733 95 : nDebugReportCount++;
734 : }
735 :
736 477 : CPLFree( pszProj4Defn );
737 :
738 477 : if( psPJTarget == NULL )
739 0 : return FALSE;
740 :
741 477 : return TRUE;
742 : }
743 :
744 : /************************************************************************/
745 : /* GetSourceCS() */
746 : /************************************************************************/
747 :
748 38 : OGRSpatialReference *OGRProj4CT::GetSourceCS()
749 :
750 : {
751 38 : return poSRSSource;
752 : }
753 :
754 : /************************************************************************/
755 : /* GetTargetCS() */
756 : /************************************************************************/
757 :
758 372 : OGRSpatialReference *OGRProj4CT::GetTargetCS()
759 :
760 : {
761 372 : return poSRSTarget;
762 : }
763 :
764 : /************************************************************************/
765 : /* Transform() */
766 : /* */
767 : /* This is a small wrapper for the extended transform version. */
768 : /************************************************************************/
769 :
770 1247 : int OGRProj4CT::Transform( int nCount, double *x, double *y, double *z )
771 :
772 : {
773 1247 : int *pabSuccess = (int *) CPLMalloc(sizeof(int) * nCount );
774 : int bOverallSuccess, i;
775 :
776 1247 : bOverallSuccess = TransformEx( nCount, x, y, z, pabSuccess );
777 :
778 2501 : for( i = 0; i < nCount; i++ )
779 : {
780 1254 : if( !pabSuccess[i] )
781 : {
782 0 : bOverallSuccess = FALSE;
783 0 : break;
784 : }
785 : }
786 :
787 1247 : CPLFree( pabSuccess );
788 :
789 1247 : return bOverallSuccess;
790 : }
791 :
792 : /************************************************************************/
793 : /* OCTTransform() */
794 : /************************************************************************/
795 :
796 344 : int CPL_STDCALL OCTTransform( OGRCoordinateTransformationH hTransform,
797 : int nCount, double *x, double *y, double *z )
798 :
799 : {
800 344 : VALIDATE_POINTER1( hTransform, "OCTTransform", FALSE );
801 :
802 : return ((OGRCoordinateTransformation*) hTransform)->
803 344 : Transform( nCount, x, y,z );
804 : }
805 :
806 : /************************************************************************/
807 : /* TransformEx() */
808 : /************************************************************************/
809 :
810 39098 : int OGRProj4CT::TransformEx( int nCount, double *x, double *y, double *z,
811 : int *pabSuccess )
812 :
813 : {
814 : int err, i;
815 :
816 : /* -------------------------------------------------------------------- */
817 : /* Potentially transform to radians. */
818 : /* -------------------------------------------------------------------- */
819 39098 : if( bSourceLatLong )
820 : {
821 10672 : if( bSourceWrap )
822 : {
823 14826 : for( i = 0; i < nCount; i++ )
824 : {
825 10511 : if( x[i] != HUGE_VAL && y[i] != HUGE_VAL )
826 : {
827 9854 : if( x[i] < dfSourceWrapLong - 180.0 )
828 0 : x[i] += 360.0;
829 9854 : else if( x[i] > dfSourceWrapLong + 180 )
830 0 : x[i] -= 360.0;
831 : }
832 : }
833 : }
834 :
835 3797789 : for( i = 0; i < nCount; i++ )
836 : {
837 3787117 : if( x[i] != HUGE_VAL )
838 : {
839 3786460 : x[i] *= dfSourceToRadians;
840 3786460 : y[i] *= dfSourceToRadians;
841 : }
842 : }
843 : }
844 :
845 : /* -------------------------------------------------------------------- */
846 : /* Do the transformation using PROJ.4. */
847 : /* -------------------------------------------------------------------- */
848 39098 : if (pjctx == NULL)
849 : {
850 : /* The mutex has already been created */
851 0 : CPLAssert(hPROJMutex != NULL);
852 0 : CPLAcquireMutex(hPROJMutex, 1000.0);
853 : }
854 :
855 39098 : if (bCheckWithInvertProj)
856 : {
857 : /* For some projections, we cannot detect if we are trying to reproject */
858 : /* coordinates outside the validity area of the projection. So let's do */
859 : /* the reverse reprojection and compare with the source coordinates */
860 5617 : if (nCount > nMaxCount)
861 : {
862 21 : nMaxCount = nCount;
863 21 : padfOriX = (double*) CPLRealloc(padfOriX, sizeof(double)*nCount);
864 21 : padfOriY = (double*) CPLRealloc(padfOriY, sizeof(double)*nCount);
865 21 : padfOriZ = (double*) CPLRealloc(padfOriZ, sizeof(double)*nCount);
866 21 : padfTargetX = (double*) CPLRealloc(padfTargetX, sizeof(double)*nCount);
867 21 : padfTargetY = (double*) CPLRealloc(padfTargetY, sizeof(double)*nCount);
868 21 : padfTargetZ = (double*) CPLRealloc(padfTargetZ, sizeof(double)*nCount);
869 : }
870 5617 : memcpy(padfOriX, x, sizeof(double)*nCount);
871 5617 : memcpy(padfOriY, y, sizeof(double)*nCount);
872 5617 : if (z)
873 : {
874 5617 : memcpy(padfOriZ, z, sizeof(double)*nCount);
875 : }
876 5617 : err = pfn_pj_transform( psPJSource, psPJTarget, nCount, 1, x, y, z );
877 5617 : if (err == 0)
878 : {
879 5617 : memcpy(padfTargetX, x, sizeof(double)*nCount);
880 5617 : memcpy(padfTargetY, y, sizeof(double)*nCount);
881 5617 : if (z)
882 : {
883 5617 : memcpy(padfTargetZ, z, sizeof(double)*nCount);
884 : }
885 :
886 : err = pfn_pj_transform( psPJTarget, psPJSource , nCount, 1,
887 5617 : padfTargetX, padfTargetY, (z) ? padfTargetZ : NULL);
888 5617 : if (err == 0)
889 : {
890 750440 : for( i = 0; i < nCount; i++ )
891 : {
892 2075054 : if ( x[i] != HUGE_VAL && y[i] != HUGE_VAL &&
893 744081 : (fabs(padfTargetX[i] - padfOriX[i]) > dfThreshold ||
894 586150 : fabs(padfTargetY[i] - padfOriY[i]) > dfThreshold) )
895 : {
896 158274 : x[i] = HUGE_VAL;
897 158274 : y[i] = HUGE_VAL;
898 : }
899 : }
900 : }
901 : }
902 : }
903 : else
904 : {
905 33481 : err = pfn_pj_transform( psPJSource, psPJTarget, nCount, 1, x, y, z );
906 : }
907 :
908 : /* -------------------------------------------------------------------- */
909 : /* Try to report an error through CPL. Get proj.4 error string */
910 : /* if possible. Try to avoid reporting thousands of error */
911 : /* ... supress further error reporting on this OGRProj4CT if we */
912 : /* have already reported 20 errors. */
913 : /* -------------------------------------------------------------------- */
914 39097 : if( err != 0 )
915 : {
916 0 : if( pabSuccess )
917 0 : memset( pabSuccess, 0, sizeof(int) * nCount );
918 :
919 0 : if( ++nErrorCount < 20 )
920 : {
921 0 : if (pjctx != NULL)
922 : /* pfn_pj_strerrno not yet thread-safe in PROJ 4.8.0 */
923 0 : CPLAcquireMutex(hPROJMutex, 1000.0);
924 :
925 0 : const char *pszError = NULL;
926 0 : if( pfn_pj_strerrno != NULL )
927 0 : pszError = pfn_pj_strerrno( err );
928 :
929 0 : if( pszError == NULL )
930 : CPLError( CE_Failure, CPLE_AppDefined,
931 : "Reprojection failed, err = %d",
932 0 : err );
933 : else
934 0 : CPLError( CE_Failure, CPLE_AppDefined, "%s", pszError );
935 :
936 0 : if (pjctx != NULL)
937 : /* pfn_pj_strerrno not yet thread-safe in PROJ 4.8.0 */
938 0 : CPLReleaseMutex(hPROJMutex);
939 : }
940 0 : else if( nErrorCount == 20 )
941 : {
942 : CPLError( CE_Failure, CPLE_AppDefined,
943 : "Reprojection failed, err = %d, further errors will be supressed on the transform object.",
944 0 : err );
945 : }
946 :
947 0 : if (pjctx == NULL)
948 0 : CPLReleaseMutex(hPROJMutex);
949 0 : return FALSE;
950 : }
951 :
952 39097 : if (pjctx == NULL)
953 0 : CPLReleaseMutex(hPROJMutex);
954 :
955 : /* -------------------------------------------------------------------- */
956 : /* Potentially transform back to degrees. */
957 : /* -------------------------------------------------------------------- */
958 39098 : if( bTargetLatLong )
959 : {
960 3691429 : for( i = 0; i < nCount; i++ )
961 : {
962 3677761 : if( x[i] != HUGE_VAL && y[i] != HUGE_VAL )
963 : {
964 3520005 : x[i] *= dfTargetFromRadians;
965 3520005 : y[i] *= dfTargetFromRadians;
966 : }
967 : }
968 :
969 13668 : if( bTargetWrap )
970 : {
971 3683853 : for( i = 0; i < nCount; i++ )
972 : {
973 3671728 : if( x[i] != HUGE_VAL && y[i] != HUGE_VAL )
974 : {
975 3513972 : if( x[i] < dfTargetWrapLong - 180.0 )
976 808 : x[i] += 360.0;
977 3513164 : else if( x[i] > dfTargetWrapLong + 180 )
978 0 : x[i] -= 360.0;
979 : }
980 : }
981 : }
982 : }
983 :
984 : /* -------------------------------------------------------------------- */
985 : /* Establish error information if pabSuccess provided. */
986 : /* -------------------------------------------------------------------- */
987 39098 : if( pabSuccess )
988 : {
989 4818705 : for( i = 0; i < nCount; i++ )
990 : {
991 4939591 : if( x[i] == HUGE_VAL || y[i] == HUGE_VAL )
992 159983 : pabSuccess[i] = FALSE;
993 : else
994 4619625 : pabSuccess[i] = TRUE;
995 : }
996 : }
997 :
998 39098 : return TRUE;
999 : }
1000 :
1001 : /************************************************************************/
1002 : /* OCTTransformEx() */
1003 : /************************************************************************/
1004 :
1005 0 : int CPL_STDCALL OCTTransformEx( OGRCoordinateTransformationH hTransform,
1006 : int nCount, double *x, double *y, double *z,
1007 : int *pabSuccess )
1008 :
1009 : {
1010 0 : VALIDATE_POINTER1( hTransform, "OCTTransformEx", FALSE );
1011 :
1012 : return ((OGRCoordinateTransformation*) hTransform)->
1013 0 : TransformEx( nCount, x, y, z, pabSuccess );
1014 : }
1015 :
|
