1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // Projection_Proj4 //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // crs_transform.cpp //
14 // //
15 // Copyright (C) 2010 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "crs_transform.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
59 #if PROJ_VERSION_MAJOR < 6
60 extern "C" {
61 #include <projects.h>
62 }
63
64 #define PROJ4_FREE(p) if( p ) { pj_free((PJ *)p); p = NULL; }
65
66 //---------------------------------------------------------
67 #else
68 #include <proj.h>
69
70 #if PROJ_VERSION_MINOR < 2
71 #define PROJ4_FREE(p) if( p ) { proj_destroy((PJ *)p); p = NULL; }
72 #else
73 #define PROJ4_FREE(p) if( p ) { proj_destroy((PJ *)p); p = NULL; proj_cleanup(); }
74 #endif
75 #endif
76
77
78 ///////////////////////////////////////////////////////////
79 // //
80 ///////////////////////////////////////////////////////////
81
82 //---------------------------------------------------------
CSG_CRSProjector(void)83 CSG_CRSProjector::CSG_CRSProjector(void)
84 {
85 _On_Construction();
86 }
87
88 //---------------------------------------------------------
CSG_CRSProjector(const CSG_CRSProjector & Projector)89 CSG_CRSProjector::CSG_CRSProjector(const CSG_CRSProjector &Projector)
90 {
91 _On_Construction();
92
93 Create(Projector);
94 }
95
96 //---------------------------------------------------------
~CSG_CRSProjector(void)97 CSG_CRSProjector::~CSG_CRSProjector(void)
98 {
99 Destroy();
100
101 #if PROJ_VERSION_MAJOR >= 6
102 proj_context_destroy((PJ_CONTEXT *)m_pContext);
103
104 #if PROJ_VERSION_MINOR >= 2
105 proj_cleanup();
106 #endif
107 #endif
108 }
109
110 //---------------------------------------------------------
_On_Construction(void)111 void CSG_CRSProjector::_On_Construction(void)
112 {
113 m_pSource = NULL;
114 m_pTarget = NULL;
115 m_pGCS = NULL;
116
117 m_bInverse = false;
118
119 m_Copies = NULL;
120 m_nCopies = 0;
121
122 #if PROJ_VERSION_MAJOR < 6
123 m_pContext = NULL;
124 #else
125 m_pContext = proj_context_create();
126 #endif
127 }
128
129 //---------------------------------------------------------
Create(const CSG_CRSProjector & Projector)130 bool CSG_CRSProjector::Create(const CSG_CRSProjector &Projector)
131 {
132 Destroy();
133
134 Set_Source(Projector.m_Source);
135 Set_Target(Projector.m_Target);
136
137 Set_Inverse(Projector.m_bInverse);
138
139 Set_Precise_Mode(Projector.Get_Precise_Mode());
140
141 return( true );
142 }
143
144 //---------------------------------------------------------
Destroy(void)145 bool CSG_CRSProjector::Destroy(void)
146 {
147 m_bInverse = false;
148
149 PROJ4_FREE(m_pSource);
150 PROJ4_FREE(m_pTarget);
151 PROJ4_FREE(m_pGCS );
152
153 Set_Copies();
154
155 return( true );
156 }
157
158
159 ///////////////////////////////////////////////////////////
160 // //
161 ///////////////////////////////////////////////////////////
162
163 //---------------------------------------------------------
Set_Copies(int nCopies)164 bool CSG_CRSProjector::Set_Copies(int nCopies)
165 {
166 if( m_Copies )
167 {
168 delete[](m_Copies);
169
170 m_Copies = NULL;
171 m_nCopies = 0;
172 }
173
174 if( nCopies > 1 )
175 {
176 m_Copies = new CSG_CRSProjector[m_nCopies = nCopies - 1];
177
178 for(int i=0; i<m_nCopies; i++)
179 {
180 m_Copies[i].Create(*this);
181 }
182 }
183
184 return( true );
185 }
186
187 //---------------------------------------------------------
operator [](int iCopy)188 CSG_CRSProjector & CSG_CRSProjector::operator [] (int iCopy)
189 {
190 if( iCopy > 0 && iCopy <= m_nCopies )
191 {
192 return( m_Copies[iCopy - 1] );
193 }
194
195 return( *this );
196 }
197
198
199 ///////////////////////////////////////////////////////////
200 // //
201 ///////////////////////////////////////////////////////////
202
203 //---------------------------------------------------------
Get_Version(void)204 CSG_String CSG_CRSProjector::Get_Version(void)
205 {
206 #if PROJ_VERSION_MAJOR < 6
207 return( pj_release );
208 #else
209 return( CSG_String::Format("%d.%d.%d", PROJ_VERSION_MAJOR, PROJ_VERSION_MINOR, PROJ_VERSION_PATCH) );
210 #endif
211 }
212
213 //---------------------------------------------------------
Get_Description(void)214 CSG_String CSG_CRSProjector::Get_Description(void)
215 {
216 CSG_String s;
217
218 s += _TL("Projection routines make use of the Proj.4 Cartographic Projections library.");
219 s += "\n";
220 s += _TW("Proj.4 was originally developed by Gerald Evenden and later continued by the "
221 "United States Department of the Interior, Geological Survey (USGS).");
222 s += "\n";
223 s += _TL("Proj.4 Version is ") + Get_Version();
224 s += "\n";
225 s += "<a target=\"_blank\" href=\"http://trac.osgeo.org/proj/\">Proj.4 Homepage</a>";
226
227 return( s );
228 }
229
230
231 ///////////////////////////////////////////////////////////
232 // //
233 ///////////////////////////////////////////////////////////
234
235 //---------------------------------------------------------
_Set_Projection(const CSG_Projection & Projection,void ** ppProjection,bool bInverse)236 bool CSG_CRSProjector::_Set_Projection(const CSG_Projection &Projection, void **ppProjection, bool bInverse)
237 {
238 PROJ4_FREE(*ppProjection);
239
240 //-------------------------------------------------
241 #if PROJ_VERSION_MAJOR < 6
242 if( (*ppProjection = pj_init_plus(Projection.Get_Proj4())) == NULL )
243 {
244 CSG_String Error(pj_strerrno(pj_errno));
245 #else
246 if( (*ppProjection = proj_create((PJ_CONTEXT *)m_pContext, Projection.Get_Proj4())) == NULL )
247 {
248 CSG_String Error(proj_errno_string(proj_errno((PJ *)(*ppProjection))));
249 #endif
250
251 SG_UI_Msg_Add_Error(CSG_String::Format("Proj4 [%s]: %s", _TL("initialization"), Error.c_str()));
252
253 return( false );
254 }
255
256 //-------------------------------------------------
257 #if PROJ_VERSION_MAJOR < 6
258 if( bInverse && ((PJ *)(*ppProjection))->inv == NULL )
259 #else
260 if( bInverse && !proj_pj_info((PJ *)(*ppProjection)).has_inverse )
261 #endif
262 {
263 SG_UI_Msg_Add_Error(CSG_String::Format("Proj4 [%s]: %s", _TL("initialization"), _TL("inverse transformation not available")));
264
265 return( false );
266 }
267
268 return( true );
269 }
270
271 //---------------------------------------------------------
272 bool CSG_CRSProjector::Set_Source(const CSG_Projection &Projection)
273 {
274 return( Projection.is_Okay() && _Set_Projection(Projection, &m_pSource, true) && m_Source.Create(Projection) );
275 }
276
277 //---------------------------------------------------------
278 bool CSG_CRSProjector::Set_Target(const CSG_Projection &Projection)
279 {
280 return( Projection.is_Okay() && _Set_Projection(Projection, &m_pTarget, false) && m_Target.Create(Projection) );
281 }
282
283
284 ///////////////////////////////////////////////////////////
285 // //
286 ///////////////////////////////////////////////////////////
287
288 //---------------------------------------------------------
289 bool CSG_CRSProjector::Set_Inverse(bool bOn)
290 {
291 if( m_bInverse == bOn )
292 {
293 return( true );
294 }
295
296 #if PROJ_VERSION_MAJOR < 6
297 if( m_pTarget && ((PJ *)m_pTarget)->inv )
298 #else
299 if( m_pTarget && proj_pj_info((PJ *)m_pTarget).has_inverse )
300 #endif
301 {
302 m_bInverse = bOn;
303
304 void *pTMP = m_pSource;
305 m_pSource = m_pTarget;
306 m_pTarget = pTMP;
307
308 return( true );
309 }
310
311 SG_UI_Msg_Add_Error(CSG_String::Format("Proj4 [%s]: %s", _TL("initialization"), _TL("inverse transformation not available")));
312
313 return( false );
314 }
315
316 //---------------------------------------------------------
317 bool CSG_CRSProjector::Set_Precise_Mode(bool bOn)
318 {
319 if( bOn )
320 {
321 if( m_pGCS == NULL )
322 {
323 #if PROJ_VERSION_MAJOR < 6
324 return( (m_pGCS = pj_init_plus("+proj=longlat +datum=WGS84")) != NULL );
325 #else
326 return( (m_pGCS = proj_create((PJ_CONTEXT *)m_pContext, "+proj=longlat +datum=WGS84")) != NULL );
327 #endif
328 }
329 }
330 else
331 {
332 PROJ4_FREE(m_pGCS);
333 }
334
335 return( true );
336 }
337
338
339 ///////////////////////////////////////////////////////////
340 // //
341 ///////////////////////////////////////////////////////////
342
343 //---------------------------------------------------------
344 bool CSG_CRSProjector::Get_Projection(double &x, double &y) const
345 {
346 if( !m_pSource || !m_pTarget )
347 {
348 return( false );
349 }
350
351 #if PROJ_VERSION_MAJOR < 6
352 if( pj_is_latlong((PJ *)m_pSource) )
353 #else
354 if( proj_angular_output((PJ *)m_pSource, PJ_FWD) )
355 #endif
356 {
357 x *= M_DEG_TO_RAD;
358 y *= M_DEG_TO_RAD;
359 }
360
361 #if PROJ_VERSION_MAJOR < 6
362 if( m_pGCS ) // precise datum conversion
363 {
364 if( pj_transform((PJ *)m_pSource, (PJ *)m_pGCS , 1, 0, &x, &y, NULL) != 0
365 || pj_transform((PJ *)m_pGCS , (PJ *)m_pTarget, 1, 0, &x, &y, NULL) != 0 )
366 {
367 return( false );
368 }
369 }
370 else // direct projection
371 {
372 if( pj_transform((PJ *)m_pSource, (PJ *)m_pTarget, 1, 0, &x, &y, NULL) != 0 )
373 {
374 return( false );
375 }
376 }
377 #else
378 PJ_COORD c = proj_coord(x, y, 0, 0);
379
380 c = proj_trans((PJ *)m_pSource, PJ_INV, c); if( proj_errno((PJ *)m_pSource) ) { proj_errno_reset((PJ *)m_pSource); return( false ); }
381 c = proj_trans((PJ *)m_pTarget, PJ_FWD, c); if( proj_errno((PJ *)m_pTarget) ) { proj_errno_reset((PJ *)m_pTarget); return( false ); }
382
383 x = c.v[0];
384 y = c.v[1];
385 #endif
386
387 #if PROJ_VERSION_MAJOR < 6
388 if( pj_is_latlong((PJ *)m_pTarget) )
389 #else
390 if( proj_angular_output((PJ *)m_pTarget, PJ_FWD) )
391 #endif
392 {
393 x *= M_RAD_TO_DEG;
394 y *= M_RAD_TO_DEG;
395 }
396
397 return( true );
398 }
399
400 //---------------------------------------------------------
401 bool CSG_CRSProjector::Get_Projection(TSG_Point &Point) const
402 {
403 return( Get_Projection(Point.x, Point.y) );
404 }
405
406 //---------------------------------------------------------
407 bool CSG_CRSProjector::Get_Projection(CSG_Point &Point) const
408 {
409 double x = Point.Get_X();
410 double y = Point.Get_Y();
411
412 if( Get_Projection(x, y) )
413 {
414 Point.Assign(x, y);
415
416 return( true );
417 }
418
419 return( false );
420 }
421
422 //---------------------------------------------------------
423 bool CSG_CRSProjector::Get_Projection(double &x, double &y, double &z) const
424 {
425 if( !m_pSource || !m_pTarget )
426 {
427 return( false );
428 }
429
430 #if PROJ_VERSION_MAJOR < 6
431 if( pj_is_latlong((PJ *)m_pSource) )
432 #else
433 if( proj_angular_output((PJ *)m_pSource, PJ_FWD) )
434 #endif
435 {
436 x *= M_DEG_TO_RAD;
437 y *= M_DEG_TO_RAD;
438 }
439
440 #if PROJ_VERSION_MAJOR < 6
441 if( m_pGCS ) // precise datum conversion
442 {
443 if( pj_transform((PJ *)m_pSource, (PJ *)m_pGCS , 1, 0, &x, &y, &z) != 0
444 || pj_transform((PJ *)m_pGCS , (PJ *)m_pTarget, 1, 0, &x, &y, &z) != 0 )
445 {
446 return( false );
447 }
448 }
449 else // direct projection
450 {
451 if( pj_transform((PJ *)m_pSource, (PJ *)m_pTarget, 1, 0, &x, &y, &z) != 0 )
452 {
453 return( false );
454 }
455 }
456 #else
457 PJ_COORD c = proj_coord(x, y, z, 0);
458
459 c = proj_trans((PJ *)m_pSource, PJ_INV, c); if( proj_errno((PJ *)m_pSource) ) { proj_errno_reset((PJ *)m_pSource); return( false ); }
460 c = proj_trans((PJ *)m_pTarget, PJ_FWD, c); if( proj_errno((PJ *)m_pTarget) ) { proj_errno_reset((PJ *)m_pTarget); return( false ); }
461
462 x = c.v[0];
463 y = c.v[1];
464 z = c.v[2];
465 #endif
466
467 #if PROJ_VERSION_MAJOR < 6
468 if( pj_is_latlong((PJ *)m_pTarget) )
469 #else
470 if( proj_angular_output((PJ *)m_pTarget, PJ_FWD) )
471 #endif
472 {
473 x *= M_RAD_TO_DEG;
474 y *= M_RAD_TO_DEG;
475 }
476
477 return( true );
478 }
479
480 //---------------------------------------------------------
481 bool CSG_CRSProjector::Get_Projection(TSG_Point_Z &Point) const
482 {
483 return( Get_Projection(Point.x, Point.y, Point.z) );
484 }
485
486 //---------------------------------------------------------
487 bool CSG_CRSProjector::Get_Projection(CSG_Point_Z &Point) const
488 {
489 double x = Point.Get_X();
490 double y = Point.Get_Y();
491 double z = Point.Get_Z();
492
493 if( Get_Projection(x, y, z) )
494 {
495 Point.Assign(x, y, z);
496
497 return( true );
498 }
499
500 return( false );
501 }
502
503 ///////////////////////////////////////////////////////////
504 // //
505 // //
506 // //
507 ///////////////////////////////////////////////////////////
508
509 //---------------------------------------------------------
510