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