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