Line data Source code
1 : /******************************************************************************
2 : * $Id$
3 : *
4 : * Project: libgeotiff
5 : * Purpose: Code to normalize PCS and other composite codes in a GeoTIFF file.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, Frank Warmerdam
10 : * Copyright (c) 2018, Even Rouault <even.rouault at spatialys.com>
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : *****************************************************************************/
30 :
31 : #include <assert.h>
32 : #include <stdio.h>
33 : #include <stdlib.h>
34 : #include <string.h>
35 :
36 : #include "cpl_serv.h"
37 : #include "geo_tiffp.h"
38 : #include "geovalues.h"
39 : #include "geo_normalize.h"
40 : #include "geo_keyp.h"
41 :
42 : #include "proj.h"
43 :
44 : #ifndef KvUserDefined
45 : # define KvUserDefined 32767
46 : #endif
47 :
48 : #ifndef M_PI
49 : # define M_PI 3.14159265358979323846
50 : #endif
51 :
52 : /* EPSG Codes for projection parameters. Unfortunately, these bear no
53 : relationship to the GeoTIFF codes even though the names are so similar. */
54 :
55 : #define EPSGNatOriginLat 8801
56 : #define EPSGNatOriginLong 8802
57 : #define EPSGNatOriginScaleFactor 8805
58 : #define EPSGFalseEasting 8806
59 : #define EPSGFalseNorthing 8807
60 : #define EPSGProjCenterLat 8811
61 : #define EPSGProjCenterLong 8812
62 : #define EPSGAzimuth 8813
63 : #define EPSGAngleRectifiedToSkewedGrid 8814
64 : #define EPSGInitialLineScaleFactor 8815
65 : #define EPSGProjCenterEasting 8816
66 : #define EPSGProjCenterNorthing 8817
67 : #define EPSGPseudoStdParallelLat 8818
68 : #define EPSGPseudoStdParallelScaleFactor 8819
69 : #define EPSGFalseOriginLat 8821
70 : #define EPSGFalseOriginLong 8822
71 : #define EPSGStdParallel1Lat 8823
72 : #define EPSGStdParallel2Lat 8824
73 : #define EPSGFalseOriginEasting 8826
74 : #define EPSGFalseOriginNorthing 8827
75 : #define EPSGSphericalOriginLat 8828
76 : #define EPSGSphericalOriginLong 8829
77 : #define EPSGInitialLongitude 8830
78 : #define EPSGZoneWidth 8831
79 : #define EPSGLatOfStdParallel 8832
80 : #define EPSGOriginLong 8833
81 : #define EPSGTopocentricOriginLat 8834
82 : #define EPSGTopocentricOriginLong 8835
83 : #define EPSGTopocentricOriginHeight 8836
84 :
85 : #define CT_Ext_Mercator_2SP -CT_Mercator
86 :
87 : #ifndef CPL_INLINE
88 : # if (defined(__GNUC__) && !defined(__NO_INLINE__)) || defined(_MSC_VER)
89 : # define HAS_CPL_INLINE 1
90 : # define CPL_INLINE __inline
91 : # elif defined(__SUNPRO_CC)
92 : # define HAS_CPL_INLINE 1
93 : # define CPL_INLINE inline
94 : # else
95 : # define CPL_INLINE
96 : # endif
97 : #endif
98 :
99 : #ifndef CPL_UNUSED
100 : #if defined(__GNUC__) && __GNUC__ >= 4
101 : # define CPL_UNUSED __attribute((__unused__))
102 : #else
103 : /* TODO: add cases for other compilers */
104 : # define CPL_UNUSED
105 : #endif
106 : #endif
107 :
108 13556 : CPL_INLINE static void CPL_IGNORE_RET_VAL_INT(CPL_UNUSED int unused) {}
109 :
110 : /************************************************************************/
111 : /* GTIFKeyGetSSHORT() */
112 : /************************************************************************/
113 :
114 : // Geotiff SHORT keys are supposed to be unsigned, but geo_normalize interface
115 : // uses signed short...
116 57597 : static int GTIFKeyGetSSHORT( GTIF *gtif, geokey_t key, short* pnVal )
117 : {
118 : unsigned short sVal;
119 57597 : if( GTIFKeyGetSHORT(gtif, key, &sVal, 0, 1) == 1 )
120 : {
121 17696 : memcpy(pnVal, &sVal, 2);
122 17696 : return 1;
123 : }
124 39901 : return 0;
125 : }
126 :
127 : /************************************************************************/
128 : /* GTIFGetPCSInfo() */
129 : /************************************************************************/
130 :
131 5894 : int GTIFGetPCSInfoEx( void* ctxIn,
132 : int nPCSCode, char **ppszEPSGName,
133 : short *pnProjOp, short *pnUOMLengthCode,
134 : short *pnGeogCS )
135 :
136 : {
137 : int nDatum;
138 : int nZone;
139 :
140 : /* Deal with a few well known CRS */
141 5894 : int Proj = GTIFPCSToMapSys( nPCSCode, &nDatum, &nZone );
142 5894 : if ((Proj == MapSys_UTM_North || Proj == MapSys_UTM_South) &&
143 5570 : nDatum != KvUserDefined)
144 : {
145 5570 : const char* pszDatumName = NULL;
146 5570 : switch (nDatum)
147 : {
148 5264 : case GCS_NAD27: pszDatumName = "NAD27"; break;
149 10 : case GCS_NAD83: pszDatumName = "NAD83"; break;
150 0 : case GCS_WGS_72: pszDatumName = "WGS 72"; break;
151 0 : case GCS_WGS_72BE: pszDatumName = "WGS 72BE"; break;
152 296 : case GCS_WGS_84: pszDatumName = "WGS 84"; break;
153 0 : default: break;
154 : }
155 :
156 5570 : if (pszDatumName)
157 : {
158 5570 : if (ppszEPSGName)
159 : {
160 : char szEPSGName[64];
161 2785 : snprintf(
162 : szEPSGName, sizeof(szEPSGName), "%s / UTM zone %d%c",
163 : pszDatumName, nZone,
164 : (Proj == MapSys_UTM_North) ? 'N' : 'S');
165 2785 : *ppszEPSGName = CPLStrdup(szEPSGName);
166 : }
167 :
168 5570 : if (pnProjOp)
169 2785 : *pnProjOp = (short) (((Proj == MapSys_UTM_North) ? Proj_UTM_zone_1N - 1 : Proj_UTM_zone_1S - 1) + nZone);
170 :
171 5570 : if (pnUOMLengthCode)
172 2785 : *pnUOMLengthCode = 9001; /* Linear_Meter */
173 :
174 5570 : if (pnGeogCS)
175 2785 : *pnGeogCS = (short) nDatum;
176 :
177 5570 : return TRUE;
178 : }
179 : }
180 :
181 324 : if( nPCSCode == KvUserDefined )
182 0 : return FALSE;
183 :
184 324 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
185 : {
186 : char szCode[12];
187 :
188 324 : snprintf(szCode, sizeof(szCode), "%d", nPCSCode);
189 324 : PJ* proj_crs = proj_create_from_database(
190 : ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
191 324 : if( !proj_crs )
192 : {
193 0 : return FALSE;
194 : }
195 :
196 324 : if( proj_get_type(proj_crs) != PJ_TYPE_PROJECTED_CRS )
197 : {
198 0 : proj_destroy(proj_crs);
199 0 : return FALSE;
200 : }
201 :
202 324 : if( ppszEPSGName )
203 : {
204 162 : const char* pszName = proj_get_name(proj_crs);
205 162 : if( !pszName )
206 : {
207 : // shouldn't happen
208 0 : proj_destroy(proj_crs);
209 0 : return FALSE;
210 : }
211 162 : *ppszEPSGName = CPLStrdup(pszName);
212 : }
213 :
214 324 : if( pnProjOp )
215 : {
216 162 : PJ* conversion = proj_crs_get_coordoperation(
217 : ctx, proj_crs);
218 162 : if( !conversion )
219 : {
220 : // shouldn't happen except out of memory
221 0 : proj_destroy(proj_crs);
222 0 : return FALSE;
223 : }
224 :
225 : {
226 162 : const char* pszConvCode = proj_get_id_code(conversion, 0);
227 162 : assert( pszConvCode );
228 162 : *pnProjOp = (short) atoi(pszConvCode);
229 : }
230 :
231 162 : proj_destroy(conversion);
232 : }
233 :
234 324 : if( pnUOMLengthCode )
235 : {
236 162 : PJ* coordSys = proj_crs_get_coordinate_system(
237 : ctx, proj_crs);
238 162 : if( !coordSys )
239 : {
240 : // shouldn't happen except out of memory
241 0 : proj_destroy(proj_crs);
242 0 : return FALSE;
243 : }
244 :
245 : {
246 162 : const char* pszUnitCode = NULL;
247 162 : if( !proj_cs_get_axis_info(
248 : ctx, coordSys, 0,
249 : NULL, /* name */
250 : NULL, /* abbreviation*/
251 : NULL, /* direction */
252 : NULL, /* conversion factor */
253 : NULL, /* unit name */
254 : NULL, /* unit auth name (should be EPSG) */
255 162 : &pszUnitCode) || pszUnitCode == NULL )
256 : {
257 0 : proj_destroy(coordSys);
258 0 : return FALSE;
259 : }
260 162 : *pnUOMLengthCode = (short) atoi(pszUnitCode);
261 162 : proj_destroy(coordSys);
262 : }
263 : }
264 :
265 324 : if( pnGeogCS )
266 : {
267 162 : PJ* geod_crs = proj_crs_get_geodetic_crs(ctx, proj_crs);
268 162 : if( !geod_crs )
269 : {
270 : // shouldn't happen except out of memory
271 0 : proj_destroy(proj_crs);
272 0 : return FALSE;
273 : }
274 :
275 : {
276 162 : const char* pszGeodCode = proj_get_id_code(geod_crs, 0);
277 162 : assert( pszGeodCode );
278 162 : *pnGeogCS = (short) atoi(pszGeodCode);
279 : }
280 :
281 162 : proj_destroy(geod_crs);
282 : }
283 :
284 :
285 324 : proj_destroy(proj_crs);
286 324 : return TRUE;
287 : }
288 : }
289 :
290 :
291 0 : int GTIFGetPCSInfo( int nPCSCode, char **ppszEPSGName,
292 : short *pnProjOp, short *pnUOMLengthCode,
293 : short *pnGeogCS )
294 :
295 : {
296 0 : PJ_CONTEXT* ctx = proj_context_create();
297 0 : int ret = GTIFGetPCSInfoEx(ctx, nPCSCode, ppszEPSGName, pnProjOp,
298 : pnUOMLengthCode, pnGeogCS);
299 0 : proj_context_destroy(ctx);
300 0 : return ret;
301 : }
302 :
303 : /************************************************************************/
304 : /* GTIFAngleToDD() */
305 : /* */
306 : /* Convert a numeric angle to decimal degrees. */
307 : /************************************************************************/
308 :
309 406 : double GTIFAngleToDD( double dfAngle, int nUOMAngle )
310 :
311 : {
312 406 : if( nUOMAngle == 9110 ) /* DDD.MMSSsss */
313 : {
314 0 : if( dfAngle > -999.9 && dfAngle < 999.9 )
315 : {
316 : char szAngleString[32];
317 :
318 0 : snprintf(szAngleString, sizeof(szAngleString), "%12.7f", dfAngle);
319 0 : dfAngle = GTIFAngleStringToDD( szAngleString, nUOMAngle );
320 : }
321 : }
322 406 : else if ( nUOMAngle != KvUserDefined )
323 : {
324 378 : double dfInDegrees = 1.0;
325 :
326 378 : GTIFGetUOMAngleInfo( nUOMAngle, NULL, &dfInDegrees );
327 378 : dfAngle = dfAngle * dfInDegrees;
328 : }
329 :
330 406 : return dfAngle;
331 : }
332 :
333 : /************************************************************************/
334 : /* GTIFAngleStringToDD() */
335 : /* */
336 : /* Convert an angle in the specified units to decimal degrees. */
337 : /************************************************************************/
338 :
339 0 : double GTIFAngleStringToDD( const char * pszAngle, int nUOMAngle )
340 :
341 : {
342 : double dfAngle;
343 :
344 0 : if( nUOMAngle == 9110 ) /* DDD.MMSSsss */
345 : {
346 : char *pszDecimal;
347 :
348 0 : dfAngle = ABS(atoi(pszAngle));
349 0 : pszDecimal = strchr(pszAngle,'.');
350 0 : if( pszDecimal != NULL && strlen(pszDecimal) > 1 )
351 : {
352 : char szMinutes[3];
353 : char szSeconds[64];
354 :
355 0 : szMinutes[0] = pszDecimal[1];
356 0 : if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' )
357 0 : szMinutes[1] = pszDecimal[2];
358 : else
359 0 : szMinutes[1] = '0';
360 :
361 0 : szMinutes[2] = '\0';
362 0 : dfAngle += atoi(szMinutes) / 60.0;
363 :
364 0 : if( strlen(pszDecimal) > 3 )
365 : {
366 0 : szSeconds[0] = pszDecimal[3];
367 0 : if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' )
368 : {
369 0 : szSeconds[1] = pszDecimal[4];
370 0 : szSeconds[2] = '.';
371 0 : strncpy( szSeconds+3, pszDecimal + 5, sizeof(szSeconds) - 3 );
372 0 : szSeconds[sizeof(szSeconds) - 1] = 0;
373 : }
374 : else
375 : {
376 0 : szSeconds[1] = '0';
377 0 : szSeconds[2] = '\0';
378 : }
379 0 : dfAngle += GTIFAtof(szSeconds) / 3600.0;
380 : }
381 : }
382 :
383 0 : if( pszAngle[0] == '-' )
384 0 : dfAngle *= -1;
385 : }
386 0 : else if( nUOMAngle == 9105 || nUOMAngle == 9106 ) /* grad */
387 : {
388 0 : dfAngle = 180 * (GTIFAtof(pszAngle ) / 200);
389 : }
390 0 : else if( nUOMAngle == 9101 ) /* radians */
391 : {
392 0 : dfAngle = 180 * (GTIFAtof(pszAngle ) / M_PI);
393 : }
394 0 : else if( nUOMAngle == 9103 ) /* arc-minute */
395 : {
396 0 : dfAngle = GTIFAtof(pszAngle) / 60;
397 : }
398 0 : else if( nUOMAngle == 9104 ) /* arc-second */
399 : {
400 0 : dfAngle = GTIFAtof(pszAngle) / 3600;
401 : }
402 : else /* decimal degrees ... some cases missing but seemingly never used */
403 : {
404 0 : CPLAssert( nUOMAngle == 9102 || nUOMAngle == KvUserDefined
405 : || nUOMAngle == 0 );
406 :
407 0 : dfAngle = GTIFAtof(pszAngle );
408 : }
409 :
410 0 : return dfAngle;
411 : }
412 :
413 : /************************************************************************/
414 : /* GTIFGetGCSInfo() */
415 : /* */
416 : /* Fetch the datum, and prime meridian related to a particular */
417 : /* GCS. */
418 : /************************************************************************/
419 :
420 10283 : int GTIFGetGCSInfoEx( void* ctxIn,
421 : int nGCSCode, char ** ppszName,
422 : short * pnDatum, short * pnPM, short *pnUOMAngle )
423 :
424 : {
425 10283 : int nDatum=0;
426 :
427 : /* -------------------------------------------------------------------- */
428 : /* Handle some "well known" GCS codes directly */
429 : /* -------------------------------------------------------------------- */
430 10283 : const char * pszName = NULL;
431 10283 : const int nPM = PM_Greenwich;
432 10283 : const int nUOMAngle = Angular_DMS_Hemisphere;
433 10283 : if( nGCSCode == GCS_NAD27 )
434 : {
435 5314 : nDatum = Datum_North_American_Datum_1927;
436 5314 : pszName = "NAD27";
437 : }
438 4969 : else if( nGCSCode == GCS_NAD83 )
439 : {
440 26 : nDatum = Datum_North_American_Datum_1983;
441 26 : pszName = "NAD83";
442 : }
443 4943 : else if( nGCSCode == GCS_WGS_84 )
444 : {
445 4670 : nDatum = Datum_WGS84;
446 4670 : pszName = "WGS 84";
447 : }
448 273 : else if( nGCSCode == GCS_WGS_72 )
449 : {
450 10 : nDatum = Datum_WGS72;
451 10 : pszName = "WGS 72";
452 : }
453 263 : else if ( nGCSCode == KvUserDefined )
454 : {
455 73 : return FALSE;
456 : }
457 :
458 10210 : if (pszName != NULL)
459 : {
460 10020 : if( ppszName != NULL )
461 5006 : *ppszName = CPLStrdup( pszName );
462 10020 : if( pnDatum != NULL )
463 5014 : *pnDatum = (short) nDatum;
464 10020 : if( pnPM != NULL )
465 5014 : *pnPM = (short) nPM;
466 10020 : if( pnUOMAngle != NULL )
467 5014 : *pnUOMAngle = (short) nUOMAngle;
468 :
469 10020 : return TRUE;
470 : }
471 :
472 : /* -------------------------------------------------------------------- */
473 : /* Search the database. */
474 : /* -------------------------------------------------------------------- */
475 :
476 190 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
477 : {
478 : char szCode[12];
479 :
480 190 : snprintf(szCode, sizeof(szCode), "%d", nGCSCode);
481 190 : PJ* geod_crs = proj_create_from_database(
482 : ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
483 190 : if( !geod_crs )
484 : {
485 0 : return FALSE;
486 : }
487 :
488 : {
489 190 : const int objType = proj_get_type(geod_crs);
490 190 : if( objType != PJ_TYPE_GEODETIC_CRS &&
491 190 : objType != PJ_TYPE_GEOCENTRIC_CRS &&
492 2 : objType != PJ_TYPE_GEOGRAPHIC_2D_CRS &&
493 : objType != PJ_TYPE_GEOGRAPHIC_3D_CRS )
494 : {
495 0 : proj_destroy(geod_crs);
496 0 : return FALSE;
497 : }
498 : }
499 :
500 190 : if( ppszName )
501 : {
502 95 : pszName = proj_get_name(geod_crs);
503 95 : if( !pszName )
504 : {
505 : // shouldn't happen
506 0 : proj_destroy(geod_crs);
507 0 : return FALSE;
508 : }
509 95 : *ppszName = CPLStrdup(pszName);
510 : }
511 :
512 190 : if( pnDatum )
513 : {
514 : #if PROJ_VERSION_MAJOR >= 8
515 : PJ* datum = proj_crs_get_datum_forced(ctx, geod_crs);
516 : #else
517 95 : PJ* datum = proj_crs_get_datum(ctx, geod_crs);
518 : #endif
519 95 : if( !datum )
520 : {
521 0 : proj_destroy(geod_crs);
522 0 : return FALSE;
523 : }
524 :
525 : {
526 95 : const char* pszDatumCode = proj_get_id_code(datum, 0);
527 95 : assert( pszDatumCode );
528 95 : *pnDatum = (short) atoi(pszDatumCode);
529 : }
530 :
531 95 : proj_destroy(datum);
532 : }
533 :
534 190 : if( pnPM )
535 : {
536 95 : PJ* pm = proj_get_prime_meridian(ctx, geod_crs);
537 95 : if( !pm )
538 : {
539 0 : proj_destroy(geod_crs);
540 0 : return FALSE;
541 : }
542 :
543 : {
544 95 : const char* pszPMCode = proj_get_id_code(pm, 0);
545 95 : assert( pszPMCode );
546 95 : *pnPM = (short) atoi(pszPMCode);
547 : }
548 :
549 95 : proj_destroy(pm);
550 : }
551 :
552 190 : if( pnUOMAngle )
553 : {
554 95 : PJ* coordSys = proj_crs_get_coordinate_system(
555 : ctx, geod_crs);
556 95 : if( !coordSys )
557 : {
558 : // shouldn't happen except out of memory
559 0 : proj_destroy(geod_crs);
560 0 : return FALSE;
561 : }
562 :
563 : {
564 95 : const char* pszUnitCode = NULL;
565 95 : if( !proj_cs_get_axis_info(
566 : ctx, coordSys, 0,
567 : NULL, /* name */
568 : NULL, /* abbreviation*/
569 : NULL, /* direction */
570 : NULL, /* conversion factor */
571 : NULL, /* unit name */
572 : NULL, /* unit auth name (should be EPSG) */
573 95 : &pszUnitCode) || pszUnitCode == NULL )
574 : {
575 0 : proj_destroy(coordSys);
576 0 : return FALSE;
577 : }
578 95 : *pnUOMAngle = (short) atoi(pszUnitCode);
579 95 : proj_destroy(coordSys);
580 : }
581 : }
582 :
583 190 : proj_destroy(geod_crs);
584 190 : return TRUE;
585 : }
586 : }
587 :
588 20 : int GTIFGetGCSInfo( int nGCSCode, char ** ppszName,
589 : short * pnDatum, short * pnPM, short *pnUOMAngle )
590 :
591 : {
592 20 : PJ_CONTEXT* ctx = proj_context_create();
593 20 : const int ret = GTIFGetGCSInfoEx(ctx, nGCSCode, ppszName, pnDatum,
594 : pnPM, pnUOMAngle);
595 20 : proj_context_destroy(ctx);
596 20 : return ret;
597 : }
598 :
599 : /************************************************************************/
600 : /* GTIFGetEllipsoidInfo() */
601 : /* */
602 : /* Fetch info about an ellipsoid. Axes are always returned in */
603 : /* meters. SemiMajor computed based on inverse flattening */
604 : /* where that is provided. */
605 : /************************************************************************/
606 :
607 10226 : int GTIFGetEllipsoidInfoEx( void* ctxIn,
608 : int nEllipseCode, char ** ppszName,
609 : double * pdfSemiMajor, double * pdfSemiMinor )
610 :
611 : {
612 10226 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
613 : /* -------------------------------------------------------------------- */
614 : /* Try some well known ellipsoids. */
615 : /* -------------------------------------------------------------------- */
616 10226 : double dfSemiMajor=0.0;
617 10226 : double dfInvFlattening=0.0, dfSemiMinor=0.0;
618 10226 : const char *pszName = NULL;
619 :
620 10226 : if( nEllipseCode == Ellipse_Clarke_1866 )
621 : {
622 5316 : pszName = "Clarke 1866";
623 5316 : dfSemiMajor = 6378206.4;
624 5316 : dfSemiMinor = 6356583.8;
625 5316 : dfInvFlattening = 0.0;
626 : }
627 4910 : else if( nEllipseCode == Ellipse_GRS_1980 )
628 : {
629 135 : pszName = "GRS 1980";
630 135 : dfSemiMajor = 6378137.0;
631 135 : dfSemiMinor = 0.0;
632 135 : dfInvFlattening = 298.257222101;
633 : }
634 4775 : else if( nEllipseCode == Ellipse_WGS_84 )
635 : {
636 4669 : pszName = "WGS 84";
637 4669 : dfSemiMajor = 6378137.0;
638 4669 : dfSemiMinor = 0.0;
639 4669 : dfInvFlattening = 298.257223563;
640 : }
641 106 : else if( nEllipseCode == 7043 )
642 : {
643 10 : pszName = "WGS 72";
644 10 : dfSemiMajor = 6378135.0;
645 10 : dfSemiMinor = 0.0;
646 10 : dfInvFlattening = 298.26;
647 : }
648 :
649 10226 : if (pszName != NULL)
650 : {
651 10130 : if( dfSemiMinor == 0.0 )
652 4814 : dfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening);
653 :
654 10130 : if( pdfSemiMinor != NULL )
655 5069 : *pdfSemiMinor = dfSemiMinor;
656 10130 : if( pdfSemiMajor != NULL )
657 5069 : *pdfSemiMajor = dfSemiMajor;
658 10130 : if( ppszName != NULL )
659 5061 : *ppszName = CPLStrdup( pszName );
660 :
661 10130 : return TRUE;
662 : }
663 :
664 96 : if( nEllipseCode == KvUserDefined )
665 8 : return FALSE;
666 :
667 : /* -------------------------------------------------------------------- */
668 : /* Search the database. */
669 : /* -------------------------------------------------------------------- */
670 : {
671 : char szCode[12];
672 :
673 88 : snprintf(szCode, sizeof(szCode), "%d", nEllipseCode);
674 88 : PJ* ellipsoid = proj_create_from_database(
675 : ctx, "EPSG", szCode, PJ_CATEGORY_ELLIPSOID, 0, NULL);
676 88 : if( !ellipsoid )
677 : {
678 3 : return FALSE;
679 : }
680 :
681 85 : if( ppszName )
682 : {
683 42 : pszName = proj_get_name(ellipsoid);
684 42 : if( !pszName )
685 : {
686 : // shouldn't happen
687 0 : proj_destroy(ellipsoid);
688 0 : return FALSE;
689 : }
690 42 : *ppszName = CPLStrdup(pszName);
691 : }
692 :
693 85 : proj_ellipsoid_get_parameters(
694 : ctx, ellipsoid, pdfSemiMajor, pdfSemiMinor, NULL, NULL);
695 :
696 85 : proj_destroy(ellipsoid);
697 :
698 85 : return TRUE;
699 : }
700 : }
701 :
702 20 : int GTIFGetEllipsoidInfo( int nEllipseCode, char ** ppszName,
703 : double * pdfSemiMajor, double * pdfSemiMinor )
704 :
705 : {
706 20 : PJ_CONTEXT* ctx = proj_context_create();
707 20 : const int ret = GTIFGetEllipsoidInfoEx(ctx, nEllipseCode, ppszName, pdfSemiMajor,
708 : pdfSemiMinor);
709 20 : proj_context_destroy(ctx);
710 20 : return ret;
711 : }
712 :
713 : /************************************************************************/
714 : /* GTIFGetPMInfo() */
715 : /* */
716 : /* Get the offset between a given prime meridian and Greenwich */
717 : /* in degrees. */
718 : /************************************************************************/
719 :
720 10219 : int GTIFGetPMInfoEx( void* ctxIn,
721 : int nPMCode, char ** ppszName, double *pdfOffset )
722 :
723 : {
724 : /* -------------------------------------------------------------------- */
725 : /* Use a special short cut for Greenwich, since it is so common. */
726 : /* -------------------------------------------------------------------- */
727 10219 : if( nPMCode == PM_Greenwich )
728 : {
729 10198 : if( pdfOffset != NULL )
730 5103 : *pdfOffset = 0.0;
731 10198 : if( ppszName != NULL )
732 5095 : *ppszName = CPLStrdup( "Greenwich" );
733 10198 : return TRUE;
734 : }
735 :
736 :
737 21 : if( nPMCode == KvUserDefined )
738 7 : return FALSE;
739 :
740 : /* -------------------------------------------------------------------- */
741 : /* Search the database. */
742 : /* -------------------------------------------------------------------- */
743 14 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
744 : {
745 : char szCode[12];
746 :
747 14 : snprintf(szCode, sizeof(szCode), "%d", nPMCode);
748 14 : PJ* pm = proj_create_from_database(
749 : ctx, "EPSG", szCode, PJ_CATEGORY_PRIME_MERIDIAN, 0, NULL);
750 14 : if( !pm )
751 : {
752 0 : return FALSE;
753 : }
754 :
755 14 : if( ppszName )
756 : {
757 7 : const char* pszName = proj_get_name(pm);
758 7 : if( !pszName )
759 : {
760 : // shouldn't happen
761 0 : proj_destroy(pm);
762 0 : return FALSE;
763 : }
764 7 : *ppszName = CPLStrdup(pszName);
765 : }
766 :
767 14 : if( pdfOffset )
768 : {
769 7 : double conv_factor = 0;
770 7 : proj_prime_meridian_get_parameters(
771 : ctx, pm, pdfOffset, &conv_factor, NULL);
772 7 : *pdfOffset *= conv_factor * 180.0 / M_PI;
773 : }
774 :
775 14 : proj_destroy(pm);
776 :
777 14 : return TRUE;
778 : }
779 : }
780 :
781 20 : int GTIFGetPMInfo( int nPMCode, char ** ppszName, double *pdfOffset )
782 :
783 : {
784 20 : PJ_CONTEXT* ctx = proj_context_create();
785 20 : const int ret = GTIFGetPMInfoEx(ctx, nPMCode, ppszName, pdfOffset);
786 20 : proj_context_destroy(ctx);
787 20 : return ret;
788 : }
789 :
790 : /************************************************************************/
791 : /* GTIFGetDatumInfo() */
792 : /* */
793 : /* Fetch the ellipsoid, and name for a datum. */
794 : /************************************************************************/
795 :
796 10222 : int GTIFGetDatumInfoEx( void* ctxIn,
797 : int nDatumCode, char ** ppszName, short * pnEllipsoid )
798 :
799 : {
800 10222 : const char* pszName = NULL;
801 10222 : int nEllipsoid = 0;
802 :
803 : /* -------------------------------------------------------------------- */
804 : /* Handle a few built-in datums. */
805 : /* -------------------------------------------------------------------- */
806 10222 : if( nDatumCode == Datum_North_American_Datum_1927 )
807 : {
808 5314 : nEllipsoid = Ellipse_Clarke_1866;
809 5314 : pszName = "North American Datum 1927";
810 : }
811 4908 : else if( nDatumCode == Datum_North_American_Datum_1983 )
812 : {
813 26 : nEllipsoid = Ellipse_GRS_1980;
814 26 : pszName = "North American Datum 1983";
815 : }
816 4882 : else if( nDatumCode == Datum_WGS84 )
817 : {
818 4668 : nEllipsoid = Ellipse_WGS_84;
819 4668 : pszName = "World Geodetic System 1984";
820 : }
821 214 : else if( nDatumCode == Datum_WGS72 )
822 : {
823 10 : nEllipsoid = 7043; /* WGS72 */
824 10 : pszName = "World Geodetic System 1972";
825 : }
826 :
827 10222 : if (pszName != NULL)
828 : {
829 10018 : if( pnEllipsoid != NULL )
830 5013 : *pnEllipsoid = (short) nEllipsoid;
831 :
832 10018 : if( ppszName != NULL )
833 5005 : *ppszName = CPLStrdup( pszName );
834 :
835 10018 : return TRUE;
836 : }
837 :
838 204 : if( nDatumCode == KvUserDefined )
839 8 : return FALSE;
840 :
841 : /* -------------------------------------------------------------------- */
842 : /* Search the database. */
843 : /* -------------------------------------------------------------------- */
844 196 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
845 : {
846 : char szCode[12];
847 :
848 196 : snprintf(szCode, sizeof(szCode), "%d", nDatumCode);
849 196 : PJ* datum = proj_create_from_database(
850 : ctx, "EPSG", szCode, PJ_CATEGORY_DATUM, 0, NULL);
851 196 : if( !datum )
852 : {
853 2 : return FALSE;
854 : }
855 :
856 194 : const PJ_TYPE pjType = proj_get_type(datum);
857 194 : if( pjType != PJ_TYPE_GEODETIC_REFERENCE_FRAME &&
858 : pjType != PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME )
859 : {
860 0 : proj_destroy(datum);
861 0 : return FALSE;
862 : }
863 :
864 194 : if( ppszName )
865 : {
866 97 : pszName = proj_get_name(datum);
867 97 : if( !pszName )
868 : {
869 : // shouldn't happen
870 0 : proj_destroy(datum);
871 0 : return FALSE;
872 : }
873 97 : *ppszName = CPLStrdup(pszName);
874 : }
875 :
876 194 : if( pnEllipsoid )
877 : {
878 97 : PJ* ellipsoid = proj_get_ellipsoid(ctx, datum);
879 97 : if( !ellipsoid )
880 : {
881 0 : proj_destroy(datum);
882 0 : return FALSE;
883 : }
884 :
885 : {
886 97 : const char* pszEllipsoidCode = proj_get_id_code(
887 : ellipsoid, 0);
888 97 : assert( pszEllipsoidCode );
889 97 : *pnEllipsoid = (short) atoi(pszEllipsoidCode);
890 : }
891 :
892 97 : proj_destroy(ellipsoid);
893 : }
894 :
895 194 : proj_destroy(datum);
896 :
897 194 : return TRUE;
898 : }
899 : }
900 :
901 20 : int GTIFGetDatumInfo( int nDatumCode, char ** ppszName, short * pnEllipsoid )
902 :
903 : {
904 20 : PJ_CONTEXT* ctx = proj_context_create();
905 20 : const int ret = GTIFGetDatumInfoEx(ctx, nDatumCode, ppszName, pnEllipsoid);
906 20 : proj_context_destroy(ctx);
907 20 : return ret;
908 : }
909 :
910 : /************************************************************************/
911 : /* GTIFGetUOMLengthInfo() */
912 : /* */
913 : /* Note: This function should eventually also know how to */
914 : /* lookup length aliases in the UOM_LE_ALIAS table. */
915 : /************************************************************************/
916 :
917 6810 : int GTIFGetUOMLengthInfoEx( void* ctxIn,
918 : int nUOMLengthCode,
919 : char **ppszUOMName,
920 : double * pdfInMeters )
921 :
922 : {
923 : /* -------------------------------------------------------------------- */
924 : /* We short cut meter to save work and avoid failure for missing */
925 : /* in the most common cases. */
926 : /* -------------------------------------------------------------------- */
927 6810 : if( nUOMLengthCode == 9001 )
928 : {
929 6757 : if( ppszUOMName != NULL )
930 3406 : *ppszUOMName = CPLStrdup( "metre" );
931 6757 : if( pdfInMeters != NULL )
932 3351 : *pdfInMeters = 1.0;
933 :
934 6757 : return TRUE;
935 : }
936 :
937 53 : if( nUOMLengthCode == 9002 )
938 : {
939 0 : if( ppszUOMName != NULL )
940 0 : *ppszUOMName = CPLStrdup( "foot" );
941 0 : if( pdfInMeters != NULL )
942 0 : *pdfInMeters = 0.3048;
943 :
944 0 : return TRUE;
945 : }
946 :
947 53 : if( nUOMLengthCode == 9003 )
948 : {
949 46 : if( ppszUOMName != NULL )
950 28 : *ppszUOMName = CPLStrdup( "US survey foot" );
951 46 : if( pdfInMeters != NULL )
952 18 : *pdfInMeters = 12.0 / 39.37;
953 :
954 46 : return TRUE;
955 : }
956 :
957 7 : if( nUOMLengthCode == KvUserDefined )
958 5 : return FALSE;
959 :
960 : /* -------------------------------------------------------------------- */
961 : /* Search the units database for this unit. If we don't find */
962 : /* it return failure. */
963 : /* -------------------------------------------------------------------- */
964 : char szCode[12];
965 2 : const char* pszName = NULL;
966 :
967 2 : snprintf(szCode, sizeof(szCode), "%d", nUOMLengthCode);
968 2 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
969 2 : if( !proj_uom_get_info_from_database(
970 : ctx, "EPSG", szCode, &pszName, pdfInMeters, NULL) )
971 : {
972 0 : return FALSE;
973 : }
974 2 : if( ppszUOMName )
975 : {
976 1 : *ppszUOMName = CPLStrdup(pszName);
977 : }
978 2 : return TRUE;
979 : }
980 :
981 105 : int GTIFGetUOMLengthInfo( int nUOMLengthCode,
982 : char **ppszUOMName,
983 : double * pdfInMeters )
984 :
985 : {
986 105 : PJ_CONTEXT* ctx = proj_context_create();
987 105 : const int ret = GTIFGetUOMLengthInfoEx(
988 : ctx, nUOMLengthCode, ppszUOMName, pdfInMeters);
989 105 : proj_context_destroy(ctx);
990 105 : return ret;
991 : }
992 :
993 : /************************************************************************/
994 : /* GTIFGetUOMAngleInfo() */
995 : /************************************************************************/
996 :
997 11026 : int GTIFGetUOMAngleInfoEx( void* ctxIn,
998 : int nUOMAngleCode,
999 : char **ppszUOMName,
1000 : double * pdfInDegrees )
1001 :
1002 : {
1003 11026 : const char *pszUOMName = NULL;
1004 11026 : double dfInDegrees = 1.0;
1005 :
1006 11026 : switch( nUOMAngleCode )
1007 : {
1008 4 : case 9101:
1009 4 : pszUOMName = "radian";
1010 4 : dfInDegrees = 180.0 / M_PI;
1011 4 : break;
1012 :
1013 10989 : case 9102:
1014 : case 9107:
1015 : case 9108:
1016 : case 9110:
1017 : case 9122:
1018 10989 : pszUOMName = "degree";
1019 10989 : dfInDegrees = 1.0;
1020 10989 : break;
1021 :
1022 4 : case 9103:
1023 4 : pszUOMName = "arc-minute";
1024 4 : dfInDegrees = 1 / 60.0;
1025 4 : break;
1026 :
1027 4 : case 9104:
1028 4 : pszUOMName = "arc-second";
1029 4 : dfInDegrees = 1 / 3600.0;
1030 4 : break;
1031 :
1032 14 : case 9105:
1033 14 : pszUOMName = "grad";
1034 14 : dfInDegrees = 180.0 / 200.0;
1035 14 : break;
1036 :
1037 4 : case 9106:
1038 4 : pszUOMName = "gon";
1039 4 : dfInDegrees = 180.0 / 200.0;
1040 4 : break;
1041 :
1042 0 : case 9109:
1043 0 : pszUOMName = "microradian";
1044 0 : dfInDegrees = 180.0 / (M_PI * 1000000.0);
1045 0 : break;
1046 :
1047 7 : default:
1048 7 : break;
1049 : }
1050 :
1051 11026 : if (pszUOMName)
1052 : {
1053 11019 : if( ppszUOMName != NULL )
1054 : {
1055 5153 : *ppszUOMName = CPLStrdup( pszUOMName );
1056 : }
1057 :
1058 11019 : if( pdfInDegrees != NULL )
1059 11009 : *pdfInDegrees = dfInDegrees;
1060 :
1061 11019 : return TRUE;
1062 : }
1063 :
1064 7 : if( nUOMAngleCode == KvUserDefined )
1065 7 : return FALSE;
1066 :
1067 : /* -------------------------------------------------------------------- */
1068 : /* Search the units database for this unit. If we don't find */
1069 : /* it return failure. */
1070 : /* -------------------------------------------------------------------- */
1071 : char szCode[12];
1072 0 : const char* pszName = NULL;
1073 0 : double dfConvFactorToRadians = 0;
1074 :
1075 0 : snprintf(szCode, sizeof(szCode), "%d", nUOMAngleCode);
1076 0 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
1077 0 : if( !proj_uom_get_info_from_database(
1078 : ctx, "EPSG", szCode, &pszName, &dfConvFactorToRadians, NULL) )
1079 : {
1080 0 : return FALSE;
1081 : }
1082 0 : if( ppszUOMName )
1083 : {
1084 0 : *ppszUOMName = CPLStrdup(pszName);
1085 : }
1086 0 : if( pdfInDegrees )
1087 : {
1088 0 : *pdfInDegrees = dfConvFactorToRadians * 180.0 / M_PI;
1089 : }
1090 0 : return TRUE;
1091 : }
1092 :
1093 398 : int GTIFGetUOMAngleInfo( int nUOMAngleCode,
1094 : char **ppszUOMName,
1095 : double * pdfInDegrees )
1096 :
1097 : {
1098 398 : PJ_CONTEXT* ctx = proj_context_create();
1099 398 : const int ret = GTIFGetUOMAngleInfoEx(
1100 : ctx, nUOMAngleCode, ppszUOMName, pdfInDegrees);
1101 398 : proj_context_destroy(ctx);
1102 398 : return ret;
1103 : }
1104 :
1105 : /************************************************************************/
1106 : /* EPSGProjMethodToCTProjMethod() */
1107 : /* */
1108 : /* Convert between the EPSG enumeration for projection methods, */
1109 : /* and the GeoTIFF CT codes. */
1110 : /************************************************************************/
1111 :
1112 6070 : static int EPSGProjMethodToCTProjMethod( int nEPSG, int bReturnExtendedCTCode )
1113 :
1114 : {
1115 6070 : switch( nEPSG )
1116 : {
1117 15 : case 9801:
1118 15 : return CT_LambertConfConic_1SP;
1119 :
1120 12 : case 9802:
1121 12 : return CT_LambertConfConic_2SP;
1122 :
1123 0 : case 9803:
1124 0 : return CT_LambertConfConic_2SP; /* Belgian variant not supported */
1125 :
1126 3 : case 9804:
1127 3 : return CT_Mercator; /* 1SP and 2SP not differentiated */
1128 :
1129 3 : case 9805:
1130 3 : if( bReturnExtendedCTCode )
1131 2 : return CT_Ext_Mercator_2SP;
1132 : else
1133 1 : return CT_Mercator; /* 1SP and 2SP not differentiated */
1134 :
1135 : /* Mercator 1SP (Spherical) For EPSG:3785 */
1136 0 : case 9841:
1137 0 : return CT_Mercator; /* 1SP and 2SP not differentiated */
1138 :
1139 : /* Google Mercator For EPSG:3857 */
1140 171 : case 1024:
1141 171 : return CT_Mercator; /* 1SP and 2SP not differentiated */
1142 :
1143 3 : case 9806:
1144 3 : return CT_CassiniSoldner;
1145 :
1146 5788 : case 9807:
1147 5788 : return CT_TransverseMercator;
1148 :
1149 0 : case 9808:
1150 0 : return CT_TransvMercator_SouthOriented;
1151 :
1152 12 : case 9809:
1153 12 : return CT_ObliqueStereographic;
1154 :
1155 18 : case 9810:
1156 : case 9829: /* variant B not quite the same - not sure how to handle */
1157 18 : return CT_PolarStereographic;
1158 :
1159 3 : case 9811:
1160 3 : return CT_NewZealandMapGrid;
1161 :
1162 3 : case 9812:
1163 3 : return CT_ObliqueMercator; /* is hotine actually different? */
1164 :
1165 0 : case 9813:
1166 0 : return CT_ObliqueMercator_Laborde;
1167 :
1168 0 : case 9814:
1169 0 : return CT_ObliqueMercator_Rosenmund; /* swiss */
1170 :
1171 6 : case 9815:
1172 6 : return CT_HotineObliqueMercatorAzimuthCenter;
1173 :
1174 0 : case 9816: /* tunesia mining grid has no counterpart */
1175 0 : return KvUserDefined;
1176 :
1177 3 : case 9818:
1178 3 : return CT_Polyconic;
1179 :
1180 3 : case 9820:
1181 : case 1027:
1182 3 : return CT_LambertAzimEqualArea;
1183 :
1184 3 : case 9822:
1185 3 : return CT_AlbersEqualArea;
1186 :
1187 0 : case 9834:
1188 0 : return CT_CylindricalEqualArea;
1189 :
1190 9 : case 1028:
1191 : case 1029:
1192 : case 9823: /* spherical */
1193 : case 9842: /* elliptical */
1194 9 : return CT_Equirectangular;
1195 :
1196 15 : default: /* use the EPSG code for other methods */
1197 15 : return nEPSG;
1198 : }
1199 : }
1200 :
1201 : /************************************************************************/
1202 : /* SetGTParamIds() */
1203 : /* */
1204 : /* This is hardcoded logic to set the GeoTIFF parameter */
1205 : /* identifiers for all the EPSG supported projections. As new */
1206 : /* projection methods are added, this code will need to be updated */
1207 : /************************************************************************/
1208 :
1209 3110 : static int SetGTParamIds( int nCTProjection,
1210 : int nEPSGProjMethod,
1211 : int *panProjParamId,
1212 : int *panEPSGCodes )
1213 :
1214 : {
1215 : int anWorkingDummy[7];
1216 :
1217 3110 : if( panEPSGCodes == NULL )
1218 2960 : panEPSGCodes = anWorkingDummy;
1219 3110 : if( panProjParamId == NULL )
1220 150 : panProjParamId = anWorkingDummy;
1221 :
1222 3110 : memset( panEPSGCodes, 0, sizeof(int) * 7 );
1223 :
1224 : /* psDefn->nParms = 7; */
1225 :
1226 3110 : switch( nCTProjection )
1227 : {
1228 6 : case CT_CassiniSoldner:
1229 : case CT_NewZealandMapGrid:
1230 : case CT_Polyconic:
1231 6 : panProjParamId[0] = ProjNatOriginLatGeoKey;
1232 6 : panProjParamId[1] = ProjNatOriginLongGeoKey;
1233 6 : panProjParamId[5] = ProjFalseEastingGeoKey;
1234 6 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1235 :
1236 6 : panEPSGCodes[0] = EPSGNatOriginLat;
1237 6 : panEPSGCodes[1] = EPSGNatOriginLong;
1238 6 : panEPSGCodes[5] = EPSGFalseEasting;
1239 6 : panEPSGCodes[6] = EPSGFalseNorthing;
1240 6 : return TRUE;
1241 :
1242 6 : case CT_ObliqueMercator:
1243 : case CT_HotineObliqueMercatorAzimuthCenter:
1244 6 : panProjParamId[0] = ProjCenterLatGeoKey;
1245 6 : panProjParamId[1] = ProjCenterLongGeoKey;
1246 6 : panProjParamId[2] = ProjAzimuthAngleGeoKey;
1247 6 : panProjParamId[3] = ProjRectifiedGridAngleGeoKey;
1248 6 : panProjParamId[4] = ProjScaleAtCenterGeoKey;
1249 6 : panProjParamId[5] = ProjFalseEastingGeoKey;
1250 6 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1251 :
1252 6 : panEPSGCodes[0] = EPSGProjCenterLat;
1253 6 : panEPSGCodes[1] = EPSGProjCenterLong;
1254 6 : panEPSGCodes[2] = EPSGAzimuth;
1255 6 : panEPSGCodes[3] = EPSGAngleRectifiedToSkewedGrid;
1256 6 : panEPSGCodes[4] = EPSGInitialLineScaleFactor;
1257 6 : panEPSGCodes[5] = EPSGProjCenterEasting; /* EPSG proj method 9812 uses EPSGFalseEasting, but 9815 uses EPSGProjCenterEasting */
1258 6 : panEPSGCodes[6] = EPSGProjCenterNorthing; /* EPSG proj method 9812 uses EPSGFalseNorthing, but 9815 uses EPSGProjCenterNorthing */
1259 6 : return TRUE;
1260 :
1261 0 : case CT_ObliqueMercator_Laborde:
1262 0 : panProjParamId[0] = ProjCenterLatGeoKey;
1263 0 : panProjParamId[1] = ProjCenterLongGeoKey;
1264 0 : panProjParamId[2] = ProjAzimuthAngleGeoKey;
1265 0 : panProjParamId[4] = ProjScaleAtCenterGeoKey;
1266 0 : panProjParamId[5] = ProjFalseEastingGeoKey;
1267 0 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1268 :
1269 0 : panEPSGCodes[0] = EPSGProjCenterLat;
1270 0 : panEPSGCodes[1] = EPSGProjCenterLong;
1271 0 : panEPSGCodes[2] = EPSGAzimuth;
1272 0 : panEPSGCodes[4] = EPSGInitialLineScaleFactor;
1273 0 : panEPSGCodes[5] = EPSGFalseEasting;
1274 0 : panEPSGCodes[6] = EPSGFalseNorthing;
1275 0 : return TRUE;
1276 :
1277 3068 : case CT_LambertConfConic_1SP:
1278 : case CT_Mercator:
1279 : case CT_ObliqueStereographic:
1280 : case CT_PolarStereographic:
1281 : case CT_TransverseMercator:
1282 : case CT_TransvMercator_SouthOriented:
1283 3068 : panProjParamId[0] = ProjNatOriginLatGeoKey;
1284 3068 : if( nCTProjection == CT_PolarStereographic )
1285 : {
1286 12 : panProjParamId[1] = ProjStraightVertPoleLongGeoKey;
1287 : }
1288 : else
1289 : {
1290 3056 : panProjParamId[1] = ProjNatOriginLongGeoKey;
1291 : }
1292 3068 : if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
1293 : {
1294 0 : panProjParamId[2] = ProjStdParallel1GeoKey;
1295 : }
1296 3068 : panProjParamId[4] = ProjScaleAtNatOriginGeoKey;
1297 3068 : panProjParamId[5] = ProjFalseEastingGeoKey;
1298 3068 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1299 :
1300 3068 : panEPSGCodes[0] = EPSGNatOriginLat;
1301 3068 : panEPSGCodes[1] = EPSGNatOriginLong;
1302 3068 : if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
1303 : {
1304 0 : panEPSGCodes[2] = EPSGStdParallel1Lat;
1305 : }
1306 3068 : panEPSGCodes[4] = EPSGNatOriginScaleFactor;
1307 3068 : panEPSGCodes[5] = EPSGFalseEasting;
1308 3068 : panEPSGCodes[6] = EPSGFalseNorthing;
1309 3068 : return TRUE;
1310 :
1311 8 : case CT_LambertConfConic_2SP:
1312 8 : panProjParamId[0] = ProjFalseOriginLatGeoKey;
1313 8 : panProjParamId[1] = ProjFalseOriginLongGeoKey;
1314 8 : panProjParamId[2] = ProjStdParallel1GeoKey;
1315 8 : panProjParamId[3] = ProjStdParallel2GeoKey;
1316 8 : panProjParamId[5] = ProjFalseEastingGeoKey;
1317 8 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1318 :
1319 8 : panEPSGCodes[0] = EPSGFalseOriginLat;
1320 8 : panEPSGCodes[1] = EPSGFalseOriginLong;
1321 8 : panEPSGCodes[2] = EPSGStdParallel1Lat;
1322 8 : panEPSGCodes[3] = EPSGStdParallel2Lat;
1323 8 : panEPSGCodes[5] = EPSGFalseOriginEasting;
1324 8 : panEPSGCodes[6] = EPSGFalseOriginNorthing;
1325 8 : return TRUE;
1326 :
1327 2 : case CT_AlbersEqualArea:
1328 2 : panProjParamId[0] = ProjStdParallel1GeoKey;
1329 2 : panProjParamId[1] = ProjStdParallel2GeoKey;
1330 2 : panProjParamId[2] = ProjNatOriginLatGeoKey;
1331 2 : panProjParamId[3] = ProjNatOriginLongGeoKey;
1332 2 : panProjParamId[5] = ProjFalseEastingGeoKey;
1333 2 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1334 :
1335 2 : panEPSGCodes[0] = EPSGStdParallel1Lat;
1336 2 : panEPSGCodes[1] = EPSGStdParallel2Lat;
1337 2 : panEPSGCodes[2] = EPSGFalseOriginLat;
1338 2 : panEPSGCodes[3] = EPSGFalseOriginLong;
1339 2 : panEPSGCodes[5] = EPSGFalseOriginEasting;
1340 2 : panEPSGCodes[6] = EPSGFalseOriginNorthing;
1341 2 : return TRUE;
1342 :
1343 0 : case CT_SwissObliqueCylindrical:
1344 0 : panProjParamId[0] = ProjCenterLatGeoKey;
1345 0 : panProjParamId[1] = ProjCenterLongGeoKey;
1346 0 : panProjParamId[5] = ProjFalseEastingGeoKey;
1347 0 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1348 :
1349 : /* EPSG codes? */
1350 0 : return TRUE;
1351 :
1352 2 : case CT_LambertAzimEqualArea:
1353 2 : panProjParamId[0] = ProjCenterLatGeoKey;
1354 2 : panProjParamId[1] = ProjCenterLongGeoKey;
1355 2 : panProjParamId[5] = ProjFalseEastingGeoKey;
1356 2 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1357 :
1358 2 : panEPSGCodes[0] = EPSGNatOriginLat;
1359 2 : panEPSGCodes[1] = EPSGNatOriginLong;
1360 2 : panEPSGCodes[5] = EPSGFalseEasting;
1361 2 : panEPSGCodes[6] = EPSGFalseNorthing;
1362 2 : return TRUE;
1363 :
1364 0 : case CT_CylindricalEqualArea:
1365 0 : panProjParamId[0] = ProjStdParallel1GeoKey;
1366 0 : panProjParamId[1] = ProjNatOriginLongGeoKey;
1367 0 : panProjParamId[5] = ProjFalseEastingGeoKey;
1368 0 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1369 :
1370 0 : panEPSGCodes[0] = EPSGStdParallel1Lat;
1371 0 : panEPSGCodes[1] = EPSGFalseOriginLong;
1372 0 : panEPSGCodes[5] = EPSGFalseOriginEasting;
1373 0 : panEPSGCodes[6] = EPSGFalseOriginNorthing;
1374 0 : return TRUE;
1375 :
1376 6 : case CT_Equirectangular:
1377 6 : panProjParamId[0] = ProjCenterLatGeoKey;
1378 6 : panProjParamId[1] = ProjCenterLongGeoKey;
1379 6 : panProjParamId[2] = ProjStdParallel1GeoKey;
1380 6 : panProjParamId[5] = ProjFalseEastingGeoKey;
1381 6 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1382 :
1383 6 : panEPSGCodes[0] = EPSGNatOriginLat;
1384 6 : panEPSGCodes[1] = EPSGNatOriginLong;
1385 6 : panEPSGCodes[2] = EPSGStdParallel1Lat;
1386 6 : panEPSGCodes[5] = EPSGFalseEasting;
1387 6 : panEPSGCodes[6] = EPSGFalseNorthing;
1388 6 : return TRUE;
1389 :
1390 2 : case CT_Ext_Mercator_2SP:
1391 2 : panProjParamId[0] = ProjNatOriginLatGeoKey;
1392 2 : panProjParamId[1] = ProjNatOriginLongGeoKey;
1393 2 : panProjParamId[2] = ProjStdParallel1GeoKey;
1394 2 : panProjParamId[5] = ProjFalseEastingGeoKey;
1395 2 : panProjParamId[6] = ProjFalseNorthingGeoKey;
1396 :
1397 2 : panEPSGCodes[0] = EPSGNatOriginLat;
1398 2 : panEPSGCodes[1] = EPSGNatOriginLong;
1399 2 : panEPSGCodes[2] = EPSGStdParallel1Lat;
1400 2 : panEPSGCodes[5] = EPSGFalseEasting;
1401 2 : panEPSGCodes[6] = EPSGFalseNorthing;
1402 2 : return TRUE;
1403 :
1404 10 : default:
1405 10 : return FALSE;
1406 : }
1407 : }
1408 :
1409 : /************************************************************************/
1410 : /* GTIFGetProjTRFInfo() */
1411 : /* */
1412 : /* Transform a PROJECTION_TRF_CODE into a projection method, */
1413 : /* and a set of parameters. The parameters identify will */
1414 : /* depend on the returned method, but they will all have been */
1415 : /* normalized into degrees and meters. */
1416 : /************************************************************************/
1417 :
1418 2960 : int GTIFGetProjTRFInfoEx( void* ctxIn,
1419 : int nProjTRFCode,
1420 : char **ppszProjTRFName,
1421 : short * pnProjMethod,
1422 : double * padfProjParams )
1423 :
1424 : {
1425 2960 : if ((nProjTRFCode >= Proj_UTM_zone_1N && nProjTRFCode <= Proj_UTM_zone_60N) ||
1426 48 : (nProjTRFCode >= Proj_UTM_zone_1S && nProjTRFCode <= Proj_UTM_zone_60S))
1427 : {
1428 : int bNorth;
1429 : int nZone;
1430 2810 : if (nProjTRFCode <= Proj_UTM_zone_60N)
1431 : {
1432 2807 : bNorth = TRUE;
1433 2807 : nZone = nProjTRFCode - Proj_UTM_zone_1N + 1;
1434 : }
1435 : else
1436 : {
1437 3 : bNorth = FALSE;
1438 3 : nZone = nProjTRFCode - Proj_UTM_zone_1S + 1;
1439 : }
1440 :
1441 2810 : if (ppszProjTRFName)
1442 : {
1443 : char szProjTRFName[64];
1444 0 : snprintf(szProjTRFName, sizeof(szProjTRFName), "UTM zone %d%c",
1445 : nZone, (bNorth) ? 'N' : 'S');
1446 0 : *ppszProjTRFName = CPLStrdup(szProjTRFName);
1447 : }
1448 :
1449 2810 : if (pnProjMethod)
1450 2810 : *pnProjMethod = 9807;
1451 :
1452 2810 : if (padfProjParams)
1453 : {
1454 2810 : padfProjParams[0] = 0;
1455 2810 : padfProjParams[1] = -183 + 6 * nZone;
1456 2810 : padfProjParams[2] = 0;
1457 2810 : padfProjParams[3] = 0;
1458 2810 : padfProjParams[4] = 0.9996;
1459 2810 : padfProjParams[5] = 500000;
1460 2810 : padfProjParams[6] = (bNorth) ? 0 : 10000000;
1461 : }
1462 :
1463 2810 : return TRUE;
1464 : }
1465 :
1466 150 : if( nProjTRFCode == KvUserDefined )
1467 0 : return FALSE;
1468 :
1469 : {
1470 : char szCode[12];
1471 150 : snprintf(szCode, sizeof(szCode), "%d", nProjTRFCode);
1472 150 : PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
1473 150 : PJ *transf = proj_create_from_database(
1474 : ctx, "EPSG", szCode, PJ_CATEGORY_COORDINATE_OPERATION, 0, NULL);
1475 150 : if( !transf )
1476 : {
1477 0 : return FALSE;
1478 : }
1479 :
1480 150 : if( proj_get_type(transf) != PJ_TYPE_CONVERSION )
1481 : {
1482 0 : proj_destroy(transf);
1483 0 : return FALSE;
1484 : }
1485 :
1486 : /* Get the projection method code */
1487 150 : const char* pszMethodCode = NULL;
1488 150 : proj_coordoperation_get_method_info(ctx, transf,
1489 : NULL, /* method name */
1490 : NULL, /* method auth name (should be EPSG) */
1491 : &pszMethodCode);
1492 150 : assert( pszMethodCode );
1493 150 : const int nProjMethod = atoi(pszMethodCode);
1494 :
1495 : /* -------------------------------------------------------------------- */
1496 : /* Initialize a definition of what EPSG codes need to be loaded */
1497 : /* into what fields in adfProjParams. */
1498 : /* -------------------------------------------------------------------- */
1499 150 : const int nCTProjMethod = EPSGProjMethodToCTProjMethod( nProjMethod, TRUE );
1500 : int anEPSGCodes[7];
1501 150 : SetGTParamIds( nCTProjMethod, nProjMethod, NULL, anEPSGCodes );
1502 :
1503 :
1504 : /* -------------------------------------------------------------------- */
1505 : /* Get the parameters for this projection. */
1506 : /* -------------------------------------------------------------------- */
1507 :
1508 : double adfProjParams[7];
1509 :
1510 1200 : for( int i = 0; i < 7; i++ )
1511 : {
1512 1050 : int nEPSGCode = anEPSGCodes[i];
1513 :
1514 : /* Establish default */
1515 1050 : if( nEPSGCode == EPSGAngleRectifiedToSkewedGrid )
1516 3 : adfProjParams[i] = 90.0;
1517 1047 : else if( nEPSGCode == EPSGNatOriginScaleFactor
1518 918 : || nEPSGCode == EPSGInitialLineScaleFactor
1519 915 : || nEPSGCode == EPSGPseudoStdParallelScaleFactor )
1520 132 : adfProjParams[i] = 1.0;
1521 : else
1522 915 : adfProjParams[i] = 0.0;
1523 :
1524 : /* If there is no parameter, skip */
1525 1050 : if( nEPSGCode == 0 )
1526 381 : continue;
1527 :
1528 732 : const int nParamCount = proj_coordoperation_get_param_count(ctx, transf);
1529 :
1530 : /* Find the matching parameter */
1531 732 : const char *pszUOMCategory = NULL;
1532 732 : double dfValue = 0.0;
1533 732 : double dfUnitConvFactor = 0.0;
1534 : int iEPSG; /* Used after for */
1535 2241 : for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
1536 : {
1537 2172 : const char* pszParamCode = NULL;
1538 2172 : proj_coordoperation_get_param(
1539 : ctx, transf, iEPSG,
1540 : NULL, /* name */
1541 : NULL, /* auth name */
1542 : &pszParamCode,
1543 : &dfValue,
1544 : NULL, /* value (string) */
1545 : &dfUnitConvFactor, /* unit conv factor */
1546 : NULL, /* unit name */
1547 : NULL, /* unit auth name */
1548 : NULL, /* unit code */
1549 : &pszUOMCategory /* unit category */);
1550 2172 : assert(pszParamCode);
1551 2172 : if( atoi(pszParamCode) == nEPSGCode )
1552 : {
1553 663 : break;
1554 : }
1555 : }
1556 :
1557 : /* not found, accept the default */
1558 732 : if( iEPSG == nParamCount )
1559 : {
1560 : /* for CT_ObliqueMercator try alternate parameter codes first */
1561 : /* because EPSG proj method 9812 uses EPSGFalseXXXXX, but 9815 uses EPSGProjCenterXXXXX */
1562 69 : if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterEasting )
1563 1 : nEPSGCode = EPSGFalseEasting;
1564 68 : else if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterNorthing )
1565 1 : nEPSGCode = EPSGFalseNorthing;
1566 : /* for CT_PolarStereographic try alternate parameter codes first */
1567 : /* because EPSG proj method 9829 uses EPSGLatOfStdParallel instead of EPSGNatOriginLat */
1568 : /* and EPSGOriginLong instead of EPSGNatOriginLong */
1569 67 : else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLat )
1570 2 : nEPSGCode = EPSGLatOfStdParallel;
1571 65 : else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLong )
1572 2 : nEPSGCode = EPSGOriginLong;
1573 : else
1574 63 : continue;
1575 :
1576 19 : for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
1577 : {
1578 19 : const char* pszParamCode = NULL;
1579 19 : proj_coordoperation_get_param(
1580 : ctx, transf, iEPSG,
1581 : NULL, /* name */
1582 : NULL, /* auth name */
1583 : &pszParamCode,
1584 : &dfValue,
1585 : NULL, /* value (string) */
1586 : &dfUnitConvFactor, /* unit conv factor */
1587 : NULL, /* unit name */
1588 : NULL, /* unit auth name */
1589 : NULL, /* unit code */
1590 : &pszUOMCategory /* unit category */);
1591 19 : assert(pszParamCode);
1592 19 : if( atoi(pszParamCode) == nEPSGCode )
1593 : {
1594 6 : break;
1595 : }
1596 : }
1597 :
1598 6 : if( iEPSG == nParamCount )
1599 0 : continue;
1600 : }
1601 :
1602 669 : assert(pszUOMCategory);
1603 :
1604 669 : adfProjParams[i] = dfValue * dfUnitConvFactor;
1605 669 : if( strcmp(pszUOMCategory, "angular") == 0.0 )
1606 : {
1607 : /* Convert from radians to degrees */
1608 306 : adfProjParams[i] *= 180 / M_PI;
1609 : }
1610 : }
1611 :
1612 : /* -------------------------------------------------------------------- */
1613 : /* Get the name, if requested. */
1614 : /* -------------------------------------------------------------------- */
1615 150 : if( ppszProjTRFName != NULL )
1616 : {
1617 0 : const char* pszName = proj_get_name(transf);
1618 0 : if( !pszName )
1619 : {
1620 : // shouldn't happen
1621 0 : proj_destroy(transf);
1622 0 : return FALSE;
1623 : }
1624 0 : *ppszProjTRFName = CPLStrdup(pszName);
1625 : }
1626 :
1627 : /* -------------------------------------------------------------------- */
1628 : /* Transfer requested data into passed variables. */
1629 : /* -------------------------------------------------------------------- */
1630 150 : if( pnProjMethod != NULL )
1631 150 : *pnProjMethod = (short) nProjMethod;
1632 :
1633 150 : if( padfProjParams != NULL )
1634 : {
1635 1200 : for( int i = 0; i < 7; i++ )
1636 1050 : padfProjParams[i] = adfProjParams[i];
1637 : }
1638 :
1639 150 : proj_destroy(transf);
1640 :
1641 150 : return TRUE;
1642 : }
1643 : }
1644 :
1645 0 : int GTIFGetProjTRFInfo( /* Conversion code */
1646 : int nProjTRFCode,
1647 : char **ppszProjTRFName,
1648 : short * pnProjMethod,
1649 : double * padfProjParams )
1650 : {
1651 0 : PJ_CONTEXT* ctx = proj_context_create();
1652 0 : const int ret = GTIFGetProjTRFInfoEx(
1653 : ctx, nProjTRFCode, ppszProjTRFName, pnProjMethod, padfProjParams);
1654 0 : proj_context_destroy(ctx);
1655 0 : return ret;
1656 : }
1657 :
1658 : /************************************************************************/
1659 : /* GTIFFetchProjParms() */
1660 : /* */
1661 : /* Fetch the projection parameters for a particular projection */
1662 : /* from a GeoTIFF file, and fill the GTIFDefn structure out */
1663 : /* with them. */
1664 : /************************************************************************/
1665 :
1666 74 : static void GTIFFetchProjParms( GTIF * psGTIF, GTIFDefn * psDefn )
1667 : {
1668 :
1669 : /* -------------------------------------------------------------------- */
1670 : /* Get the false easting, and northing if available. */
1671 : /* -------------------------------------------------------------------- */
1672 74 : double dfFalseEasting = 0.0;
1673 74 : if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseEastingGeoKey, &dfFalseEasting, 0, 1)
1674 12 : && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterEastingGeoKey,
1675 : &dfFalseEasting, 0, 1)
1676 12 : && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginEastingGeoKey,
1677 : &dfFalseEasting, 0, 1) )
1678 0 : dfFalseEasting = 0.0;
1679 :
1680 74 : double dfFalseNorthing = 0.0;
1681 74 : if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseNorthingGeoKey, &dfFalseNorthing,0,1)
1682 12 : && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterNorthingGeoKey,
1683 : &dfFalseNorthing, 0, 1)
1684 12 : && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginNorthingGeoKey,
1685 : &dfFalseNorthing, 0, 1) )
1686 0 : dfFalseNorthing = 0.0;
1687 :
1688 74 : double dfNatOriginLong = 0.0, dfNatOriginLat = 0.0, dfRectGridAngle = 0.0;
1689 74 : double dfNatOriginScale = 1.0;
1690 74 : double dfStdParallel1 = 0.0, dfStdParallel2 = 0.0, dfAzimuth = 0.0;
1691 :
1692 74 : switch( psDefn->CTProjection )
1693 : {
1694 : /* -------------------------------------------------------------------- */
1695 0 : case CT_Stereographic:
1696 : /* -------------------------------------------------------------------- */
1697 0 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1698 : &dfNatOriginLong, 0, 1 ) == 0
1699 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1700 : &dfNatOriginLong, 0, 1 ) == 0
1701 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1702 : &dfNatOriginLong, 0, 1 ) == 0 )
1703 0 : dfNatOriginLong = 0.0;
1704 :
1705 0 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1706 : &dfNatOriginLat, 0, 1 ) == 0
1707 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1708 : &dfNatOriginLat, 0, 1 ) == 0
1709 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1710 : &dfNatOriginLat, 0, 1 ) == 0 )
1711 0 : dfNatOriginLat = 0.0;
1712 :
1713 0 : if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1714 : &dfNatOriginScale, 0, 1 ) == 0 )
1715 0 : dfNatOriginScale = 1.0;
1716 :
1717 : /* notdef: should transform to decimal degrees at this point */
1718 :
1719 0 : psDefn->ProjParm[0] = dfNatOriginLat;
1720 0 : psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1721 0 : psDefn->ProjParm[1] = dfNatOriginLong;
1722 0 : psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1723 0 : psDefn->ProjParm[4] = dfNatOriginScale;
1724 0 : psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1725 0 : psDefn->ProjParm[5] = dfFalseEasting;
1726 0 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1727 0 : psDefn->ProjParm[6] = dfFalseNorthing;
1728 0 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1729 :
1730 0 : psDefn->nParms = 7;
1731 0 : break;
1732 :
1733 : /* -------------------------------------------------------------------- */
1734 8 : case CT_Mercator:
1735 : /* -------------------------------------------------------------------- */
1736 8 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1737 : &dfNatOriginLong, 0, 1 ) == 0
1738 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1739 : &dfNatOriginLong, 0, 1 ) == 0
1740 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1741 : &dfNatOriginLong, 0, 1 ) == 0 )
1742 0 : dfNatOriginLong = 0.0;
1743 :
1744 8 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1745 : &dfNatOriginLat, 0, 1 ) == 0
1746 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1747 : &dfNatOriginLat, 0, 1 ) == 0
1748 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1749 : &dfNatOriginLat, 0, 1 ) == 0 )
1750 0 : dfNatOriginLat = 0.0;
1751 :
1752 :
1753 8 : const int bHaveSP1 = GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
1754 : &dfStdParallel1, 0, 1 );
1755 :
1756 8 : int bHaveNOS = GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1757 : &dfNatOriginScale, 0, 1 );
1758 :
1759 : /* Default scale only if dfStdParallel1 isn't defined either */
1760 8 : if( !bHaveNOS && !bHaveSP1)
1761 : {
1762 0 : bHaveNOS = TRUE;
1763 0 : dfNatOriginScale = 1.0;
1764 : }
1765 :
1766 : /* notdef: should transform to decimal degrees at this point */
1767 :
1768 8 : psDefn->ProjParm[0] = dfNatOriginLat;
1769 8 : psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1770 8 : psDefn->ProjParm[1] = dfNatOriginLong;
1771 8 : psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1772 8 : if( bHaveSP1 )
1773 : {
1774 2 : psDefn->ProjParm[2] = dfStdParallel1;
1775 2 : psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
1776 : }
1777 8 : if( bHaveNOS )
1778 : {
1779 7 : psDefn->ProjParm[4] = dfNatOriginScale;
1780 7 : psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1781 : }
1782 8 : psDefn->ProjParm[5] = dfFalseEasting;
1783 8 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1784 8 : psDefn->ProjParm[6] = dfFalseNorthing;
1785 8 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1786 :
1787 8 : psDefn->nParms = 7;
1788 8 : break;
1789 :
1790 : /* -------------------------------------------------------------------- */
1791 6 : case CT_LambertConfConic_1SP:
1792 : case CT_ObliqueStereographic:
1793 : case CT_TransverseMercator:
1794 : case CT_TransvMercator_SouthOriented:
1795 : /* -------------------------------------------------------------------- */
1796 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1797 : &dfNatOriginLong, 0, 1 ) == 0
1798 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1799 : &dfNatOriginLong, 0, 1 ) == 0
1800 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1801 : &dfNatOriginLong, 0, 1 ) == 0 )
1802 0 : dfNatOriginLong = 0.0;
1803 :
1804 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1805 : &dfNatOriginLat, 0, 1 ) == 0
1806 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1807 : &dfNatOriginLat, 0, 1 ) == 0
1808 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1809 : &dfNatOriginLat, 0, 1 ) == 0 )
1810 0 : dfNatOriginLat = 0.0;
1811 :
1812 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1813 : &dfNatOriginScale, 0, 1 ) == 0
1814 : /* See https://github.com/OSGeo/gdal/files/1665718/lasinfo.txt */
1815 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1816 : &dfNatOriginScale, 0, 1 ) == 0 )
1817 0 : dfNatOriginScale = 1.0;
1818 :
1819 : /* notdef: should transform to decimal degrees at this point */
1820 :
1821 6 : psDefn->ProjParm[0] = dfNatOriginLat;
1822 6 : psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1823 6 : psDefn->ProjParm[1] = dfNatOriginLong;
1824 6 : psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1825 6 : psDefn->ProjParm[4] = dfNatOriginScale;
1826 6 : psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1827 6 : psDefn->ProjParm[5] = dfFalseEasting;
1828 6 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1829 6 : psDefn->ProjParm[6] = dfFalseNorthing;
1830 6 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1831 :
1832 6 : psDefn->nParms = 7;
1833 6 : break;
1834 :
1835 : /* -------------------------------------------------------------------- */
1836 3 : case CT_ObliqueMercator: /* hotine */
1837 : case CT_HotineObliqueMercatorAzimuthCenter:
1838 : /* -------------------------------------------------------------------- */
1839 3 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1840 : &dfNatOriginLong, 0, 1 ) == 0
1841 3 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1842 : &dfNatOriginLong, 0, 1 ) == 0
1843 3 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1844 : &dfNatOriginLong, 0, 1 ) == 0 )
1845 0 : dfNatOriginLong = 0.0;
1846 :
1847 3 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1848 : &dfNatOriginLat, 0, 1 ) == 0
1849 3 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1850 : &dfNatOriginLat, 0, 1 ) == 0
1851 3 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1852 : &dfNatOriginLat, 0, 1 ) == 0 )
1853 0 : dfNatOriginLat = 0.0;
1854 :
1855 3 : if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
1856 : &dfAzimuth, 0, 1 ) == 0 )
1857 0 : dfAzimuth = 0.0;
1858 :
1859 3 : if( GTIFKeyGetDOUBLE(psGTIF, ProjRectifiedGridAngleGeoKey,
1860 : &dfRectGridAngle, 0, 1 ) == 0 )
1861 0 : dfRectGridAngle = 90.0;
1862 :
1863 3 : if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1864 : &dfNatOriginScale, 0, 1 ) == 0
1865 3 : && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1866 : &dfNatOriginScale, 0, 1 ) == 0 )
1867 0 : dfNatOriginScale = 1.0;
1868 :
1869 : /* notdef: should transform to decimal degrees at this point */
1870 :
1871 3 : psDefn->ProjParm[0] = dfNatOriginLat;
1872 3 : psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1873 3 : psDefn->ProjParm[1] = dfNatOriginLong;
1874 3 : psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1875 3 : psDefn->ProjParm[2] = dfAzimuth;
1876 3 : psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
1877 3 : psDefn->ProjParm[3] = dfRectGridAngle;
1878 3 : psDefn->ProjParmId[3] = ProjRectifiedGridAngleGeoKey;
1879 3 : psDefn->ProjParm[4] = dfNatOriginScale;
1880 3 : psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
1881 3 : psDefn->ProjParm[5] = dfFalseEasting;
1882 3 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1883 3 : psDefn->ProjParm[6] = dfFalseNorthing;
1884 3 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1885 :
1886 3 : psDefn->nParms = 7;
1887 3 : break;
1888 :
1889 : /* -------------------------------------------------------------------- */
1890 0 : case CT_ObliqueMercator_Laborde:
1891 : /* -------------------------------------------------------------------- */
1892 0 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1893 : &dfNatOriginLong, 0, 1 ) == 0
1894 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1895 : &dfNatOriginLong, 0, 1 ) == 0
1896 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1897 : &dfNatOriginLong, 0, 1 ) == 0 )
1898 0 : dfNatOriginLong = 0.0;
1899 :
1900 0 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1901 : &dfNatOriginLat, 0, 1 ) == 0
1902 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1903 : &dfNatOriginLat, 0, 1 ) == 0
1904 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1905 : &dfNatOriginLat, 0, 1 ) == 0 )
1906 0 : dfNatOriginLat = 0.0;
1907 :
1908 0 : if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
1909 : &dfAzimuth, 0, 1 ) == 0 )
1910 0 : dfAzimuth = 0.0;
1911 :
1912 0 : if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1913 : &dfNatOriginScale, 0, 1 ) == 0
1914 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1915 : &dfNatOriginScale, 0, 1 ) == 0 )
1916 0 : dfNatOriginScale = 1.0;
1917 :
1918 : /* notdef: should transform to decimal degrees at this point */
1919 :
1920 0 : psDefn->ProjParm[0] = dfNatOriginLat;
1921 0 : psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1922 0 : psDefn->ProjParm[1] = dfNatOriginLong;
1923 0 : psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1924 0 : psDefn->ProjParm[2] = dfAzimuth;
1925 0 : psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
1926 0 : psDefn->ProjParm[4] = dfNatOriginScale;
1927 0 : psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
1928 0 : psDefn->ProjParm[5] = dfFalseEasting;
1929 0 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1930 0 : psDefn->ProjParm[6] = dfFalseNorthing;
1931 0 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1932 :
1933 0 : psDefn->nParms = 7;
1934 0 : break;
1935 :
1936 : /* -------------------------------------------------------------------- */
1937 2 : case CT_CassiniSoldner:
1938 : case CT_Polyconic:
1939 : /* -------------------------------------------------------------------- */
1940 2 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1941 : &dfNatOriginLong, 0, 1 ) == 0
1942 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1943 : &dfNatOriginLong, 0, 1 ) == 0
1944 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1945 : &dfNatOriginLong, 0, 1 ) == 0 )
1946 0 : dfNatOriginLong = 0.0;
1947 :
1948 2 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1949 : &dfNatOriginLat, 0, 1 ) == 0
1950 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1951 : &dfNatOriginLat, 0, 1 ) == 0
1952 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1953 : &dfNatOriginLat, 0, 1 ) == 0 )
1954 0 : dfNatOriginLat = 0.0;
1955 :
1956 2 : if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1957 : &dfNatOriginScale, 0, 1 ) == 0
1958 1 : && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1959 : &dfNatOriginScale, 0, 1 ) == 0 )
1960 1 : dfNatOriginScale = 1.0;
1961 :
1962 : /* notdef: should transform to decimal degrees at this point */
1963 :
1964 2 : psDefn->ProjParm[0] = dfNatOriginLat;
1965 2 : psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1966 2 : psDefn->ProjParm[1] = dfNatOriginLong;
1967 2 : psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1968 2 : psDefn->ProjParm[4] = dfNatOriginScale;
1969 2 : psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1970 2 : psDefn->ProjParm[5] = dfFalseEasting;
1971 2 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1972 2 : psDefn->ProjParm[6] = dfFalseNorthing;
1973 2 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1974 :
1975 2 : psDefn->nParms = 7;
1976 2 : break;
1977 :
1978 : /* -------------------------------------------------------------------- */
1979 18 : case CT_AzimuthalEquidistant:
1980 : case CT_MillerCylindrical:
1981 : case CT_Gnomonic:
1982 : case CT_LambertAzimEqualArea:
1983 : case CT_Orthographic:
1984 : case CT_NewZealandMapGrid:
1985 : /* -------------------------------------------------------------------- */
1986 18 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1987 : &dfNatOriginLong, 0, 1 ) == 0
1988 17 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1989 : &dfNatOriginLong, 0, 1 ) == 0
1990 17 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1991 : &dfNatOriginLong, 0, 1 ) == 0 )
1992 0 : dfNatOriginLong = 0.0;
1993 :
1994 18 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1995 : &dfNatOriginLat, 0, 1 ) == 0
1996 17 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1997 : &dfNatOriginLat, 0, 1 ) == 0
1998 17 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1999 : &dfNatOriginLat, 0, 1 ) == 0 )
2000 0 : dfNatOriginLat = 0.0;
2001 :
2002 : /* notdef: should transform to decimal degrees at this point */
2003 :
2004 18 : psDefn->ProjParm[0] = dfNatOriginLat;
2005 18 : psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
2006 18 : psDefn->ProjParm[1] = dfNatOriginLong;
2007 18 : psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
2008 18 : psDefn->ProjParm[5] = dfFalseEasting;
2009 18 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2010 18 : psDefn->ProjParm[6] = dfFalseNorthing;
2011 18 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2012 :
2013 18 : psDefn->nParms = 7;
2014 18 : break;
2015 :
2016 : /* -------------------------------------------------------------------- */
2017 5 : case CT_Equirectangular:
2018 : /* -------------------------------------------------------------------- */
2019 5 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2020 : &dfNatOriginLong, 0, 1 ) == 0
2021 5 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2022 : &dfNatOriginLong, 0, 1 ) == 0
2023 5 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2024 : &dfNatOriginLong, 0, 1 ) == 0 )
2025 0 : dfNatOriginLong = 0.0;
2026 :
2027 5 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2028 : &dfNatOriginLat, 0, 1 ) == 0
2029 5 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2030 : &dfNatOriginLat, 0, 1 ) == 0
2031 5 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2032 : &dfNatOriginLat, 0, 1 ) == 0 )
2033 0 : dfNatOriginLat = 0.0;
2034 :
2035 5 : if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2036 : &dfStdParallel1, 0, 1 ) == 0 )
2037 0 : dfStdParallel1 = 0.0;
2038 :
2039 : /* notdef: should transform to decimal degrees at this point */
2040 :
2041 5 : psDefn->ProjParm[0] = dfNatOriginLat;
2042 5 : psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
2043 5 : psDefn->ProjParm[1] = dfNatOriginLong;
2044 5 : psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
2045 5 : psDefn->ProjParm[2] = dfStdParallel1;
2046 5 : psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
2047 5 : psDefn->ProjParm[5] = dfFalseEasting;
2048 5 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2049 5 : psDefn->ProjParm[6] = dfFalseNorthing;
2050 5 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2051 :
2052 5 : psDefn->nParms = 7;
2053 5 : break;
2054 :
2055 : /* -------------------------------------------------------------------- */
2056 7 : case CT_Robinson:
2057 : case CT_Sinusoidal:
2058 : case CT_VanDerGrinten:
2059 : /* -------------------------------------------------------------------- */
2060 7 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2061 : &dfNatOriginLong, 0, 1 ) == 0
2062 7 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2063 : &dfNatOriginLong, 0, 1 ) == 0
2064 7 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2065 : &dfNatOriginLong, 0, 1 ) == 0 )
2066 0 : dfNatOriginLong = 0.0;
2067 :
2068 : /* notdef: should transform to decimal degrees at this point */
2069 :
2070 7 : psDefn->ProjParm[1] = dfNatOriginLong;
2071 7 : psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
2072 7 : psDefn->ProjParm[5] = dfFalseEasting;
2073 7 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2074 7 : psDefn->ProjParm[6] = dfFalseNorthing;
2075 7 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2076 :
2077 7 : psDefn->nParms = 7;
2078 7 : break;
2079 :
2080 : /* -------------------------------------------------------------------- */
2081 6 : case CT_PolarStereographic:
2082 : /* -------------------------------------------------------------------- */
2083 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjStraightVertPoleLongGeoKey,
2084 : &dfNatOriginLong, 0, 1 ) == 0
2085 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2086 : &dfNatOriginLong, 0, 1 ) == 0
2087 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2088 : &dfNatOriginLong, 0, 1 ) == 0
2089 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2090 : &dfNatOriginLong, 0, 1 ) == 0 )
2091 0 : dfNatOriginLong = 0.0;
2092 :
2093 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2094 : &dfNatOriginLat, 0, 1 ) == 0
2095 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2096 : &dfNatOriginLat, 0, 1 ) == 0
2097 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2098 : &dfNatOriginLat, 0, 1 ) == 0 )
2099 0 : dfNatOriginLat = 0.0;
2100 :
2101 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
2102 : &dfNatOriginScale, 0, 1 ) == 0
2103 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
2104 : &dfNatOriginScale, 0, 1 ) == 0 )
2105 0 : dfNatOriginScale = 1.0;
2106 :
2107 : /* notdef: should transform to decimal degrees at this point */
2108 :
2109 6 : psDefn->ProjParm[0] = dfNatOriginLat;
2110 6 : psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;;
2111 6 : psDefn->ProjParm[1] = dfNatOriginLong;
2112 6 : psDefn->ProjParmId[1] = ProjStraightVertPoleLongGeoKey;
2113 6 : psDefn->ProjParm[4] = dfNatOriginScale;
2114 6 : psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
2115 6 : psDefn->ProjParm[5] = dfFalseEasting;
2116 6 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2117 6 : psDefn->ProjParm[6] = dfFalseNorthing;
2118 6 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2119 :
2120 6 : psDefn->nParms = 7;
2121 6 : break;
2122 :
2123 : /* -------------------------------------------------------------------- */
2124 12 : case CT_LambertConfConic_2SP:
2125 : /* -------------------------------------------------------------------- */
2126 12 : if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2127 : &dfStdParallel1, 0, 1 ) == 0 )
2128 0 : dfStdParallel1 = 0.0;
2129 :
2130 12 : if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
2131 : &dfStdParallel2, 0, 1 ) == 0 )
2132 0 : dfStdParallel1 = 0.0;
2133 :
2134 12 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2135 : &dfNatOriginLong, 0, 1 ) == 0
2136 12 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2137 : &dfNatOriginLong, 0, 1 ) == 0
2138 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2139 : &dfNatOriginLong, 0, 1 ) == 0 )
2140 0 : dfNatOriginLong = 0.0;
2141 :
2142 12 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2143 : &dfNatOriginLat, 0, 1 ) == 0
2144 12 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2145 : &dfNatOriginLat, 0, 1 ) == 0
2146 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2147 : &dfNatOriginLat, 0, 1 ) == 0 )
2148 0 : dfNatOriginLat = 0.0;
2149 :
2150 : /* notdef: should transform to decimal degrees at this point */
2151 :
2152 12 : psDefn->ProjParm[0] = dfNatOriginLat;
2153 12 : psDefn->ProjParmId[0] = ProjFalseOriginLatGeoKey;
2154 12 : psDefn->ProjParm[1] = dfNatOriginLong;
2155 12 : psDefn->ProjParmId[1] = ProjFalseOriginLongGeoKey;
2156 12 : psDefn->ProjParm[2] = dfStdParallel1;
2157 12 : psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
2158 12 : psDefn->ProjParm[3] = dfStdParallel2;
2159 12 : psDefn->ProjParmId[3] = ProjStdParallel2GeoKey;
2160 12 : psDefn->ProjParm[5] = dfFalseEasting;
2161 12 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2162 12 : psDefn->ProjParm[6] = dfFalseNorthing;
2163 12 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2164 :
2165 12 : psDefn->nParms = 7;
2166 12 : break;
2167 :
2168 : /* -------------------------------------------------------------------- */
2169 6 : case CT_AlbersEqualArea:
2170 : case CT_EquidistantConic:
2171 : /* -------------------------------------------------------------------- */
2172 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2173 : &dfStdParallel1, 0, 1 ) == 0 )
2174 0 : dfStdParallel1 = 0.0;
2175 :
2176 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
2177 : &dfStdParallel2, 0, 1 ) == 0 )
2178 0 : dfStdParallel2 = 0.0;
2179 :
2180 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2181 : &dfNatOriginLong, 0, 1 ) == 0
2182 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2183 : &dfNatOriginLong, 0, 1 ) == 0
2184 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2185 : &dfNatOriginLong, 0, 1 ) == 0 )
2186 0 : dfNatOriginLong = 0.0;
2187 :
2188 6 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2189 : &dfNatOriginLat, 0, 1 ) == 0
2190 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2191 : &dfNatOriginLat, 0, 1 ) == 0
2192 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2193 : &dfNatOriginLat, 0, 1 ) == 0 )
2194 0 : dfNatOriginLat = 0.0;
2195 :
2196 : /* notdef: should transform to decimal degrees at this point */
2197 :
2198 6 : psDefn->ProjParm[0] = dfStdParallel1;
2199 6 : psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
2200 6 : psDefn->ProjParm[1] = dfStdParallel2;
2201 6 : psDefn->ProjParmId[1] = ProjStdParallel2GeoKey;
2202 6 : psDefn->ProjParm[2] = dfNatOriginLat;
2203 6 : psDefn->ProjParmId[2] = ProjNatOriginLatGeoKey;
2204 6 : psDefn->ProjParm[3] = dfNatOriginLong;
2205 6 : psDefn->ProjParmId[3] = ProjNatOriginLongGeoKey;
2206 6 : psDefn->ProjParm[5] = dfFalseEasting;
2207 6 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2208 6 : psDefn->ProjParm[6] = dfFalseNorthing;
2209 6 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2210 :
2211 6 : psDefn->nParms = 7;
2212 6 : break;
2213 :
2214 : /* -------------------------------------------------------------------- */
2215 1 : case CT_CylindricalEqualArea:
2216 : /* -------------------------------------------------------------------- */
2217 1 : if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2218 : &dfStdParallel1, 0, 1 ) == 0 )
2219 0 : dfStdParallel1 = 0.0;
2220 :
2221 1 : if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2222 : &dfNatOriginLong, 0, 1 ) == 0
2223 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2224 : &dfNatOriginLong, 0, 1 ) == 0
2225 0 : && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2226 : &dfNatOriginLong, 0, 1 ) == 0 )
2227 0 : dfNatOriginLong = 0.0;
2228 :
2229 : /* notdef: should transform to decimal degrees at this point */
2230 :
2231 1 : psDefn->ProjParm[0] = dfStdParallel1;
2232 1 : psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
2233 1 : psDefn->ProjParm[1] = dfNatOriginLong;
2234 1 : psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
2235 1 : psDefn->ProjParm[5] = dfFalseEasting;
2236 1 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2237 1 : psDefn->ProjParm[6] = dfFalseNorthing;
2238 1 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2239 :
2240 1 : psDefn->nParms = 7;
2241 1 : break;
2242 : }
2243 :
2244 : /* -------------------------------------------------------------------- */
2245 : /* Normalize any linear parameters into meters. In GeoTIFF */
2246 : /* the linear projection parameter tags are normally in the */
2247 : /* units of the coordinate system described. */
2248 : /* -------------------------------------------------------------------- */
2249 592 : for( int iParam = 0; iParam < psDefn->nParms; iParam++ )
2250 : {
2251 518 : switch( psDefn->ProjParmId[iParam] )
2252 : {
2253 148 : case ProjFalseEastingGeoKey:
2254 : case ProjFalseNorthingGeoKey:
2255 : case ProjFalseOriginEastingGeoKey:
2256 : case ProjFalseOriginNorthingGeoKey:
2257 : case ProjCenterEastingGeoKey:
2258 : case ProjCenterNorthingGeoKey:
2259 148 : if( psDefn->UOMLengthInMeters != 0
2260 148 : && psDefn->UOMLengthInMeters != 1.0 )
2261 : {
2262 26 : psDefn->ProjParm[iParam] *= psDefn->UOMLengthInMeters;
2263 : }
2264 148 : break;
2265 :
2266 370 : default:
2267 370 : break;
2268 : }
2269 : }
2270 74 : }
2271 :
2272 : /************************************************************************/
2273 : /* GTIFGetDefn() */
2274 : /************************************************************************/
2275 :
2276 : /**
2277 : @param psGTIF GeoTIFF information handle as returned by GTIFNew.
2278 : @param psDefn Pointer to an existing GTIFDefn structure allocated by GTIFAllocDefn().
2279 :
2280 : @return TRUE if the function has been successful, otherwise FALSE.
2281 :
2282 : This function reads the coordinate system definition from a GeoTIFF file,
2283 : and <i>normalizes</i> it into a set of component information using
2284 : definitions from the EPSG database as provided by the PROJ library.
2285 : This function is intended to simplify correct support for
2286 : reading files with defined PCS (Projected Coordinate System) codes that
2287 : wouldn't otherwise be directly known by application software by reducing
2288 : it to the underlying projection method, parameters, datum, ellipsoid,
2289 : prime meridian and units.<p>
2290 :
2291 : The application should pass a pointer to an existing uninitialized
2292 : GTIFDefn structure, and GTIFGetDefn() will fill it in. The function
2293 : currently always returns TRUE but in the future will return FALSE if
2294 : the database is not found. In any event, all geokeys actually found in the
2295 : file will be copied into the GTIFDefn. However, if the database isn't
2296 : found, codes implied by other codes will not be set properly.<p>
2297 :
2298 : GTIFGetDefn() will not generally work if the EPSG derived database cannot
2299 : be found.<p>
2300 :
2301 : The normalization methodology operates by fetching tags from the GeoTIFF
2302 : file, and then setting all other tags implied by them in the structure. The
2303 : implied relationships are worked out by reading definitions from the
2304 : various EPSG derived database tables.<p>
2305 :
2306 : For instance, if a PCS (ProjectedCSTypeGeoKey) is found in the GeoTIFF file
2307 : this code is used to lookup a record in the database.
2308 : For example given the PCS 26746 we can find the name
2309 : (NAD27 / California zone VI), the GCS 4257 (NAD27), and the ProjectionCode
2310 : 10406 (California CS27 zone VI). The GCS, and ProjectionCode can in turn
2311 : be looked up in other tables until all the details of units, ellipsoid,
2312 : prime meridian, datum, projection (LambertConfConic_2SP) and projection
2313 : parameters are established. A full listgeo dump of a file
2314 : for this result might look like the following, all based on a single PCS
2315 : value:<p>
2316 :
2317 : <pre>
2318 : % listgeo -norm ~/data/geotiff/pci_eg/spaf27.tif
2319 : Geotiff_Information:
2320 : Version: 1
2321 : Key_Revision: 1.0
2322 : Tagged_Information:
2323 : ModelTiepointTag (2,3):
2324 : 0 0 0
2325 : 1577139.71 634349.176 0
2326 : ModelPixelScaleTag (1,3):
2327 : 195.509321 198.32184 0
2328 : End_Of_Tags.
2329 : Keyed_Information:
2330 : GTModelTypeGeoKey (Short,1): ModelTypeProjected
2331 : GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
2332 : ProjectedCSTypeGeoKey (Short,1): PCS_NAD27_California_VI
2333 : End_Of_Keys.
2334 : End_Of_Geotiff.
2335 :
2336 : PCS = 26746 (NAD27 / California zone VI)
2337 : Projection = 10406 (California CS27 zone VI)
2338 : Projection Method: CT_LambertConfConic_2SP
2339 : ProjStdParallel1GeoKey: 33.883333
2340 : ProjStdParallel2GeoKey: 32.766667
2341 : ProjFalseOriginLatGeoKey: 32.166667
2342 : ProjFalseOriginLongGeoKey: -116.233333
2343 : ProjFalseEastingGeoKey: 609601.219202
2344 : ProjFalseNorthingGeoKey: 0.000000
2345 : GCS: 4267/NAD27
2346 : Datum: 6267/North American Datum 1927
2347 : Ellipsoid: 7008/Clarke 1866 (6378206.40,6356583.80)
2348 : Prime Meridian: 8901/Greenwich (0.000000)
2349 : Projection Linear Units: 9003/US survey foot (0.304801m)
2350 : </pre>
2351 :
2352 : Note that GTIFGetDefn() does not inspect or return the tiepoints and scale.
2353 : This must be handled separately as it normally would. It is intended to
2354 : simplify capture and normalization of the coordinate system definition.
2355 : Note that GTIFGetDefn() also does the following things:
2356 :
2357 : <ol>
2358 : <li> Convert all angular values to decimal degrees.
2359 : <li> Convert all linear values to meters.
2360 : <li> Return the linear units and conversion to meters for the tiepoints and
2361 : scale (though the tiepoints and scale remain in their native units).
2362 : <li> When reading projection parameters a variety of differences between
2363 : different GeoTIFF generators are handled, and a normalized set of parameters
2364 : for each projection are always returned.
2365 : </ol>
2366 :
2367 : Code fields in the GTIFDefn are filled with KvUserDefined if there is not
2368 : value to assign. The parameter lists for each of the underlying projection
2369 : transform methods can be found at the
2370 : <a href="http://www.remotesensing.org/geotiff/proj_list">Projections</a>
2371 : page. Note that nParms will be set based on the maximum parameter used.
2372 : Some of the parameters may not be used in which case the
2373 : GTIFDefn::ProjParmId[] will
2374 : be zero. This is done to retain correspondence to the EPSG parameter
2375 : numbering scheme.<p>
2376 :
2377 : The
2378 : <a href="http://www.remotesensing.org/cgi-bin/cvsweb.cgi/~checkout~/osrs/geotiff/libgeotiff/geotiff_proj4.c">geotiff_proj4.c</a> module distributed with libgeotiff can
2379 : be used as an example of code that converts a GTIFDefn into another projection
2380 : system.<p>
2381 :
2382 : @see GTIFKeySet()
2383 :
2384 : */
2385 :
2386 7785 : int GTIFGetDefn( GTIF * psGTIF, GTIFDefn * psDefn )
2387 :
2388 : {
2389 7785 : if( !GTIFGetPROJContext(psGTIF, TRUE, NULL) )
2390 : {
2391 0 : return FALSE;
2392 : }
2393 :
2394 : /* -------------------------------------------------------------------- */
2395 : /* Initially we default all the information we can. */
2396 : /* -------------------------------------------------------------------- */
2397 7785 : psDefn->DefnSet = 1;
2398 7785 : psDefn->Model = KvUserDefined;
2399 7785 : psDefn->PCS = KvUserDefined;
2400 7785 : psDefn->GCS = KvUserDefined;
2401 7785 : psDefn->UOMLength = KvUserDefined;
2402 7785 : psDefn->UOMLengthInMeters = 1.0;
2403 7785 : psDefn->UOMAngle = KvUserDefined;
2404 7785 : psDefn->UOMAngleInDegrees = 1.0;
2405 7785 : psDefn->Datum = KvUserDefined;
2406 7785 : psDefn->Ellipsoid = KvUserDefined;
2407 7785 : psDefn->SemiMajor = 0.0;
2408 7785 : psDefn->SemiMinor = 0.0;
2409 7785 : psDefn->PM = KvUserDefined;
2410 7785 : psDefn->PMLongToGreenwich = 0.0;
2411 : #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2412 7785 : psDefn->TOWGS84Count = 0;
2413 7785 : memset( psDefn->TOWGS84, 0, sizeof(psDefn->TOWGS84) );
2414 : #endif
2415 :
2416 7785 : psDefn->ProjCode = KvUserDefined;
2417 7785 : psDefn->Projection = KvUserDefined;
2418 7785 : psDefn->CTProjection = KvUserDefined;
2419 :
2420 7785 : psDefn->nParms = 0;
2421 85635 : for( int i = 0; i < MAX_GTIF_PROJPARMS; i++ )
2422 : {
2423 77850 : psDefn->ProjParm[i] = 0.0;
2424 77850 : psDefn->ProjParmId[i] = 0;
2425 : }
2426 :
2427 7785 : psDefn->MapSys = KvUserDefined;
2428 7785 : psDefn->Zone = 0;
2429 :
2430 : /* -------------------------------------------------------------------- */
2431 : /* Do we have any geokeys? */
2432 : /* -------------------------------------------------------------------- */
2433 : {
2434 7785 : int nKeyCount = 0;
2435 : int anVersion[3];
2436 7785 : GTIFDirectoryInfo( psGTIF, anVersion, &nKeyCount );
2437 :
2438 7785 : if( nKeyCount == 0 )
2439 : {
2440 2281 : psDefn->DefnSet = 0;
2441 2281 : return FALSE;
2442 : }
2443 : }
2444 :
2445 : /* -------------------------------------------------------------------- */
2446 : /* Try to get the overall model type. */
2447 : /* -------------------------------------------------------------------- */
2448 5504 : GTIFKeyGetSSHORT(psGTIF,GTModelTypeGeoKey,&(psDefn->Model));
2449 :
2450 : /* -------------------------------------------------------------------- */
2451 : /* Extract the Geog units. */
2452 : /* -------------------------------------------------------------------- */
2453 5504 : short nGeogUOMLinear = 9001; /* Linear_Meter */
2454 5504 : if( GTIFKeyGetSSHORT(psGTIF, GeogLinearUnitsGeoKey, &nGeogUOMLinear) == 1 )
2455 : {
2456 23 : psDefn->UOMLength = nGeogUOMLinear;
2457 : }
2458 :
2459 : /* -------------------------------------------------------------------- */
2460 : /* Try to get a PCS. */
2461 : /* -------------------------------------------------------------------- */
2462 5504 : if( GTIFKeyGetSSHORT(psGTIF,ProjectedCSTypeGeoKey, &(psDefn->PCS)) == 1
2463 3035 : && psDefn->PCS != KvUserDefined )
2464 : {
2465 : /*
2466 : * Translate this into useful information.
2467 : */
2468 2947 : GTIFGetPCSInfoEx( psGTIF->pj_context,
2469 2947 : psDefn->PCS, NULL, &(psDefn->ProjCode),
2470 : &(psDefn->UOMLength), &(psDefn->GCS) );
2471 : }
2472 :
2473 : /* -------------------------------------------------------------------- */
2474 : /* If we have the PCS code, but didn't find it in the database */
2475 : /* (likely because we can't find them) we will try some ``jiffy */
2476 : /* rules'' for UTM and state plane. */
2477 : /* -------------------------------------------------------------------- */
2478 5504 : if( psDefn->PCS != KvUserDefined && psDefn->ProjCode == KvUserDefined )
2479 : {
2480 : int nZone;
2481 0 : int nGCS = psDefn->GCS;
2482 :
2483 0 : const int nMapSys = GTIFPCSToMapSys( psDefn->PCS, &nGCS, &nZone );
2484 0 : if( nMapSys != KvUserDefined )
2485 : {
2486 0 : psDefn->ProjCode = (short) GTIFMapSysToProj( nMapSys, nZone );
2487 0 : psDefn->GCS = (short) nGCS;
2488 : }
2489 : }
2490 :
2491 : /* -------------------------------------------------------------------- */
2492 : /* If the Proj_ code is specified directly, use that. */
2493 : /* -------------------------------------------------------------------- */
2494 5504 : if( psDefn->ProjCode == KvUserDefined )
2495 2557 : GTIFKeyGetSSHORT(psGTIF, ProjectionGeoKey, &(psDefn->ProjCode));
2496 :
2497 5504 : if( psDefn->ProjCode != KvUserDefined )
2498 : {
2499 : /*
2500 : * We have an underlying projection transformation value. Look
2501 : * this up. For a PCS of ``WGS 84 / UTM 11'' the transformation
2502 : * would be Transverse Mercator, with a particular set of options.
2503 : * The nProjTRFCode itself would correspond to the name
2504 : * ``UTM zone 11N'', and doesn't include datum info.
2505 : */
2506 2960 : GTIFGetProjTRFInfoEx( psGTIF->pj_context,
2507 2960 : psDefn->ProjCode, NULL, &(psDefn->Projection),
2508 2960 : psDefn->ProjParm );
2509 :
2510 : /*
2511 : * Set the GeoTIFF identity of the parameters.
2512 : */
2513 2960 : psDefn->CTProjection = (short)
2514 2960 : EPSGProjMethodToCTProjMethod( psDefn->Projection, FALSE );
2515 :
2516 2960 : SetGTParamIds( EPSGProjMethodToCTProjMethod(psDefn->Projection, TRUE),
2517 2960 : psDefn->Projection,
2518 2960 : psDefn->ProjParmId, NULL);
2519 2960 : psDefn->nParms = 7;
2520 : }
2521 :
2522 : /* -------------------------------------------------------------------- */
2523 : /* Try to get a GCS. If found, it will override any implied by */
2524 : /* the PCS. */
2525 : /* -------------------------------------------------------------------- */
2526 5504 : GTIFKeyGetSSHORT(psGTIF, GeographicTypeGeoKey, &(psDefn->GCS));
2527 5504 : if( psDefn->GCS < 1 || psDefn->GCS >= KvUserDefined )
2528 405 : psDefn->GCS = KvUserDefined;
2529 :
2530 : /* -------------------------------------------------------------------- */
2531 : /* Derive the datum, and prime meridian from the GCS. */
2532 : /* -------------------------------------------------------------------- */
2533 5504 : if( psDefn->GCS != KvUserDefined )
2534 : {
2535 5099 : GTIFGetGCSInfoEx( psGTIF->pj_context,
2536 5099 : psDefn->GCS, NULL, &(psDefn->Datum), &(psDefn->PM),
2537 : &(psDefn->UOMAngle) );
2538 : }
2539 :
2540 : /* -------------------------------------------------------------------- */
2541 : /* Handle the GCS angular units. GeogAngularUnitsGeoKey */
2542 : /* overrides the GCS or PCS setting. */
2543 : /* -------------------------------------------------------------------- */
2544 5504 : GTIFKeyGetSSHORT(psGTIF, GeogAngularUnitsGeoKey, &(psDefn->UOMAngle));
2545 5504 : if( psDefn->UOMAngle != KvUserDefined )
2546 : {
2547 5478 : GTIFGetUOMAngleInfoEx( psGTIF->pj_context,
2548 5478 : psDefn->UOMAngle, NULL,
2549 : &(psDefn->UOMAngleInDegrees) );
2550 : }
2551 :
2552 : /* -------------------------------------------------------------------- */
2553 : /* Check for a datum setting, and then use the datum to derive */
2554 : /* an ellipsoid. */
2555 : /* -------------------------------------------------------------------- */
2556 5504 : GTIFKeyGetSSHORT(psGTIF, GeogGeodeticDatumGeoKey, &(psDefn->Datum));
2557 :
2558 5504 : if( psDefn->Datum != KvUserDefined )
2559 : {
2560 5101 : GTIFGetDatumInfoEx( psGTIF->pj_context,
2561 5101 : psDefn->Datum, NULL, &(psDefn->Ellipsoid) );
2562 : }
2563 :
2564 : /* -------------------------------------------------------------------- */
2565 : /* Check for an explicit ellipsoid. Use the ellipsoid to */
2566 : /* derive the ellipsoid characteristics, if possible. */
2567 : /* -------------------------------------------------------------------- */
2568 5504 : GTIFKeyGetSSHORT(psGTIF, GeogEllipsoidGeoKey, &(psDefn->Ellipsoid));
2569 :
2570 5504 : if( psDefn->Ellipsoid != KvUserDefined )
2571 : {
2572 5104 : GTIFGetEllipsoidInfoEx( psGTIF->pj_context,
2573 5104 : psDefn->Ellipsoid, NULL,
2574 : &(psDefn->SemiMajor), &(psDefn->SemiMinor) );
2575 : }
2576 :
2577 : /* -------------------------------------------------------------------- */
2578 : /* Check for overridden ellipsoid parameters. It would be nice */
2579 : /* to warn if they conflict with provided information, but for */
2580 : /* now we just override. */
2581 : /* -------------------------------------------------------------------- */
2582 5504 : CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMajorAxisGeoKey, &(psDefn->SemiMajor), 0, 1 ));
2583 5504 : CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMinorAxisGeoKey, &(psDefn->SemiMinor), 0, 1 ));
2584 :
2585 : double dfInvFlattening;
2586 5504 : if( GTIFKeyGetDOUBLE(psGTIF, GeogInvFlatteningGeoKey, &dfInvFlattening,
2587 : 0, 1 ) == 1 )
2588 : {
2589 1826 : if( dfInvFlattening != 0.0 )
2590 1826 : psDefn->SemiMinor =
2591 1826 : psDefn->SemiMajor * (1 - 1.0/dfInvFlattening);
2592 : else
2593 0 : psDefn->SemiMinor = psDefn->SemiMajor;
2594 : }
2595 :
2596 : /* -------------------------------------------------------------------- */
2597 : /* Get the prime meridian info. */
2598 : /* -------------------------------------------------------------------- */
2599 5504 : GTIFKeyGetSSHORT(psGTIF, GeogPrimeMeridianGeoKey, &(psDefn->PM));
2600 :
2601 5504 : if( psDefn->PM != KvUserDefined )
2602 : {
2603 5100 : GTIFGetPMInfoEx( psGTIF->pj_context,
2604 5100 : psDefn->PM, NULL, &(psDefn->PMLongToGreenwich) );
2605 : }
2606 : else
2607 : {
2608 404 : CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogPrimeMeridianLongGeoKey,
2609 : &(psDefn->PMLongToGreenwich), 0, 1 ));
2610 :
2611 404 : psDefn->PMLongToGreenwich =
2612 404 : GTIFAngleToDD( psDefn->PMLongToGreenwich,
2613 404 : psDefn->UOMAngle );
2614 : }
2615 :
2616 : /* -------------------------------------------------------------------- */
2617 : /* Get the TOWGS84 parameters. */
2618 : /* -------------------------------------------------------------------- */
2619 : #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2620 5504 : psDefn->TOWGS84Count =
2621 5504 : (short)GTIFKeyGetDOUBLE(psGTIF, GeogTOWGS84GeoKey, psDefn->TOWGS84, 0, 7 );
2622 : #endif
2623 :
2624 : /* -------------------------------------------------------------------- */
2625 : /* Have the projection units of measure been overridden? We */
2626 : /* should likely be doing something about angular units too, */
2627 : /* but these are very rarely not decimal degrees for actual */
2628 : /* file coordinates. */
2629 : /* -------------------------------------------------------------------- */
2630 5504 : GTIFKeyGetSSHORT(psGTIF,ProjLinearUnitsGeoKey,&(psDefn->UOMLength));
2631 :
2632 5504 : if( psDefn->UOMLength != KvUserDefined )
2633 : {
2634 3360 : GTIFGetUOMLengthInfoEx( psGTIF->pj_context,
2635 3360 : psDefn->UOMLength, NULL,
2636 : &(psDefn->UOMLengthInMeters) );
2637 : }
2638 : else
2639 : {
2640 2144 : CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF,ProjLinearUnitSizeGeoKey,&(psDefn->UOMLengthInMeters),0,1));
2641 : }
2642 :
2643 : /* -------------------------------------------------------------------- */
2644 : /* Handle a variety of user defined transform types. */
2645 : /* -------------------------------------------------------------------- */
2646 5504 : if( GTIFKeyGetSSHORT(psGTIF,ProjCoordTransGeoKey,
2647 : &(psDefn->CTProjection)) == 1)
2648 : {
2649 74 : GTIFFetchProjParms( psGTIF, psDefn );
2650 : }
2651 :
2652 : /* -------------------------------------------------------------------- */
2653 : /* Try to set the zoned map system information. */
2654 : /* -------------------------------------------------------------------- */
2655 5504 : psDefn->MapSys = GTIFProjToMapSys( psDefn->ProjCode, &(psDefn->Zone) );
2656 :
2657 : /* -------------------------------------------------------------------- */
2658 : /* If this is UTM, and we were unable to extract the projection */
2659 : /* parameters from the database just set them directly now, */
2660 : /* since it's pretty easy, and a common case. */
2661 : /* -------------------------------------------------------------------- */
2662 5504 : if( (psDefn->MapSys == MapSys_UTM_North
2663 2697 : || psDefn->MapSys == MapSys_UTM_South)
2664 2810 : && psDefn->CTProjection == KvUserDefined )
2665 : {
2666 0 : psDefn->CTProjection = CT_TransverseMercator;
2667 0 : psDefn->nParms = 7;
2668 0 : psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
2669 0 : psDefn->ProjParm[0] = 0.0;
2670 :
2671 0 : psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
2672 0 : psDefn->ProjParm[1] = psDefn->Zone*6 - 183.0;
2673 :
2674 0 : psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
2675 0 : psDefn->ProjParm[4] = 0.9996;
2676 :
2677 0 : psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2678 0 : psDefn->ProjParm[5] = 500000.0;
2679 :
2680 0 : psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2681 :
2682 0 : if( psDefn->MapSys == MapSys_UTM_North )
2683 0 : psDefn->ProjParm[6] = 0.0;
2684 : else
2685 0 : psDefn->ProjParm[6] = 10000000.0;
2686 : }
2687 :
2688 5504 : return TRUE;
2689 : }
2690 :
2691 : /************************************************************************/
2692 : /* GTIFDecToDMS() */
2693 : /* */
2694 : /* Convenient function to translate decimal degrees to DMS */
2695 : /* format for reporting to a user. */
2696 : /************************************************************************/
2697 :
2698 0 : const char *GTIFDecToDMS( double dfAngle, const char * pszAxis,
2699 : int nPrecision )
2700 : {
2701 0 : if( !(dfAngle >= -360 && dfAngle <= 360) )
2702 0 : return "";
2703 :
2704 0 : double dfRound = 0.5/60;
2705 0 : for( int i = 0; i < nPrecision; i++ )
2706 0 : dfRound = dfRound * 0.1;
2707 :
2708 0 : int nDegrees = (int) ABS(dfAngle);
2709 0 : int nMinutes = (int) ((ABS(dfAngle) - nDegrees) * 60 + dfRound);
2710 0 : if( nMinutes == 60 )
2711 : {
2712 0 : nDegrees ++;
2713 0 : nMinutes = 0;
2714 : }
2715 0 : const double dfSeconds = ABS((ABS(dfAngle) * 3600 - nDegrees*3600 - nMinutes*60));
2716 :
2717 0 : const char *pszHemisphere = NULL;
2718 0 : if( EQUAL(pszAxis,"Long") && dfAngle < 0.0 )
2719 0 : pszHemisphere = "W";
2720 0 : else if( EQUAL(pszAxis,"Long") )
2721 0 : pszHemisphere = "E";
2722 0 : else if( dfAngle < 0.0 )
2723 0 : pszHemisphere = "S";
2724 : else
2725 0 : pszHemisphere = "N";
2726 :
2727 : char szFormat[30];
2728 0 : snprintf(szFormat, sizeof(szFormat), "%%3dd%%2d\'%%%d.%df\"%s",
2729 : nPrecision+3, nPrecision, pszHemisphere);
2730 : static char szBuffer[50];
2731 0 : snprintf(
2732 : szBuffer, sizeof(szBuffer), szFormat, nDegrees, nMinutes, dfSeconds);
2733 :
2734 0 : return szBuffer;
2735 : }
2736 :
2737 : /************************************************************************/
2738 : /* GTIFPrintDefn() */
2739 : /* */
2740 : /* Report the contents of a GTIFDefn structure ... mostly for */
2741 : /* debugging. */
2742 : /************************************************************************/
2743 :
2744 0 : void GTIFPrintDefnEx( GTIF *psGTIF, GTIFDefn * psDefn, FILE * fp )
2745 :
2746 : {
2747 0 : GTIFGetPROJContext(psGTIF, TRUE, NULL);
2748 :
2749 : /* -------------------------------------------------------------------- */
2750 : /* Do we have anything to report? */
2751 : /* -------------------------------------------------------------------- */
2752 0 : if( !psDefn->DefnSet )
2753 : {
2754 0 : fprintf( fp, "No GeoKeys found.\n" );
2755 0 : return;
2756 : }
2757 :
2758 : /* -------------------------------------------------------------------- */
2759 : /* Get the PCS name if possible. */
2760 : /* -------------------------------------------------------------------- */
2761 0 : if( psDefn->PCS != KvUserDefined )
2762 : {
2763 0 : char *pszPCSName = NULL;
2764 :
2765 0 : if( psGTIF->pj_context )
2766 : {
2767 0 : GTIFGetPCSInfoEx( psGTIF->pj_context,
2768 0 : psDefn->PCS, &pszPCSName, NULL, NULL, NULL );
2769 : }
2770 0 : if( pszPCSName == NULL )
2771 0 : pszPCSName = CPLStrdup("name unknown");
2772 :
2773 0 : fprintf( fp, "PCS = %d (%s)\n", psDefn->PCS, pszPCSName );
2774 0 : CPLFree( pszPCSName );
2775 : }
2776 :
2777 : /* -------------------------------------------------------------------- */
2778 : /* Dump the projection code if possible. */
2779 : /* -------------------------------------------------------------------- */
2780 0 : if( psDefn->ProjCode != KvUserDefined )
2781 : {
2782 0 : char *pszTRFName = NULL;
2783 :
2784 0 : if( psGTIF->pj_context )
2785 : {
2786 0 : GTIFGetProjTRFInfoEx( psGTIF->pj_context,
2787 0 : psDefn->ProjCode, &pszTRFName, NULL, NULL );
2788 : }
2789 0 : if( pszTRFName == NULL )
2790 0 : pszTRFName = CPLStrdup("");
2791 :
2792 0 : fprintf( fp, "Projection = %d (%s)\n",
2793 0 : psDefn->ProjCode, pszTRFName );
2794 :
2795 0 : CPLFree( pszTRFName );
2796 : }
2797 :
2798 : /* -------------------------------------------------------------------- */
2799 : /* Try to dump the projection method name, and parameters if possible.*/
2800 : /* -------------------------------------------------------------------- */
2801 0 : if( psDefn->CTProjection != KvUserDefined )
2802 : {
2803 : const char *pszProjectionMethodName =
2804 0 : GTIFValueNameEx(psGTIF,
2805 : ProjCoordTransGeoKey,
2806 0 : psDefn->CTProjection);
2807 :
2808 0 : if( pszProjectionMethodName == NULL )
2809 0 : pszProjectionMethodName = "(unknown)";
2810 :
2811 0 : fprintf( fp, "Projection Method: %s\n", pszProjectionMethodName );
2812 :
2813 0 : for( int i = 0; i < psDefn->nParms; i++ )
2814 : {
2815 0 : if( psDefn->ProjParmId[i] == 0 )
2816 0 : continue;
2817 :
2818 0 : char* pszName = GTIFKeyName((geokey_t) psDefn->ProjParmId[i]);
2819 0 : if( pszName == NULL )
2820 0 : pszName = "(unknown)";
2821 :
2822 0 : if( i < 4 )
2823 : {
2824 : char *pszAxisName;
2825 :
2826 0 : if( strstr(pszName,"Long") != NULL )
2827 0 : pszAxisName = "Long";
2828 0 : else if( strstr(pszName,"Lat") != NULL )
2829 0 : pszAxisName = "Lat";
2830 : else
2831 0 : pszAxisName = "?";
2832 :
2833 0 : fprintf( fp, " %s: %f (%s)\n",
2834 : pszName, psDefn->ProjParm[i],
2835 : GTIFDecToDMS( psDefn->ProjParm[i], pszAxisName, 2 ) );
2836 : }
2837 0 : else if( i == 4 )
2838 0 : fprintf( fp, " %s: %f\n", pszName, psDefn->ProjParm[i] );
2839 : else
2840 0 : fprintf( fp, " %s: %f m\n", pszName, psDefn->ProjParm[i] );
2841 : }
2842 : }
2843 :
2844 : /* -------------------------------------------------------------------- */
2845 : /* Report the GCS name, and number. */
2846 : /* -------------------------------------------------------------------- */
2847 0 : if( psDefn->GCS != KvUserDefined )
2848 : {
2849 0 : char *pszName = NULL;
2850 :
2851 0 : if( psGTIF->pj_context )
2852 : {
2853 0 : GTIFGetGCSInfoEx( psGTIF->pj_context,
2854 0 : psDefn->GCS, &pszName, NULL, NULL, NULL );
2855 : }
2856 0 : if( pszName == NULL )
2857 0 : pszName = CPLStrdup("(unknown)");
2858 :
2859 0 : fprintf( fp, "GCS: %d/%s\n", psDefn->GCS, pszName );
2860 0 : CPLFree( pszName );
2861 : }
2862 :
2863 : /* -------------------------------------------------------------------- */
2864 : /* Report the datum name. */
2865 : /* -------------------------------------------------------------------- */
2866 0 : if( psDefn->Datum != KvUserDefined )
2867 : {
2868 0 : char *pszName = NULL;
2869 :
2870 0 : if( psGTIF->pj_context )
2871 : {
2872 0 : GTIFGetDatumInfoEx( psGTIF->pj_context,
2873 0 : psDefn->Datum, &pszName, NULL );
2874 : }
2875 0 : if( pszName == NULL )
2876 0 : pszName = CPLStrdup("(unknown)");
2877 :
2878 0 : fprintf( fp, "Datum: %d/%s\n", psDefn->Datum, pszName );
2879 0 : CPLFree( pszName );
2880 : }
2881 :
2882 : /* -------------------------------------------------------------------- */
2883 : /* Report the ellipsoid. */
2884 : /* -------------------------------------------------------------------- */
2885 0 : if( psDefn->Ellipsoid != KvUserDefined )
2886 : {
2887 0 : char *pszName = NULL;
2888 :
2889 0 : if( psGTIF->pj_context )
2890 : {
2891 0 : GTIFGetEllipsoidInfoEx( psGTIF->pj_context,
2892 0 : psDefn->Ellipsoid, &pszName, NULL, NULL );
2893 : }
2894 0 : if( pszName == NULL )
2895 0 : pszName = CPLStrdup("(unknown)");
2896 :
2897 0 : fprintf( fp, "Ellipsoid: %d/%s (%.2f,%.2f)\n",
2898 0 : psDefn->Ellipsoid, pszName,
2899 : psDefn->SemiMajor, psDefn->SemiMinor );
2900 0 : CPLFree( pszName );
2901 : }
2902 :
2903 : /* -------------------------------------------------------------------- */
2904 : /* Report the prime meridian. */
2905 : /* -------------------------------------------------------------------- */
2906 0 : if( psDefn->PM != KvUserDefined )
2907 : {
2908 0 : char *pszName = NULL;
2909 :
2910 0 : if( psGTIF->pj_context )
2911 : {
2912 0 : GTIFGetPMInfoEx( psGTIF->pj_context,
2913 0 : psDefn->PM, &pszName, NULL );
2914 : }
2915 :
2916 0 : if( pszName == NULL )
2917 0 : pszName = CPLStrdup("(unknown)");
2918 :
2919 0 : fprintf( fp, "Prime Meridian: %d/%s (%f/%s)\n",
2920 0 : psDefn->PM, pszName,
2921 : psDefn->PMLongToGreenwich,
2922 : GTIFDecToDMS( psDefn->PMLongToGreenwich, "Long", 2 ) );
2923 0 : CPLFree( pszName );
2924 : }
2925 :
2926 : /* -------------------------------------------------------------------- */
2927 : /* Report TOWGS84 parameters. */
2928 : /* -------------------------------------------------------------------- */
2929 : #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2930 0 : if( psDefn->TOWGS84Count > 0 )
2931 : {
2932 0 : fprintf( fp, "TOWGS84: " );
2933 :
2934 0 : for( int i = 0; i < psDefn->TOWGS84Count; i++ )
2935 : {
2936 0 : if( i > 0 )
2937 0 : fprintf( fp, "," );
2938 0 : fprintf( fp, "%g", psDefn->TOWGS84[i] );
2939 : }
2940 :
2941 0 : fprintf( fp, "\n" );
2942 : }
2943 : #endif
2944 :
2945 : /* -------------------------------------------------------------------- */
2946 : /* Report the projection units of measure (currently just */
2947 : /* linear). */
2948 : /* -------------------------------------------------------------------- */
2949 0 : if( psDefn->UOMLength != KvUserDefined )
2950 : {
2951 0 : char *pszName = NULL;
2952 :
2953 0 : if( psGTIF->pj_context )
2954 : {
2955 0 : GTIFGetUOMLengthInfoEx(
2956 0 : psGTIF->pj_context, psDefn->UOMLength, &pszName, NULL );
2957 : }
2958 0 : if( pszName == NULL )
2959 0 : pszName = CPLStrdup( "(unknown)" );
2960 :
2961 0 : fprintf( fp, "Projection Linear Units: %d/%s (%fm)\n",
2962 0 : psDefn->UOMLength, pszName, psDefn->UOMLengthInMeters );
2963 0 : CPLFree( pszName );
2964 : }
2965 : else
2966 : {
2967 0 : fprintf( fp, "Projection Linear Units: User-Defined (%fm)\n",
2968 : psDefn->UOMLengthInMeters );
2969 : }
2970 : }
2971 :
2972 0 : void GTIFPrintDefn( GTIFDefn * psDefn, FILE * fp )
2973 : {
2974 0 : GTIF *psGTIF = GTIFNew(NULL);
2975 0 : if( psGTIF )
2976 : {
2977 0 : GTIFPrintDefnEx(psGTIF, psDefn, fp);
2978 0 : GTIFFree(psGTIF);
2979 : }
2980 0 : }
2981 :
2982 : /************************************************************************/
2983 : /* GTIFFreeMemory() */
2984 : /* */
2985 : /* Externally visible function to free memory allocated within */
2986 : /* geo_normalize.c. */
2987 : /************************************************************************/
2988 :
2989 31957 : void GTIFFreeMemory( char * pMemory )
2990 :
2991 : {
2992 31957 : if( pMemory != NULL )
2993 31943 : VSIFree( pMemory );
2994 31957 : }
2995 :
2996 : /************************************************************************/
2997 : /* GTIFAllocDefn() */
2998 : /* */
2999 : /* This allocates a GTIF structure in such a way that the */
3000 : /* calling application doesn't need to know the size and */
3001 : /* initializes it appropriately. */
3002 : /************************************************************************/
3003 :
3004 7785 : GTIFDefn *GTIFAllocDefn()
3005 : {
3006 7785 : return (GTIFDefn *) CPLCalloc(sizeof(GTIFDefn),1);
3007 : }
3008 :
3009 : /************************************************************************/
3010 : /* GTIFFreeDefn() */
3011 : /* */
3012 : /* Free a GTIF structure allocated by GTIFAllocDefn(). */
3013 : /************************************************************************/
3014 :
3015 7785 : void GTIFFreeDefn( GTIFDefn *defn )
3016 : {
3017 7785 : VSIFree( defn );
3018 7785 : }
3019 :
3020 : /************************************************************************/
3021 : /* GTIFAttachPROJContext() */
3022 : /* */
3023 : /* Attach an existing PROJ context to the GTIF handle, but */
3024 : /* ownership of the context remains to the caller. */
3025 : /************************************************************************/
3026 :
3027 21192 : void GTIFAttachPROJContext( GTIF *psGTIF, void* pjContext )
3028 : {
3029 21192 : if( psGTIF->own_pj_context )
3030 : {
3031 0 : proj_context_destroy(psGTIF->pj_context);
3032 : }
3033 21192 : psGTIF->own_pj_context = FALSE;
3034 21192 : psGTIF->pj_context = (PJ_CONTEXT*) pjContext;
3035 21192 : }
3036 :
3037 : /************************************************************************/
3038 : /* GTIFGetPROJContext() */
3039 : /* */
3040 : /* Return the PROJ context attached to the GTIF handle. */
3041 : /* If it has not yet been instantiated and instantiateIfNeeded=TRUE*/
3042 : /* then, it will be instantiated (and owned by GTIF handle). */
3043 : /************************************************************************/
3044 :
3045 18479 : void *GTIFGetPROJContext( GTIF *psGTIF, int instantiateIfNeeded,
3046 : int* out_gtif_own_pj_context )
3047 : {
3048 18479 : if( psGTIF->pj_context || !instantiateIfNeeded )
3049 : {
3050 18479 : if( out_gtif_own_pj_context )
3051 : {
3052 0 : *out_gtif_own_pj_context = psGTIF->own_pj_context;
3053 : }
3054 18479 : return psGTIF->pj_context;
3055 : }
3056 0 : psGTIF->pj_context = proj_context_create();
3057 0 : psGTIF->own_pj_context = psGTIF->pj_context != NULL;
3058 0 : if( out_gtif_own_pj_context )
3059 : {
3060 0 : *out_gtif_own_pj_context = psGTIF->own_pj_context;
3061 : }
3062 0 : return psGTIF->pj_context;
3063 : }
3064 :
3065 :
3066 0 : void GTIFDeaccessCSV( void )
3067 : {
3068 : /* No operation */
3069 0 : }
3070 :
3071 : #ifndef GDAL_COMPILATION
3072 : void SetCSVFilenameHook( const char *(*CSVFileOverride)(const char *) )
3073 : {
3074 : (void)CSVFileOverride;
3075 : /* No operation */
3076 : }
3077 : #endif
|
