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