1 : /******************************************************************************
2 : * $Id: ogrct.cpp 20079 2010-07-16 21:28:02Z 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 20079 2010-07-16 21:28:02Z 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 68 : static const char* GetProjLibraryName()
160 : {
161 68 : const char *pszLibName = LIBNAME;
162 : #if !defined(WIN32CE)
163 68 : if( CPLGetConfigOption("PROJSO",NULL) != NULL )
164 0 : pszLibName = CPLGetConfigOption("PROJSO",NULL);
165 : #endif
166 68 : return pszLibName;
167 : }
168 :
169 : /************************************************************************/
170 : /* LoadProjLibrary() */
171 : /************************************************************************/
172 :
173 1116 : static int LoadProjLibrary()
174 :
175 : {
176 1116 : CPLMutexHolderD( &hPROJMutex );
177 : static int bTriedToLoad = FALSE;
178 : const char *pszLibName;
179 :
180 1116 : if( bTriedToLoad )
181 1048 : return( pfn_pj_transform != NULL );
182 :
183 68 : bTriedToLoad = TRUE;
184 :
185 68 : 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 68 : CPLPushErrorHandler( CPLQuietErrorHandler );
206 :
207 : pfn_pj_init = (projPJ (*)(int, char**)) CPLGetSymbol( pszLibName,
208 68 : "pj_init" );
209 68 : CPLPopErrorHandler();
210 :
211 68 : if( pfn_pj_init == NULL )
212 0 : return( FALSE );
213 :
214 : pfn_pj_init_plus = (projPJ (*)(const char *))
215 68 : CPLGetSymbol( pszLibName, "pj_init_plus" );
216 : pfn_pj_free = (void (*)(projPJ))
217 68 : CPLGetSymbol( pszLibName, "pj_free" );
218 : pfn_pj_transform = (int (*)(projPJ,projPJ,long,int,double*,
219 : double*,double*))
220 68 : CPLGetSymbol( pszLibName, "pj_transform" );
221 : pfn_pj_get_errno_ref = (int *(*)(void))
222 68 : CPLGetSymbol( pszLibName, "pj_get_errno_ref" );
223 : pfn_pj_strerrno = (char *(*)(int))
224 68 : CPLGetSymbol( pszLibName, "pj_strerrno" );
225 :
226 68 : CPLPushErrorHandler( CPLQuietErrorHandler );
227 : pfn_pj_get_def = (char *(*)(projPJ,int))
228 68 : CPLGetSymbol( pszLibName, "pj_get_def" );
229 : pfn_pj_dalloc = (void (*)(void*))
230 68 : CPLGetSymbol( pszLibName, "pj_dalloc" );
231 :
232 : /* PROJ 4.8.0 symbols */
233 : pfn_pj_ctx_alloc = (projCtx (*)( void ))
234 68 : CPLGetSymbol( pszLibName, "pj_ctx_alloc" );
235 : pfn_pj_ctx_free = (void (*)( projCtx ))
236 68 : CPLGetSymbol( pszLibName, "pj_ctx_free" );
237 : pfn_pj_init_plus_ctx = (projPJ (*)( projCtx, const char * ))
238 68 : CPLGetSymbol( pszLibName, "pj_init_plus_ctx" );
239 : pfn_pj_ctx_get_errno = (int (*)( projCtx ))
240 68 : CPLGetSymbol( pszLibName, "pj_ctx_get_errno" );
241 :
242 68 : CPLPopErrorHandler();
243 68 : CPLErrorReset();
244 : #endif
245 :
246 68 : 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 68 : 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 68 : 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 68 : 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 1062 : char *OCTProj4Normalize( const char *pszProj4Src )
285 :
286 : {
287 : char *pszNewProj4Def, *pszCopy;
288 1062 : projPJ psPJSource = NULL;
289 :
290 1062 : CPLMutexHolderD( &hPROJMutex );
291 :
292 1062 : if( !LoadProjLibrary() || pfn_pj_dalloc == NULL || pfn_pj_get_def == NULL )
293 0 : return CPLStrdup( pszProj4Src );
294 :
295 1062 : psPJSource = pfn_pj_init_plus( pszProj4Src );
296 :
297 1062 : if( psPJSource == NULL )
298 816 : return CPLStrdup( pszProj4Src );
299 :
300 246 : pszNewProj4Def = pfn_pj_get_def( psPJSource, 0 );
301 :
302 246 : pfn_pj_free( psPJSource );
303 :
304 246 : if( pszNewProj4Def == NULL )
305 0 : return CPLStrdup( pszProj4Src );
306 :
307 246 : pszCopy = CPLStrdup( pszNewProj4Def );
308 246 : pfn_pj_dalloc( pszNewProj4Def );
309 :
310 246 : 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 190 : OCTDestroyCoordinateTransformation( OGRCoordinateTransformationH hCT )
327 :
328 : {
329 190 : delete (OGRCoordinateTransformation *) hCT;
330 190 : }
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 236 : void OGRCoordinateTransformation::DestroyCT(OGRCoordinateTransformation* poCT)
353 : {
354 236 : delete poCT;
355 236 : }
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 800 : OGRCreateCoordinateTransformation( OGRSpatialReference *poSource,
381 : OGRSpatialReference *poTarget )
382 :
383 : {
384 : OGRProj4CT *poCT;
385 :
386 800 : 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 800 : poCT = new OGRProj4CT();
396 :
397 800 : if( !poCT->Initialize( poSource, poTarget ) )
398 : {
399 0 : delete poCT;
400 0 : return NULL;
401 : }
402 : else
403 : {
404 800 : 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 184 : OCTNewCoordinateTransformation(
432 : OGRSpatialReferenceH hSourceSRS, OGRSpatialReferenceH hTargetSRS )
433 :
434 : {
435 : return (OGRCoordinateTransformationH)
436 : OGRCreateCoordinateTransformation(
437 : (OGRSpatialReference *) hSourceSRS,
438 184 : (OGRSpatialReference *) hTargetSRS );
439 : }
440 :
441 : /************************************************************************/
442 : /* OGRProj4CT() */
443 : /************************************************************************/
444 :
445 800 : OGRProj4CT::OGRProj4CT()
446 :
447 : {
448 800 : poSRSSource = NULL;
449 800 : poSRSTarget = NULL;
450 800 : psPJSource = NULL;
451 800 : psPJTarget = NULL;
452 :
453 800 : nErrorCount = 0;
454 :
455 800 : bCheckWithInvertProj = FALSE;
456 800 : dfThreshold = 0;
457 :
458 800 : nMaxCount = 0;
459 800 : padfOriX = NULL;
460 800 : padfOriY = NULL;
461 800 : padfOriZ = NULL;
462 800 : padfTargetX = NULL;
463 800 : padfTargetY = NULL;
464 800 : padfTargetZ = NULL;
465 :
466 800 : if (pfn_pj_ctx_alloc != NULL)
467 800 : pjctx = pfn_pj_ctx_alloc();
468 : else
469 0 : pjctx = NULL;
470 800 : }
471 :
472 : /************************************************************************/
473 : /* ~OGRProj4CT() */
474 : /************************************************************************/
475 :
476 800 : OGRProj4CT::~OGRProj4CT()
477 :
478 : {
479 800 : if( poSRSSource != NULL )
480 : {
481 800 : if( poSRSSource->Dereference() <= 0 )
482 800 : delete poSRSSource;
483 : }
484 :
485 800 : if( poSRSTarget != NULL )
486 : {
487 800 : if( poSRSTarget->Dereference() <= 0 )
488 796 : delete poSRSTarget;
489 : }
490 :
491 800 : if (pjctx != NULL)
492 : {
493 800 : pfn_pj_ctx_free(pjctx);
494 :
495 800 : if( psPJSource != NULL )
496 800 : pfn_pj_free( psPJSource );
497 :
498 800 : if( psPJTarget != NULL )
499 800 : 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 800 : CPLFree(padfOriX);
513 800 : CPLFree(padfOriY);
514 800 : CPLFree(padfOriZ);
515 800 : CPLFree(padfTargetX);
516 800 : CPLFree(padfTargetY);
517 800 : CPLFree(padfTargetZ);
518 800 : }
519 :
520 : /************************************************************************/
521 : /* Initialize() */
522 : /************************************************************************/
523 :
524 800 : int OGRProj4CT::Initialize( OGRSpatialReference * poSourceIn,
525 : OGRSpatialReference * poTargetIn )
526 :
527 : {
528 800 : if (pjctx != NULL)
529 : {
530 800 : return InitializeNoLock(poSourceIn, poTargetIn);
531 : }
532 :
533 0 : CPLMutexHolderD( &hPROJMutex );
534 0 : return InitializeNoLock(poSourceIn, poTargetIn);
535 : }
536 :
537 : /************************************************************************/
538 : /* InitializeNoLock() */
539 : /************************************************************************/
540 :
541 800 : int OGRProj4CT::InitializeNoLock( OGRSpatialReference * poSourceIn,
542 : OGRSpatialReference * poTargetIn )
543 :
544 : {
545 800 : if( poSourceIn == NULL || poTargetIn == NULL )
546 0 : return FALSE;
547 :
548 800 : poSRSSource = poSourceIn->Clone();
549 800 : poSRSTarget = poTargetIn->Clone();
550 :
551 800 : bSourceLatLong = poSRSSource->IsGeographic();
552 800 : bTargetLatLong = poSRSTarget->IsGeographic();
553 :
554 : /* -------------------------------------------------------------------- */
555 : /* Setup source and target translations to radians for lat/long */
556 : /* systems. */
557 : /* -------------------------------------------------------------------- */
558 800 : dfSourceToRadians = DEG_TO_RAD;
559 800 : dfSourceFromRadians = RAD_TO_DEG;
560 800 : bSourceWrap = FALSE;
561 800 : dfSourceWrapLong = 0.0;
562 :
563 800 : if( bSourceLatLong )
564 : {
565 466 : OGR_SRSNode *poUNITS = poSRSSource->GetAttrNode( "GEOGCS|UNIT" );
566 466 : if( poUNITS && poUNITS->GetChildCount() >= 2 )
567 : {
568 466 : dfSourceToRadians = atof(poUNITS->GetChild(1)->GetValue());
569 466 : if( dfSourceToRadians == 0.0 )
570 0 : dfSourceToRadians = DEG_TO_RAD;
571 : else
572 466 : dfSourceFromRadians = 1 / dfSourceToRadians;
573 : }
574 : }
575 :
576 800 : dfTargetToRadians = DEG_TO_RAD;
577 800 : dfTargetFromRadians = RAD_TO_DEG;
578 800 : bTargetWrap = FALSE;
579 800 : dfTargetWrapLong = 0.0;
580 :
581 800 : if( bTargetLatLong )
582 : {
583 330 : OGR_SRSNode *poUNITS = poSRSTarget->GetAttrNode( "GEOGCS|UNIT" );
584 330 : if( poUNITS && poUNITS->GetChildCount() >= 2 )
585 : {
586 330 : dfTargetToRadians = atof(poUNITS->GetChild(1)->GetValue());
587 330 : if( dfTargetToRadians == 0.0 )
588 0 : dfTargetToRadians = DEG_TO_RAD;
589 : else
590 330 : dfTargetFromRadians = 1 / dfTargetToRadians;
591 : }
592 : }
593 :
594 : /* -------------------------------------------------------------------- */
595 : /* Preliminary logic to setup wrapping. */
596 : /* -------------------------------------------------------------------- */
597 : const char *pszCENTER_LONG;
598 :
599 800 : 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 800 : pszCENTER_LONG = poSRSSource->GetExtension( "GEOGCS", "CENTER_LONG" );
608 800 : if( pszCENTER_LONG != NULL )
609 : {
610 72 : dfSourceWrapLong = atof(pszCENTER_LONG);
611 72 : bSourceWrap = TRUE;
612 72 : CPLDebug( "OGRCT", "Wrap source at %g.", dfSourceWrapLong );
613 : }
614 :
615 800 : pszCENTER_LONG = poSRSTarget->GetExtension( "GEOGCS", "CENTER_LONG" );
616 800 : if( pszCENTER_LONG != NULL )
617 : {
618 72 : dfTargetWrapLong = atof(pszCENTER_LONG);
619 72 : bTargetWrap = TRUE;
620 72 : CPLDebug( "OGRCT", "Wrap target at %g.", dfTargetWrapLong );
621 : }
622 :
623 800 : bCheckWithInvertProj = CSLTestBoolean(CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", "NO" ));
624 :
625 : /* The threshold is rather experimental... Works well with the cases of ticket #2305 */
626 800 : if (bSourceLatLong)
627 466 : 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 334 : 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 800 : char *pszProj4Defn = NULL;
642 :
643 800 : if( poSRSSource->exportToProj4( &pszProj4Defn ) != OGRERR_NONE )
644 : {
645 0 : CPLFree( pszProj4Defn );
646 0 : return FALSE;
647 : }
648 :
649 800 : 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 800 : if (pjctx)
659 800 : psPJSource = pfn_pj_init_plus_ctx( pjctx, pszProj4Defn );
660 : else
661 0 : psPJSource = pfn_pj_init_plus( pszProj4Defn );
662 :
663 800 : 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 800 : if( nDebugReportCount < 10 )
693 178 : CPLDebug( "OGRCT", "Source: %s", pszProj4Defn );
694 :
695 800 : CPLFree( pszProj4Defn );
696 :
697 800 : if( psPJSource == NULL )
698 0 : return FALSE;
699 :
700 : /* -------------------------------------------------------------------- */
701 : /* Establish PROJ.4 handle for target if projection. */
702 : /* -------------------------------------------------------------------- */
703 800 : pszProj4Defn = NULL;
704 :
705 800 : if( poSRSTarget->exportToProj4( &pszProj4Defn ) != OGRERR_NONE )
706 : {
707 0 : CPLFree( pszProj4Defn );
708 0 : return FALSE;
709 : }
710 :
711 800 : 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 800 : if (pjctx)
721 800 : psPJTarget = pfn_pj_init_plus_ctx( pjctx, pszProj4Defn );
722 : else
723 0 : psPJTarget = pfn_pj_init_plus( pszProj4Defn );
724 :
725 800 : if( psPJTarget == NULL )
726 : CPLError( CE_Failure, CPLE_NotSupported,
727 : "Failed to initialize PROJ.4 with `%s'.",
728 0 : pszProj4Defn );
729 :
730 800 : if( nDebugReportCount < 10 )
731 : {
732 178 : CPLDebug( "OGRCT", "Target: %s", pszProj4Defn );
733 178 : nDebugReportCount++;
734 : }
735 :
736 800 : CPLFree( pszProj4Defn );
737 :
738 800 : if( psPJTarget == NULL )
739 0 : return FALSE;
740 :
741 800 : return TRUE;
742 : }
743 :
744 : /************************************************************************/
745 : /* GetSourceCS() */
746 : /************************************************************************/
747 :
748 56 : OGRSpatialReference *OGRProj4CT::GetSourceCS()
749 :
750 : {
751 56 : return poSRSSource;
752 : }
753 :
754 : /************************************************************************/
755 : /* GetTargetCS() */
756 : /************************************************************************/
757 :
758 148 : OGRSpatialReference *OGRProj4CT::GetTargetCS()
759 :
760 : {
761 148 : return poSRSTarget;
762 : }
763 :
764 : /************************************************************************/
765 : /* Transform() */
766 : /* */
767 : /* This is a small wrapper for the extended transform version. */
768 : /************************************************************************/
769 :
770 1740 : int OGRProj4CT::Transform( int nCount, double *x, double *y, double *z )
771 :
772 : {
773 1740 : int *pabSuccess = (int *) CPLMalloc(sizeof(int) * nCount );
774 : int bOverallSuccess, i;
775 :
776 1740 : bOverallSuccess = TransformEx( nCount, x, y, z, pabSuccess );
777 :
778 3484 : for( i = 0; i < nCount; i++ )
779 : {
780 1744 : if( !pabSuccess[i] )
781 : {
782 0 : bOverallSuccess = FALSE;
783 0 : break;
784 : }
785 : }
786 :
787 1740 : CPLFree( pabSuccess );
788 :
789 1740 : return bOverallSuccess;
790 : }
791 :
792 : /************************************************************************/
793 : /* OCTTransform() */
794 : /************************************************************************/
795 :
796 644 : int CPL_STDCALL OCTTransform( OGRCoordinateTransformationH hTransform,
797 : int nCount, double *x, double *y, double *z )
798 :
799 : {
800 : return ((OGRCoordinateTransformation*) hTransform)->
801 644 : Transform( nCount, x, y,z );
802 : }
803 :
804 : /************************************************************************/
805 : /* TransformEx() */
806 : /************************************************************************/
807 :
808 75072 : int OGRProj4CT::TransformEx( int nCount, double *x, double *y, double *z,
809 : int *pabSuccess )
810 :
811 : {
812 : int err, i;
813 :
814 : /* -------------------------------------------------------------------- */
815 : /* Potentially transform to radians. */
816 : /* -------------------------------------------------------------------- */
817 75072 : if( bSourceLatLong )
818 : {
819 20652 : if( bSourceWrap )
820 : {
821 29652 : for( i = 0; i < nCount; i++ )
822 : {
823 21022 : if( x[i] != HUGE_VAL && y[i] != HUGE_VAL )
824 : {
825 19708 : if( x[i] < dfSourceWrapLong - 180.0 )
826 0 : x[i] += 360.0;
827 19708 : else if( x[i] > dfSourceWrapLong + 180 )
828 0 : x[i] -= 360.0;
829 : }
830 : }
831 : }
832 :
833 7587402 : for( i = 0; i < nCount; i++ )
834 : {
835 7566750 : if( x[i] != HUGE_VAL )
836 : {
837 7565436 : x[i] *= dfSourceToRadians;
838 7565436 : y[i] *= dfSourceToRadians;
839 : }
840 : }
841 : }
842 :
843 : /* -------------------------------------------------------------------- */
844 : /* Do the transformation using PROJ.4. */
845 : /* -------------------------------------------------------------------- */
846 75072 : if (pjctx == NULL)
847 : {
848 : /* The mutex has already been created */
849 0 : CPLAssert(hPROJMutex != NULL);
850 0 : CPLAcquireMutex(hPROJMutex, 1000.0);
851 : }
852 :
853 75072 : if (bCheckWithInvertProj)
854 : {
855 : /* For some projections, we cannot detect if we are trying to reproject */
856 : /* coordinates outside the validity area of the projection. So let's do */
857 : /* the reverse reprojection and compare with the source coordinates */
858 11234 : if (nCount > nMaxCount)
859 : {
860 42 : nMaxCount = nCount;
861 42 : padfOriX = (double*) CPLRealloc(padfOriX, sizeof(double)*nCount);
862 42 : padfOriY = (double*) CPLRealloc(padfOriY, sizeof(double)*nCount);
863 42 : padfOriZ = (double*) CPLRealloc(padfOriZ, sizeof(double)*nCount);
864 42 : padfTargetX = (double*) CPLRealloc(padfTargetX, sizeof(double)*nCount);
865 42 : padfTargetY = (double*) CPLRealloc(padfTargetY, sizeof(double)*nCount);
866 42 : padfTargetZ = (double*) CPLRealloc(padfTargetZ, sizeof(double)*nCount);
867 : }
868 11234 : memcpy(padfOriX, x, sizeof(double)*nCount);
869 11234 : memcpy(padfOriY, y, sizeof(double)*nCount);
870 11234 : if (z)
871 : {
872 11234 : memcpy(padfOriZ, z, sizeof(double)*nCount);
873 : }
874 11234 : err = pfn_pj_transform( psPJSource, psPJTarget, nCount, 1, x, y, z );
875 11234 : if (err == 0)
876 : {
877 11234 : memcpy(padfTargetX, x, sizeof(double)*nCount);
878 11234 : memcpy(padfTargetY, y, sizeof(double)*nCount);
879 11234 : if (z)
880 : {
881 11234 : memcpy(padfTargetZ, z, sizeof(double)*nCount);
882 : }
883 :
884 : err = pfn_pj_transform( psPJTarget, psPJSource , nCount, 1,
885 11234 : padfTargetX, padfTargetY, (z) ? padfTargetZ : NULL);
886 11234 : if (err == 0)
887 : {
888 1500880 : for( i = 0; i < nCount; i++ )
889 : {
890 4150108 : if ( x[i] != HUGE_VAL && y[i] != HUGE_VAL &&
891 1488162 : (fabs(padfTargetX[i] - padfOriX[i]) > dfThreshold ||
892 1172300 : fabs(padfTargetY[i] - padfOriY[i]) > dfThreshold) )
893 : {
894 316548 : x[i] = HUGE_VAL;
895 316548 : y[i] = HUGE_VAL;
896 : }
897 : }
898 : }
899 : }
900 : }
901 : else
902 : {
903 63838 : err = pfn_pj_transform( psPJSource, psPJTarget, nCount, 1, x, y, z );
904 : }
905 :
906 : /* -------------------------------------------------------------------- */
907 : /* Try to report an error through CPL. Get proj.4 error string */
908 : /* if possible. Try to avoid reporting thousands of error */
909 : /* ... supress further error reporting on this OGRProj4CT if we */
910 : /* have already reported 20 errors. */
911 : /* -------------------------------------------------------------------- */
912 75072 : if( err != 0 )
913 : {
914 0 : if( pabSuccess )
915 0 : memset( pabSuccess, 0, sizeof(int) * nCount );
916 :
917 0 : if( ++nErrorCount < 20 )
918 : {
919 0 : if (pjctx != NULL)
920 : /* pfn_pj_strerrno not yet thread-safe in PROJ 4.8.0 */
921 0 : CPLAcquireMutex(hPROJMutex, 1000.0);
922 :
923 0 : const char *pszError = NULL;
924 0 : if( pfn_pj_strerrno != NULL )
925 0 : pszError = pfn_pj_strerrno( err );
926 :
927 0 : if( pszError == NULL )
928 : CPLError( CE_Failure, CPLE_AppDefined,
929 : "Reprojection failed, err = %d",
930 0 : err );
931 : else
932 0 : CPLError( CE_Failure, CPLE_AppDefined, "%s", pszError );
933 :
934 0 : if (pjctx != NULL)
935 : /* pfn_pj_strerrno not yet thread-safe in PROJ 4.8.0 */
936 0 : CPLReleaseMutex(hPROJMutex);
937 : }
938 0 : else if( nErrorCount == 20 )
939 : {
940 : CPLError( CE_Failure, CPLE_AppDefined,
941 : "Reprojection failed, err = %d, further errors will be supressed on the transform object.",
942 0 : err );
943 : }
944 :
945 0 : if (pjctx == NULL)
946 0 : CPLReleaseMutex(hPROJMutex);
947 0 : return FALSE;
948 : }
949 :
950 75072 : if (pjctx == NULL)
951 0 : CPLReleaseMutex(hPROJMutex);
952 :
953 : /* -------------------------------------------------------------------- */
954 : /* Potentially transform back to degrees. */
955 : /* -------------------------------------------------------------------- */
956 75072 : if( bTargetLatLong )
957 : {
958 7374776 : for( i = 0; i < nCount; i++ )
959 : {
960 7347834 : if( x[i] != HUGE_VAL && y[i] != HUGE_VAL )
961 : {
962 7032322 : x[i] *= dfTargetFromRadians;
963 7032322 : y[i] *= dfTargetFromRadians;
964 : }
965 : }
966 :
967 26942 : if( bTargetWrap )
968 : {
969 7367706 : for( i = 0; i < nCount; i++ )
970 : {
971 7343456 : if( x[i] != HUGE_VAL && y[i] != HUGE_VAL )
972 : {
973 7027944 : if( x[i] < dfTargetWrapLong - 180.0 )
974 1616 : x[i] += 360.0;
975 7026328 : else if( x[i] > dfTargetWrapLong + 180 )
976 0 : x[i] -= 360.0;
977 : }
978 : }
979 : }
980 : }
981 :
982 : /* -------------------------------------------------------------------- */
983 : /* Establish error information if pabSuccess provided. */
984 : /* -------------------------------------------------------------------- */
985 75072 : if( pabSuccess )
986 : {
987 9616386 : for( i = 0; i < nCount; i++ )
988 : {
989 9861280 : if( x[i] == HUGE_VAL || y[i] == HUGE_VAL )
990 319966 : pabSuccess[i] = FALSE;
991 : else
992 9221348 : pabSuccess[i] = TRUE;
993 : }
994 : }
995 :
996 75072 : return TRUE;
997 : }
998 :
999 : /************************************************************************/
1000 : /* OCTTransformEx() */
1001 : /************************************************************************/
1002 :
1003 0 : int CPL_STDCALL OCTTransformEx( OGRCoordinateTransformationH hTransform,
1004 : int nCount, double *x, double *y, double *z,
1005 : int *pabSuccess )
1006 :
1007 : {
1008 : return ((OGRCoordinateTransformation*) hTransform)->
1009 0 : TransformEx( nCount, x, y, z, pabSuccess );
1010 : }
1011 :
|
