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