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