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