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