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