00001 /****************************************************************************** 00002 * $Id: geo_normalize_c-source.html 385 2001-03-05 04:58:33Z warmerda $ 00003 * 00004 * Project: libgeotiff 00005 * Purpose: Code to normalize PCS and other composite codes in a GeoTIFF file. 00006 * Author: Frank Warmerdam, warmerda@home.com 00007 * 00008 ****************************************************************************** 00009 * Copyright (c) 1999, Frank Warmerdam 00010 * 00011 * Permission is hereby granted, free of charge, to any person obtaining a 00012 * copy of this software and associated documentation files (the "Software"), 00013 * to deal in the Software without restriction, including without limitation 00014 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 00015 * and/or sell copies of the Software, and to permit persons to whom the 00016 * Software is furnished to do so, subject to the following conditions: 00017 * 00018 * The above copyright notice and this permission notice shall be included 00019 * in all copies or substantial portions of the Software. 00020 * 00021 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 00022 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 00023 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 00024 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 00025 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 00026 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 00027 * DEALINGS IN THE SOFTWARE. 00028 ****************************************************************************** 00029 * 00030 * $Log$ 00030 * Revision 1.2 2001/03/05 04:58:33 warmerda 00030 * updated 00030 * 00031 * Revision 1.24 2001/03/05 03:26:29 warmerda 00032 * fixed memory leaks in GTIFPrintDefn() 00033 * 00034 * Revision 1.23 2001/02/23 13:49:48 warmerda 00035 * Fixed GTIFPrintDefn() to use fprintf( fp ), instead of printf(). 00036 * 00037 * Revision 1.22 2000/10/13 14:30:57 warmerda 00038 * fixed LCC parm order when parameters read directly from geotiff file 00039 * 00040 * Revision 1.21 2000/09/15 19:30:14 warmerda 00041 * report units of linear proj parms 00042 * 00043 * Revision 1.20 2000/09/15 18:21:07 warmerda 00044 * Fixed order of parameters for LCC 2SP. When parameters 00045 * were read from EPSG CSV files the standard parallels and origin 00046 * were mixed up. This affects alot of state plane zones! 00047 * 00048 * Revision 1.19 2000/06/09 14:05:43 warmerda 00049 * added default knowledge of NAD27/NAD83/WGS72/WGS84 00050 * 00051 * Revision 1.18 1999/12/10 21:28:12 warmerda 00052 * fixed Stereographic to look for ProjCenterLat/Long 00053 * 00054 * Revision 1.17 1999/12/10 20:06:58 warmerda 00055 * fixed up scale geokey used for a couple of projections 00056 * 00057 * Revision 1.16 1999/12/10 19:50:21 warmerda 00058 * Added EquidistantConic support, fixed return of StdParallel2GeoKey for 00059 * LCC2, and Albers. 00060 * 00061 * Revision 1.15 1999/12/10 19:39:26 warmerda 00062 * Fixed bug setting the false northing for files with 00063 * ProjCenterNorthingGeoKey set in GTIFGetDefn(). 00064 * 00065 * Revision 1.14 1999/09/17 14:58:37 warmerda 00066 * Added ProjRectifiedGridAngleGeoKey(3096) and support for it's 00067 * use with Oblique Mercator in geo_normalize.c. 00068 * 00069 * Revision 1.13 1999/09/17 00:55:26 warmerda 00070 * added GTIFGetUOMAngleInfo(), and UOMAngle in GTIFDefn 00071 * 00072 * Revision 1.12 1999/09/15 18:51:31 warmerda 00073 * Map 9808 to TM South Oriented, not TM Modified Alaska. 00074 * 00075 * Revision 1.11 1999/09/15 16:44:06 warmerda 00076 * Change meter to metre to match EPSG database in GTIFGetUOMLengthInfo() 00077 * shortcut. 00078 * 00079 * Revision 1.10 1999/09/15 16:35:15 warmerda 00080 * Fixed the fractions of second handling properly in GTIFAngleStringToDD(). 00081 * 00082 * Revision 1.9 1999/09/15 14:24:17 warmerda 00083 * Fixed serious bug in geo_normalize.c with translation of 00084 * DD.MMSSsss values. Return value was seriously off if any 00085 * fraction of a second was included in the string. 00086 * 00087 * Revision 1.8 1999/07/13 03:12:52 warmerda 00088 * Make scale a parameter of CT_Stereographic. 00089 * 00090 * Revision 1.7 1999/05/04 03:13:22 warmerda 00091 * fixed a serious bug in parsing DMSmmss.sss values, and a bug in forming DMS strings 00092 * 00093 * Revision 1.6 1999/05/03 17:50:31 warmerda 00094 * avoid warnings on IRIX 00095 * 00096 * Revision 1.5 1999/04/28 20:04:51 warmerda 00097 * Added doxygen style documentation. 00098 * Use GTIFPCSToMapSys() and related functions to partially normalize 00099 * projections when we don't have the CSV files. 00100 * 00101 * Revision 1.4 1999/03/18 21:34:59 geotiff 00102 * added GTIFDecToDMS 00103 * 00104 * Revision 1.3 1999/03/17 19:53:15 geotiff 00105 * sys includes moved to cpl_serv.h 00106 * 00107 * Revision 1.2 1999/03/10 18:24:06 geotiff 00108 * corrected to use int' 00109 * 00110 * Revision 1.1 1999/03/09 15:57:04 geotiff 00111 * New 00112 * 00113 * Revision 1.4 1999/03/03 02:29:38 warmerda 00114 * Define PI if not already defined. 00115 * 00116 * Revision 1.3 1999/03/02 21:10:57 warmerda 00117 * added lots of projections 00118 * 00119 * Revision 1.2 1999/02/24 16:24:15 warmerda 00120 * Continuing to evolve 00121 * 00122 * Revision 1.1 1999/02/22 18:51:08 warmerda 00123 * New 00124 * 00125 */ 00126 00127 #include "cpl_csv.h" 00128 #include "geotiff.h" 00129 #include "xtiffio.h" 00130 #include "geovalues.h" 00131 #include "geo_normalize.h" 00132 00133 #ifndef KvUserDefined 00134 # define KvUserDefined 32767 00135 #endif 00136 00137 #ifndef PI 00138 # define PI 3.14159265358979323846 00139 #endif 00140 00141 /************************************************************************/ 00142 /* GTIFGetPCSInfo() */ 00143 /************************************************************************/ 00144 00145 int GTIFGetPCSInfo( int nPCSCode, char **ppszEPSGName, 00146 short *pnUOMLengthCode, short *pnUOMAngleCode, 00147 short *pnGeogCS, short *pnTRFCode ) 00148 00149 { 00150 char **papszRecord; 00151 char szSearchKey[24]; 00152 const char *pszFilename = CSVFilename( "horiz_cs.csv" ); 00153 00154 /* -------------------------------------------------------------------- */ 00155 /* Search the units database for this unit. If we don't find */ 00156 /* it return failure. */ 00157 /* -------------------------------------------------------------------- */ 00158 sprintf( szSearchKey, "%d", nPCSCode ); 00159 papszRecord = CSVScanFileByName( pszFilename, "HORIZCS_CODE", 00160 szSearchKey, CC_Integer ); 00161 00162 if( papszRecord == NULL ) 00163 return FALSE; 00164 00165 /* -------------------------------------------------------------------- */ 00166 /* Get the name, if requested. */ 00167 /* -------------------------------------------------------------------- */ 00168 if( ppszEPSGName != NULL ) 00169 { 00170 *ppszEPSGName = 00171 CPLStrdup( CSLGetField( papszRecord, 00172 CSVGetFileFieldId(pszFilename, 00173 "HORIZCS_EPSG_NAME") )); 00174 } 00175 00176 /* -------------------------------------------------------------------- */ 00177 /* Get the UOM Length code, if requested. */ 00178 /* -------------------------------------------------------------------- */ 00179 if( pnUOMLengthCode != NULL ) 00180 { 00181 const char *pszValue; 00182 00183 pszValue = 00184 CSLGetField( papszRecord, 00185 CSVGetFileFieldId(pszFilename,"UOM_LENGTH_CODE")); 00186 if( atoi(pszValue) > 0 ) 00187 *pnUOMLengthCode = atoi(pszValue); 00188 else 00189 *pnUOMLengthCode = KvUserDefined; 00190 } 00191 00192 /* -------------------------------------------------------------------- */ 00193 /* Get the UOM Angle code, if requested. */ 00194 /* -------------------------------------------------------------------- */ 00195 if( pnUOMAngleCode != NULL ) 00196 { 00197 const char *pszValue; 00198 00199 pszValue = 00200 CSLGetField( papszRecord, 00201 CSVGetFileFieldId(pszFilename,"UOM_ANGLE_CODE") ); 00202 00203 if( atoi(pszValue) > 0 ) 00204 *pnUOMAngleCode = atoi(pszValue); 00205 else 00206 *pnUOMAngleCode = KvUserDefined; 00207 } 00208 00209 /* -------------------------------------------------------------------- */ 00210 /* Get the GeogCS (Datum with PM) code, if requested. */ 00211 /* -------------------------------------------------------------------- */ 00212 if( pnGeogCS != NULL ) 00213 { 00214 const char *pszValue; 00215 00216 pszValue = 00217 CSLGetField( papszRecord, 00218 CSVGetFileFieldId(pszFilename,"SOURCE_GEOGCS_CODE") ); 00219 if( atoi(pszValue) > 0 ) 00220 *pnGeogCS = atoi(pszValue); 00221 else 00222 *pnGeogCS = KvUserDefined; 00223 } 00224 00225 /* -------------------------------------------------------------------- */ 00226 /* Get the GeogCS (Datum with PM) code, if requested. */ 00227 /* -------------------------------------------------------------------- */ 00228 if( pnTRFCode != NULL ) 00229 { 00230 const char *pszValue; 00231 00232 pszValue = 00233 CSLGetField( papszRecord, 00234 CSVGetFileFieldId(pszFilename,"PROJECTION_TRF_CODE")); 00235 00236 00237 if( atoi(pszValue) > 0 ) 00238 *pnTRFCode = atoi(pszValue); 00239 else 00240 *pnTRFCode = KvUserDefined; 00241 } 00242 00243 return TRUE; 00244 } 00245 00246 /************************************************************************/ 00247 /* GTIFAngleToDD() */ 00248 /* */ 00249 /* Convert a numeric angle to decimal degress. */ 00250 /************************************************************************/ 00251 00252 double GTIFAngleToDD( double dfAngle, int nUOMAngle ) 00253 00254 { 00255 if( nUOMAngle == 9110 ) /* DDD.MMSSsss */ 00256 { 00257 char szAngleString[32]; 00258 00259 sprintf( szAngleString, "%12.7f", dfAngle ); 00260 dfAngle = GTIFAngleStringToDD( szAngleString, nUOMAngle ); 00261 } 00262 else 00263 { 00264 double dfInDegrees; 00265 00266 GTIFGetUOMAngleInfo( nUOMAngle, NULL, &dfInDegrees ); 00267 dfAngle = dfAngle * dfInDegrees; 00268 } 00269 00270 return( dfAngle ); 00271 } 00272 00273 /************************************************************************/ 00274 /* GTIFAngleStringToDD() */ 00275 /* */ 00276 /* Convert an angle in the specified units to decimal degrees. */ 00277 /************************************************************************/ 00278 00279 double GTIFAngleStringToDD( const char * pszAngle, int nUOMAngle ) 00280 00281 { 00282 double dfAngle; 00283 00284 if( nUOMAngle == 9110 ) /* DDD.MMSSsss */ 00285 { 00286 char *pszDecimal; 00287 00288 dfAngle = ABS(atoi(pszAngle)); 00289 pszDecimal = strchr(pszAngle,'.'); 00290 if( pszDecimal != NULL && strlen(pszDecimal) > 1 ) 00291 { 00292 char szMinutes[3]; 00293 char szSeconds[64]; 00294 00295 szMinutes[0] = pszDecimal[1]; 00296 if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' ) 00297 szMinutes[1] = pszDecimal[2]; 00298 else 00299 szMinutes[1] = '0'; 00300 00301 szMinutes[2] = '\0'; 00302 dfAngle += atoi(szMinutes) / 60.0; 00303 00304 if( strlen(pszDecimal) > 3 ) 00305 { 00306 szSeconds[0] = pszDecimal[3]; 00307 if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' ) 00308 { 00309 szSeconds[1] = pszDecimal[4]; 00310 szSeconds[2] = '.'; 00311 strcpy( szSeconds+3, pszDecimal + 5 ); 00312 } 00313 else 00314 { 00315 szSeconds[1] = '0'; 00316 szSeconds[2] = '\0'; 00317 } 00318 dfAngle += atof(szSeconds) / 3600.0; 00319 } 00320 } 00321 00322 if( pszAngle[0] == '-' ) 00323 dfAngle *= -1; 00324 } 00325 else if( nUOMAngle == 9105 || nUOMAngle == 9106 ) /* grad */ 00326 { 00327 dfAngle = 180 * (atof(pszAngle ) / 200); 00328 } 00329 else if( nUOMAngle == 9101 ) /* radians */ 00330 { 00331 dfAngle = 180 * (atof(pszAngle ) / PI); 00332 } 00333 else if( nUOMAngle == 9103 ) /* arc-minute */ 00334 { 00335 dfAngle = atof(pszAngle) / 60; 00336 } 00337 else if( nUOMAngle == 9104 ) /* arc-second */ 00338 { 00339 dfAngle = atof(pszAngle) / 3600; 00340 } 00341 else /* decimal degrees ... some cases missing but seeminly never used */ 00342 { 00343 CPLAssert( nUOMAngle == 9102 || nUOMAngle == KvUserDefined 00344 || nUOMAngle == 0 ); 00345 00346 dfAngle = atof(pszAngle ); 00347 } 00348 00349 return( dfAngle ); 00350 } 00351 00352 /************************************************************************/ 00353 /* GTIFGetGCSInfo() */ 00354 /* */ 00355 /* Fetch the datum, and prime meridian related to a particular */ 00356 /* GCS. */ 00357 /************************************************************************/ 00358 00359 int GTIFGetGCSInfo( int nGCSCode, char ** ppszName, 00360 short * pnDatum, short * pnPM, short *pnUOMAngle ) 00361 00362 { 00363 char szSearchKey[24]; 00364 int nDatum, nPM, nUOMAngle; 00365 00366 /* -------------------------------------------------------------------- */ 00367 /* Search the database for the corresponding datum code. */ 00368 /* -------------------------------------------------------------------- */ 00369 sprintf( szSearchKey, "%d", nGCSCode ); 00370 00371 nDatum = atoi(CSVGetField( CSVFilename("horiz_cs.csv" ), 00372 "HORIZCS_CODE", szSearchKey, CC_Integer, 00373 "GEOD_DATUM_CODE" ) ); 00374 00375 /* -------------------------------------------------------------------- */ 00376 /* Handle some "well known" GCS codes directly if the table */ 00377 /* wasn't found. */ 00378 /* -------------------------------------------------------------------- */ 00379 if( nDatum < 1 ) 00380 { 00381 const char * pszName = NULL; 00382 nPM = PM_Greenwich; 00383 nUOMAngle = Angular_DMS_Hemisphere; 00384 if( nGCSCode == GCS_NAD27 ) 00385 { 00386 nDatum = Datum_North_American_Datum_1927; 00387 pszName = "NAD27"; 00388 } 00389 else if( nGCSCode == GCS_NAD83 ) 00390 { 00391 nDatum = Datum_North_American_Datum_1983; 00392 pszName = "NAD83"; 00393 } 00394 else if( nGCSCode == GCS_WGS_84 ) 00395 { 00396 nDatum = Datum_WGS84; 00397 pszName = "WGS 84"; 00398 } 00399 else if( nGCSCode == GCS_WGS_72 ) 00400 { 00401 nDatum = Datum_WGS72; 00402 pszName = "WGS 82"; 00403 } 00404 else 00405 return FALSE; 00406 00407 if( ppszName != NULL ) 00408 *ppszName = CPLStrdup( pszName ); 00409 if( pnDatum != NULL ) 00410 *pnDatum = nDatum; 00411 if( pnPM != NULL ) 00412 *pnPM = nPM; 00413 if( pnUOMAngle != NULL ) 00414 *pnUOMAngle = nUOMAngle; 00415 00416 return TRUE; 00417 } 00418 00419 /* -------------------------------------------------------------------- */ 00420 /* Get the PM. */ 00421 /* -------------------------------------------------------------------- */ 00422 if( pnDatum != NULL ) 00423 *pnDatum = nDatum; 00424 00425 nPM = atoi(CSVGetField( CSVFilename("horiz_cs.csv" ), 00426 "HORIZCS_CODE", szSearchKey, CC_Integer, 00427 "PRIME_MERIDIAN_CODE" ) ); 00428 00429 if( nPM < 1 ) 00430 return FALSE; 00431 00432 if( pnPM != NULL ) 00433 *pnPM = nPM; 00434 00435 /* -------------------------------------------------------------------- */ 00436 /* Get the angular units. */ 00437 /* -------------------------------------------------------------------- */ 00438 nUOMAngle = atoi(CSVGetField( CSVFilename("horiz_cs.csv" ), 00439 "HORIZCS_CODE", szSearchKey, CC_Integer, 00440 "UOM_ANGLE_CODE" ) ); 00441 00442 if( nUOMAngle < 1 ) 00443 return FALSE; 00444 00445 if( pnUOMAngle != NULL ) 00446 *pnUOMAngle = nUOMAngle; 00447 00448 /* -------------------------------------------------------------------- */ 00449 /* Get the name, if requested. */ 00450 /* -------------------------------------------------------------------- */ 00451 if( ppszName != NULL ) 00452 *ppszName = 00453 CPLStrdup(CSVGetField( CSVFilename("horiz_cs.csv" ), 00454 "HORIZCS_CODE", szSearchKey, CC_Integer, 00455 "HORIZCS_EPSG_NAME" )); 00456 00457 return( TRUE ); 00458 } 00459 00460 /************************************************************************/ 00461 /* GTIFGetEllipsoidInfo() */ 00462 /* */ 00463 /* Fetch info about an ellipsoid. Axes are always returned in */ 00464 /* meters. SemiMajor computed based on inverse flattening */ 00465 /* where that is provided. */ 00466 /************************************************************************/ 00467 00468 int GTIFGetEllipsoidInfo( int nEllipseCode, char ** ppszName, 00469 double * pdfSemiMajor, double * pdfSemiMinor ) 00470 00471 { 00472 char szSearchKey[24]; 00473 double dfSemiMajor, dfToMeters = 1.0; 00474 int nUOMLength; 00475 00476 /* -------------------------------------------------------------------- */ 00477 /* Get the semi major axis. */ 00478 /* -------------------------------------------------------------------- */ 00479 sprintf( szSearchKey, "%d", nEllipseCode ); 00480 00481 dfSemiMajor = 00482 atof(CSVGetField( CSVFilename("ellipsoid.csv" ), 00483 "ELLIPSOID_CODE", szSearchKey, CC_Integer, 00484 "SEMI_MAJOR_AXIS" ) ); 00485 00486 /* -------------------------------------------------------------------- */ 00487 /* Try some well known ellipsoids. */ 00488 /* -------------------------------------------------------------------- */ 00489 if( dfSemiMajor == 0.0 ) 00490 { 00491 double dfInvFlattening, dfSemiMinor; 00492 const char *pszName = NULL; 00493 00494 if( nEllipseCode == Ellipse_Clarke_1866 ) 00495 { 00496 pszName = "Clarke 1866"; 00497 dfSemiMajor = 6378206.4; 00498 dfSemiMinor = 6356583.8; 00499 dfInvFlattening = 0.0; 00500 } 00501 else if( nEllipseCode == Ellipse_GRS_1980 ) 00502 { 00503 pszName = "GRS 1980"; 00504 dfSemiMajor = 6378137.0; 00505 dfSemiMinor = 0.0; 00506 dfInvFlattening = 298.257222101; 00507 } 00508 else if( nEllipseCode == Ellipse_WGS_84 ) 00509 { 00510 pszName = "WGS 84"; 00511 dfSemiMajor = 6378137.0; 00512 dfSemiMinor = 0.0; 00513 dfInvFlattening = 298.257223563; 00514 } 00515 else if( nEllipseCode == 7043 ) 00516 { 00517 pszName = "WGS 72"; 00518 dfSemiMajor = 6378135.0; 00519 dfSemiMinor = 0.0; 00520 dfInvFlattening = 298.26; 00521 } 00522 else 00523 return FALSE; 00524 00525 if( dfSemiMinor == 0.0 ) 00526 dfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening); 00527 00528 if( pdfSemiMinor != NULL ) 00529 *pdfSemiMinor = dfSemiMinor; 00530 if( pdfSemiMajor != NULL ) 00531 *pdfSemiMajor = dfSemiMajor; 00532 if( ppszName != NULL ) 00533 *ppszName = CPLStrdup( pszName ); 00534 00535 return TRUE; 00536 } 00537 00538 /* -------------------------------------------------------------------- */ 00539 /* Get the translation factor into meters. */ 00540 /* -------------------------------------------------------------------- */ 00541 nUOMLength = atoi(CSVGetField( CSVFilename("ellipsoid.csv" ), 00542 "ELLIPSOID_CODE", szSearchKey, CC_Integer, 00543 "UOM_LENGTH_CODE" )); 00544 GTIFGetUOMLengthInfo( nUOMLength, NULL, &dfToMeters ); 00545 00546 dfSemiMajor *= dfToMeters; 00547 00548 if( pdfSemiMajor != NULL ) 00549 *pdfSemiMajor = dfSemiMajor; 00550 00551 /* -------------------------------------------------------------------- */ 00552 /* Get the semi-minor if requested. If the Semi-minor axis */ 00553 /* isn't available, compute it based on the inverse flattening. */ 00554 /* -------------------------------------------------------------------- */ 00555 if( pdfSemiMinor != NULL ) 00556 { 00557 *pdfSemiMinor = 00558 atof(CSVGetField( CSVFilename("ellipsoid.csv" ), 00559 "ELLIPSOID_CODE", szSearchKey, CC_Integer, 00560 "SEMI_MINOR_AXIS" )) * dfToMeters; 00561 00562 if( *pdfSemiMinor == 0.0 ) 00563 { 00564 double dfInvFlattening; 00565 00566 dfInvFlattening = 00567 atof(CSVGetField( CSVFilename("ellipsoid.csv" ), 00568 "ELLIPSOID_CODE", szSearchKey, CC_Integer, 00569 "INV_FLATTENING" )); 00570 *pdfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening); 00571 } 00572 } 00573 00574 /* -------------------------------------------------------------------- */ 00575 /* Get the name, if requested. */ 00576 /* -------------------------------------------------------------------- */ 00577 if( ppszName != NULL ) 00578 *ppszName = 00579 CPLStrdup(CSVGetField( CSVFilename("ellipsoid.csv" ), 00580 "ELLIPSOID_CODE", szSearchKey, CC_Integer, 00581 "ELLIPSOID_EPSG_NAME" )); 00582 00583 return( TRUE ); 00584 } 00585 00586 /************************************************************************/ 00587 /* GTIFGetPMInfo() */ 00588 /* */ 00589 /* Get the offset between a given prime meridian and Greenwich */ 00590 /* in degrees. */ 00591 /************************************************************************/ 00592 00593 int GTIFGetPMInfo( int nPMCode, char ** ppszName, double *pdfOffset ) 00594 00595 { 00596 char szSearchKey[24]; 00597 int nUOMAngle; 00598 00599 /* -------------------------------------------------------------------- */ 00600 /* Use a special short cut for Greenwich, since it is so common. */ 00601 /* -------------------------------------------------------------------- */ 00602 if( nPMCode == PM_Greenwich ) 00603 { 00604 if( pdfOffset != NULL ) 00605 *pdfOffset = 0.0; 00606 if( ppszName != NULL ) 00607 *ppszName = CPLStrdup( "Greenwich" ); 00608 return TRUE; 00609 } 00610 00611 /* -------------------------------------------------------------------- */ 00612 /* Search the database for the corresponding datum code. */ 00613 /* -------------------------------------------------------------------- */ 00614 sprintf( szSearchKey, "%d", nPMCode ); 00615 00616 nUOMAngle = 00617 atoi(CSVGetField( CSVFilename("p_meridian.csv" ), 00618 "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer, 00619 "UOM_ANGLE_CODE" ) ); 00620 if( nUOMAngle < 1 ) 00621 return FALSE; 00622 00623 /* -------------------------------------------------------------------- */ 00624 /* Get the PM offset. */ 00625 /* -------------------------------------------------------------------- */ 00626 if( pdfOffset != NULL ) 00627 { 00628 *pdfOffset = 00629 GTIFAngleStringToDD( 00630 CSVGetField( CSVFilename("p_meridian.csv" ), 00631 "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer, 00632 "GREENWICH_LONGITUDE" ), 00633 nUOMAngle ); 00634 } 00635 00636 /* -------------------------------------------------------------------- */ 00637 /* Get the name, if requested. */ 00638 /* -------------------------------------------------------------------- */ 00639 if( ppszName != NULL ) 00640 *ppszName = 00641 CPLStrdup( 00642 CSVGetField( CSVFilename("p_meridian.csv" ), 00643 "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer, 00644 "PRIME_MERID_EPSG_NAME" )); 00645 00646 return( TRUE ); 00647 } 00648 00649 /************************************************************************/ 00650 /* GTIFGetDatumInfo() */ 00651 /* */ 00652 /* Fetch the ellipsoid, and name for a datum. */ 00653 /************************************************************************/ 00654 00655 int GTIFGetDatumInfo( int nDatumCode, char ** ppszName, short * pnEllipsoid ) 00656 00657 { 00658 char szSearchKey[24]; 00659 int nEllipsoid; 00660 00661 /* -------------------------------------------------------------------- */ 00662 /* Search the database for the corresponding datum code. */ 00663 /* -------------------------------------------------------------------- */ 00664 sprintf( szSearchKey, "%d", nDatumCode ); 00665 00666 nEllipsoid = atoi(CSVGetField( CSVFilename("geod_datum.csv" ), 00667 "GEOD_DATUM_CODE", szSearchKey, CC_Integer, 00668 "ELLIPSOID_CODE" ) ); 00669 00670 if( pnEllipsoid != NULL ) 00671 *pnEllipsoid = nEllipsoid; 00672 00673 /* -------------------------------------------------------------------- */ 00674 /* Handle a few built-in datums. */ 00675 /* -------------------------------------------------------------------- */ 00676 if( nEllipsoid < 1 ) 00677 { 00678 const char *pszName = NULL; 00679 00680 if( nDatumCode == Datum_North_American_Datum_1927 ) 00681 { 00682 nEllipsoid = Ellipse_Clarke_1866; 00683 pszName = "North American Datum 1927"; 00684 } 00685 else if( nDatumCode == Datum_North_American_Datum_1983 ) 00686 { 00687 nEllipsoid = Ellipse_GRS_1980; 00688 pszName = "North American Datum 1983"; 00689 } 00690 else if( nDatumCode == Datum_WGS84 ) 00691 { 00692 nEllipsoid = Ellipse_WGS_84; 00693 pszName = "World Geodetic System 1984"; 00694 } 00695 else if( nDatumCode == Datum_WGS72 ) 00696 { 00697 nEllipsoid = 7043; /* WGS7 */ 00698 pszName = "World Geodetic System 1972"; 00699 } 00700 else 00701 return FALSE; 00702 00703 if( pnEllipsoid != NULL ) 00704 *pnEllipsoid = nEllipsoid; 00705 00706 if( ppszName != NULL ) 00707 *ppszName = CPLStrdup( pszName ); 00708 00709 return TRUE; 00710 } 00711 00712 /* -------------------------------------------------------------------- */ 00713 /* Get the name, if requested. */ 00714 /* -------------------------------------------------------------------- */ 00715 if( ppszName != NULL ) 00716 *ppszName = 00717 CPLStrdup(CSVGetField( CSVFilename("geod_datum.csv" ), 00718 "GEOD_DATUM_CODE", szSearchKey, CC_Integer, 00719 "GEOD_DAT_EPSG_NAME" )); 00720 00721 return( TRUE ); 00722 } 00723 00724 00725 /************************************************************************/ 00726 /* GTIFGetUOMLengthInfo() */ 00727 /* */ 00728 /* Note: This function should eventually also know how to */ 00729 /* lookup length aliases in the UOM_LE_ALIAS table. */ 00730 /************************************************************************/ 00731 00732 int GTIFGetUOMLengthInfo( int nUOMLengthCode, 00733 char **ppszUOMName, 00734 double * pdfInMeters ) 00735 00736 { 00737 char **papszUnitsRecord; 00738 char szSearchKey[24]; 00739 int iNameField; 00740 00741 /* -------------------------------------------------------------------- */ 00742 /* We short cut meter to save work in the most common case. */ 00743 /* -------------------------------------------------------------------- */ 00744 if( nUOMLengthCode == 9001 ) 00745 { 00746 if( ppszUOMName != NULL ) 00747 *ppszUOMName = CPLStrdup( "metre" ); 00748 if( pdfInMeters != NULL ) 00749 *pdfInMeters = 1.0; 00750 00751 return TRUE; 00752 } 00753 00754 /* -------------------------------------------------------------------- */ 00755 /* Search the units database for this unit. If we don't find */ 00756 /* it return failure. */ 00757 /* -------------------------------------------------------------------- */ 00758 sprintf( szSearchKey, "%d", nUOMLengthCode ); 00759 papszUnitsRecord = 00760 CSVScanFileByName( CSVFilename( "uom_length.csv" ), 00761 "UOM_LENGTH_CODE", szSearchKey, CC_Integer ); 00762 00763 if( papszUnitsRecord == NULL ) 00764 return FALSE; 00765 00766 /* -------------------------------------------------------------------- */ 00767 /* Get the name, if requested. */ 00768 /* -------------------------------------------------------------------- */ 00769 if( ppszUOMName != NULL ) 00770 { 00771 iNameField = CSVGetFileFieldId( CSVFilename( "uom_length.csv" ), 00772 "UNIT_OF_MEAS_EPSG_NAME" ); 00773 *ppszUOMName = CPLStrdup( CSLGetField(papszUnitsRecord, iNameField) ); 00774 } 00775 00776 /* -------------------------------------------------------------------- */ 00777 /* Get the A and B factor fields, and create the multiplicative */ 00778 /* factor. */ 00779 /* -------------------------------------------------------------------- */ 00780 if( pdfInMeters != NULL ) 00781 { 00782 int iBFactorField, iCFactorField; 00783 00784 iBFactorField = CSVGetFileFieldId( CSVFilename( "uom_length.csv" ), 00785 "FACTOR_B" ); 00786 iCFactorField = CSVGetFileFieldId( CSVFilename( "uom_length.csv" ), 00787 "FACTOR_C" ); 00788 00789 if( atof(CSLGetField(papszUnitsRecord, iCFactorField)) > 0.0 ) 00790 *pdfInMeters = atof(CSLGetField(papszUnitsRecord, iBFactorField)) 00791 / atof(CSLGetField(papszUnitsRecord, iCFactorField)); 00792 else 00793 *pdfInMeters = 0.0; 00794 } 00795 00796 return( TRUE ); 00797 } 00798 00799 /************************************************************************/ 00800 /* GTIFGetUOMAngleInfo() */ 00801 /************************************************************************/ 00802 00803 int GTIFGetUOMAngleInfo( int nUOMLengthCode, 00804 char **ppszUOMName, 00805 double * pdfInDegrees ) 00806 00807 { 00808 const char *pszUOMName = NULL; 00809 double dfInDegrees = 0.0; 00810 00811 switch( nUOMLengthCode ) 00812 { 00813 case 9101: 00814 pszUOMName = "radian"; 00815 dfInDegrees = 180.0 / 3.14159265358979; 00816 break; 00817 00818 case 9102: 00819 case 9107: 00820 case 9108: 00821 case 9110: 00822 pszUOMName = "degree"; 00823 dfInDegrees = 1.0; 00824 break; 00825 00826 case 9103: 00827 pszUOMName = "arc-minute"; 00828 dfInDegrees = 1 / 60.0; 00829 break; 00830 00831 case 9104: 00832 pszUOMName = "arc-second"; 00833 dfInDegrees = 1 / 3600.0; 00834 break; 00835 00836 case 9105: 00837 pszUOMName = "grad"; 00838 dfInDegrees = 180.0 / 200.0; 00839 break; 00840 00841 case 9106: 00842 pszUOMName = "gon"; 00843 dfInDegrees = 180.0 / 200.0; 00844 break; 00845 00846 case 9109: 00847 pszUOMName = "microradian"; 00848 dfInDegrees = 180.0 / (3.14159265358979 * 1000000.0); 00849 break; 00850 } 00851 00852 if( ppszUOMName != NULL ) 00853 { 00854 if( pszUOMName != NULL ) 00855 *ppszUOMName = CPLStrdup( pszUOMName ); 00856 else 00857 *ppszUOMName = NULL; 00858 } 00859 00860 if( pdfInDegrees != NULL ) 00861 *pdfInDegrees = dfInDegrees; 00862 00863 return( TRUE ); 00864 } 00865 00866 /************************************************************************/ 00867 /* EPSGProjMethodToCTProjMethod() */ 00868 /* */ 00869 /* Convert between the EPSG enumeration for projection methods, */ 00870 /* and the GeoTIFF CT codes. */ 00871 /************************************************************************/ 00872 00873 static int EPSGProjMethodToCTProjMethod( int nEPSG ) 00874 00875 { 00876 /* see trf_method.csv for list of EPSG codes */ 00877 00878 switch( nEPSG ) 00879 { 00880 case 9801: 00881 return( CT_LambertConfConic_1SP ); 00882 00883 case 9802: 00884 return( CT_LambertConfConic_2SP ); 00885 00886 case 9803: 00887 return( CT_LambertConfConic_2SP ); /* Belgian variant not supported */ 00888 00889 case 9804: 00890 return( CT_Mercator ); /* 1SP and 2SP not differentiated */ 00891 00892 case 9805: 00893 return( CT_Mercator ); /* 1SP and 2SP not differentiated */ 00894 00895 case 9806: 00896 return( CT_CassiniSoldner ); 00897 00898 case 9807: 00899 return( CT_TransverseMercator ); 00900 00901 case 9808: 00902 return( CT_TransvMercator_SouthOriented ); 00903 00904 case 9809: 00905 return( CT_ObliqueStereographic ); 00906 00907 case 9810: 00908 return( CT_PolarStereographic ); 00909 00910 case 9811: 00911 return( CT_NewZealandMapGrid ); 00912 00913 case 9812: 00914 return( CT_ObliqueMercator ); /* is hotine actually different? */ 00915 00916 case 9813: 00917 return( CT_ObliqueMercator_Laborde ); 00918 00919 case 9814: 00920 return( CT_ObliqueMercator_Rosenmund ); /* swiss */ 00921 00922 case 9815: 00923 return( CT_ObliqueMercator ); 00924 00925 case 9816: /* tunesia mining grid has no counterpart */ 00926 return( KvUserDefined ); 00927 } 00928 00929 return( KvUserDefined ); 00930 } 00931 00932 /************************************************************************/ 00933 /* GTIFGetProjTRFInfo() */ 00934 /* */ 00935 /* Transform a PROJECTION_TRF_CODE into a projection method, */ 00936 /* and a set of parameters. The parameters identify will */ 00937 /* depend on the returned method, but they will all have been */ 00938 /* normalized into degrees and meters. */ 00939 /************************************************************************/ 00940 00941 int GTIFGetProjTRFInfo( int nProjTRFCode, 00942 char **ppszProjTRFName, 00943 short * pnProjMethod, 00944 double * padfProjParms ) 00945 00946 { 00947 int nProjMethod, nUOMLinear, nUOMAngle, i; 00948 double adfProjParms[7], dfInMeters; 00949 char szTRFCode[16]; 00950 00951 /* -------------------------------------------------------------------- */ 00952 /* Get the proj method. If this fails to return a meaningful */ 00953 /* number, then the whole function fails. */ 00954 /* -------------------------------------------------------------------- */ 00955 sprintf( szTRFCode, "%d", nProjTRFCode ); 00956 nProjMethod = 00957 atoi( CSVGetField( CSVFilename( "trf_nonpolynomial.csv" ), 00958 "COORD_TRF_CODE", szTRFCode, CC_Integer, 00959 "COORD_TRF_METHOD_CODE" ) ); 00960 if( nProjMethod == 0 ) 00961 return FALSE; 00962 00963 /* -------------------------------------------------------------------- */ 00964 /* Get the linear, and angular (geog) units for the projection */ 00965 /* parameters. */ 00966 /* -------------------------------------------------------------------- */ 00967 nUOMLinear = 00968 atoi( CSVGetField( CSVFilename( "trf_nonpolynomial.csv" ), 00969 "COORD_TRF_CODE", szTRFCode, CC_Integer, 00970 "UOM_LENGTH_CODE" ) ); 00971 00972 if( !GTIFGetUOMLengthInfo( nUOMLinear, NULL, &dfInMeters ) ) 00973 dfInMeters = 1.0; 00974 00975 nUOMAngle = 00976 atoi( CSVGetField( CSVFilename( "trf_nonpolynomial.csv" ), 00977 "COORD_TRF_CODE", szTRFCode, CC_Integer, 00978 "UOM_ANGLE_CODE" ) ); 00979 00980 /* -------------------------------------------------------------------- */ 00981 /* Get the parameters for this projection. For the time being */ 00982 /* I am assuming the first four parameters are angles, the */ 00983 /* fifth is unitless (normally scale), and the remainder are */ 00984 /* linear measures. This works fine for the existing */ 00985 /* projections, but is a pretty fragile approach. */ 00986 /* -------------------------------------------------------------------- */ 00987 00988 for( i = 0; i < 7; i++ ) 00989 { 00990 char szID[16]; 00991 const char *pszValue; 00992 00993 sprintf( szID, "PARAMETER_%d", i+1 ); 00994 00995 pszValue = CSVGetField( CSVFilename( "trf_nonpolynomial.csv" ), 00996 "COORD_TRF_CODE", szTRFCode,CC_Integer, szID ); 00997 if( i < 4 ) 00998 adfProjParms[i] = GTIFAngleStringToDD( pszValue, nUOMAngle ); 00999 else if( i > 4 ) 01000 adfProjParms[i] = atof(pszValue) * dfInMeters; 01001 else 01002 adfProjParms[i] = atof( pszValue ); 01003 } 01004 01005 /* -------------------------------------------------------------------- */ 01006 /* Get the name, if requested. */ 01007 /* -------------------------------------------------------------------- */ 01008 if( ppszProjTRFName != NULL ) 01009 { 01010 *ppszProjTRFName = 01011 CPLStrdup(CSVGetField( CSVFilename("trf_nonpolynomial.csv" ), 01012 "COORD_TRF_CODE", szTRFCode, CC_Integer, 01013 "COORD_TRF_EPSG_NAME" )); 01014 } 01015 01016 /* -------------------------------------------------------------------- */ 01017 /* Transfer requested data into passed variables. */ 01018 /* -------------------------------------------------------------------- */ 01019 if( pnProjMethod != NULL ) 01020 *pnProjMethod = nProjMethod; 01021 01022 if( padfProjParms != NULL ) 01023 { 01024 for( i = 0; i < 7; i++ ) 01025 padfProjParms[i] = adfProjParms[i]; 01026 } 01027 01028 return TRUE; 01029 } 01030 01031 /************************************************************************/ 01032 /* SetGTParmIds() */ 01033 /* */ 01034 /* This is hardcoded logic to set the GeoTIFF parmaeter */ 01035 /* identifiers for all the EPSG supported projections. As the */ 01036 /* trf_method.csv table grows with new projections, this code */ 01037 /* will need to be updated. */ 01038 /************************************************************************/ 01039 01040 static int SetGTParmIds( GTIFDefn * psDefn ) 01041 01042 { 01043 psDefn->nParms = 7; 01044 01045 switch( psDefn->CTProjection ) 01046 { 01047 case CT_CassiniSoldner: 01048 case CT_NewZealandMapGrid: 01049 psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey; 01050 psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey; 01051 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01052 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01053 return TRUE; 01054 01055 case CT_ObliqueMercator: 01056 psDefn->ProjParmId[0] = ProjCenterLatGeoKey; 01057 psDefn->ProjParmId[1] = ProjCenterLongGeoKey; 01058 psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey; 01059 psDefn->ProjParmId[3] = ProjRectifiedGridAngleGeoKey; 01060 psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey; 01061 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01062 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01063 return TRUE; 01064 01065 case CT_ObliqueMercator_Laborde: 01066 psDefn->ProjParmId[0] = ProjCenterLatGeoKey; 01067 psDefn->ProjParmId[1] = ProjCenterLongGeoKey; 01068 psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey; 01069 psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey; 01070 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01071 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01072 return TRUE; 01073 01074 case CT_LambertConfConic_1SP: 01075 case CT_Mercator: 01076 case CT_ObliqueStereographic: 01077 case CT_PolarStereographic: 01078 case CT_TransverseMercator: 01079 case CT_TransvMercator_SouthOriented: 01080 psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey; 01081 psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey; 01082 psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey; 01083 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01084 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01085 return TRUE; 01086 01087 case CT_LambertConfConic_2SP: 01088 psDefn->ProjParmId[0] = ProjFalseOriginLatGeoKey; 01089 psDefn->ProjParmId[1] = ProjFalseOriginLongGeoKey; 01090 psDefn->ProjParmId[2] = ProjStdParallel1GeoKey; 01091 psDefn->ProjParmId[3] = ProjStdParallel2GeoKey; 01092 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01093 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01094 return TRUE; 01095 01096 case CT_SwissObliqueCylindrical: 01097 psDefn->ProjParmId[0] = ProjCenterLatGeoKey; 01098 psDefn->ProjParmId[1] = ProjCenterLongGeoKey; 01099 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01100 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01101 return TRUE; 01102 01103 default: 01104 return( FALSE ); 01105 } 01106 } 01107 01108 /************************************************************************/ 01109 /* GTIFFetchProjParms() */ 01110 /* */ 01111 /* Fetch the projection parameters for a particular projection */ 01112 /* from a GeoTIFF file, and fill the GTIFDefn structure out */ 01113 /* with them. */ 01114 /************************************************************************/ 01115 01116 static void GTIFFetchProjParms( GTIF * psGTIF, GTIFDefn * psDefn ) 01117 01118 { 01119 double dfNatOriginLong, dfNatOriginLat, dfRectGridAngle; 01120 double dfFalseEasting, dfFalseNorthing, dfNatOriginScale; 01121 double dfStdParallel1, dfStdParallel2, dfAzimuth; 01122 01123 /* -------------------------------------------------------------------- */ 01124 /* Get the false easting, and northing if available. */ 01125 /* -------------------------------------------------------------------- */ 01126 if( !GTIFKeyGet(psGTIF, ProjFalseEastingGeoKey, &dfFalseEasting, 0, 1) 01127 && !GTIFKeyGet(psGTIF, ProjCenterEastingGeoKey, 01128 &dfFalseEasting, 0, 1) ) 01129 dfFalseEasting = 0.0; 01130 01131 if( !GTIFKeyGet(psGTIF, ProjFalseNorthingGeoKey, &dfFalseNorthing,0,1) 01132 && !GTIFKeyGet(psGTIF, ProjCenterNorthingGeoKey, 01133 &dfFalseNorthing, 0, 1) ) 01134 dfFalseNorthing = 0.0; 01135 01136 switch( psDefn->CTProjection ) 01137 { 01138 /* -------------------------------------------------------------------- */ 01139 case CT_Stereographic: 01140 /* -------------------------------------------------------------------- */ 01141 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01142 &dfNatOriginLong, 0, 1 ) == 0 01143 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01144 &dfNatOriginLong, 0, 1 ) == 0 01145 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01146 &dfNatOriginLong, 0, 1 ) == 0 ) 01147 dfNatOriginLong = 0.0; 01148 01149 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01150 &dfNatOriginLat, 0, 1 ) == 0 01151 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01152 &dfNatOriginLat, 0, 1 ) == 0 01153 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01154 &dfNatOriginLat, 0, 1 ) == 0 ) 01155 dfNatOriginLat = 0.0; 01156 01157 if( GTIFKeyGet(psGTIF, ProjScaleAtNatOriginGeoKey, 01158 &dfNatOriginScale, 0, 1 ) == 0 ) 01159 dfNatOriginScale = 1.0; 01160 01161 /* notdef: should transform to decimal degrees at this point */ 01162 01163 psDefn->ProjParm[0] = dfNatOriginLat; 01164 psDefn->ProjParmId[0] = ProjCenterLatGeoKey; 01165 psDefn->ProjParm[1] = dfNatOriginLong; 01166 psDefn->ProjParmId[1] = ProjCenterLongGeoKey; 01167 psDefn->ProjParm[4] = dfNatOriginScale; 01168 psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey; 01169 psDefn->ProjParm[5] = dfFalseEasting; 01170 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01171 psDefn->ProjParm[6] = dfFalseNorthing; 01172 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01173 01174 psDefn->nParms = 7; 01175 break; 01176 01177 /* -------------------------------------------------------------------- */ 01178 case CT_LambertConfConic_1SP: 01179 case CT_Mercator: 01180 case CT_ObliqueStereographic: 01181 case CT_TransverseMercator: 01182 case CT_TransvMercator_SouthOriented: 01183 /* -------------------------------------------------------------------- */ 01184 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01185 &dfNatOriginLong, 0, 1 ) == 0 01186 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01187 &dfNatOriginLong, 0, 1 ) == 0 01188 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01189 &dfNatOriginLong, 0, 1 ) == 0 ) 01190 dfNatOriginLong = 0.0; 01191 01192 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01193 &dfNatOriginLat, 0, 1 ) == 0 01194 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01195 &dfNatOriginLat, 0, 1 ) == 0 01196 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01197 &dfNatOriginLat, 0, 1 ) == 0 ) 01198 dfNatOriginLat = 0.0; 01199 01200 if( GTIFKeyGet(psGTIF, ProjScaleAtNatOriginGeoKey, 01201 &dfNatOriginScale, 0, 1 ) == 0 ) 01202 dfNatOriginScale = 1.0; 01203 01204 /* notdef: should transform to decimal degrees at this point */ 01205 01206 psDefn->ProjParm[0] = dfNatOriginLat; 01207 psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey; 01208 psDefn->ProjParm[1] = dfNatOriginLong; 01209 psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey; 01210 psDefn->ProjParm[4] = dfNatOriginScale; 01211 psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey; 01212 psDefn->ProjParm[5] = dfFalseEasting; 01213 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01214 psDefn->ProjParm[6] = dfFalseNorthing; 01215 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01216 01217 psDefn->nParms = 7; 01218 break; 01219 01220 /* -------------------------------------------------------------------- */ 01221 case CT_ObliqueMercator: /* hotine */ 01222 /* -------------------------------------------------------------------- */ 01223 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01224 &dfNatOriginLong, 0, 1 ) == 0 01225 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01226 &dfNatOriginLong, 0, 1 ) == 0 01227 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01228 &dfNatOriginLong, 0, 1 ) == 0 ) 01229 dfNatOriginLong = 0.0; 01230 01231 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01232 &dfNatOriginLat, 0, 1 ) == 0 01233 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01234 &dfNatOriginLat, 0, 1 ) == 0 01235 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01236 &dfNatOriginLat, 0, 1 ) == 0 ) 01237 dfNatOriginLat = 0.0; 01238 01239 if( GTIFKeyGet(psGTIF, ProjAzimuthAngleGeoKey, 01240 &dfAzimuth, 0, 1 ) == 0 ) 01241 dfAzimuth = 0.0; 01242 01243 if( GTIFKeyGet(psGTIF, ProjRectifiedGridAngleGeoKey, 01244 &dfRectGridAngle, 0, 1 ) == 0 ) 01245 dfRectGridAngle = 90.0; 01246 01247 if( GTIFKeyGet(psGTIF, ProjScaleAtNatOriginGeoKey, 01248 &dfNatOriginScale, 0, 1 ) == 0 01249 && GTIFKeyGet(psGTIF, ProjScaleAtCenterGeoKey, 01250 &dfNatOriginScale, 0, 1 ) == 0 ) 01251 dfNatOriginScale = 1.0; 01252 01253 /* notdef: should transform to decimal degrees at this point */ 01254 01255 psDefn->ProjParm[0] = dfNatOriginLat; 01256 psDefn->ProjParmId[0] = ProjCenterLatGeoKey; 01257 psDefn->ProjParm[1] = dfNatOriginLong; 01258 psDefn->ProjParmId[1] = ProjCenterLongGeoKey; 01259 psDefn->ProjParm[2] = dfAzimuth; 01260 psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey; 01261 psDefn->ProjParm[3] = dfRectGridAngle; 01262 psDefn->ProjParmId[3] = ProjRectifiedGridAngleGeoKey; 01263 psDefn->ProjParm[4] = dfNatOriginScale; 01264 psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey; 01265 psDefn->ProjParm[5] = dfFalseEasting; 01266 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01267 psDefn->ProjParm[6] = dfFalseNorthing; 01268 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01269 01270 psDefn->nParms = 7; 01271 break; 01272 01273 /* -------------------------------------------------------------------- */ 01274 case CT_CassiniSoldner: 01275 case CT_Polyconic: 01276 /* -------------------------------------------------------------------- */ 01277 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01278 &dfNatOriginLong, 0, 1 ) == 0 01279 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01280 &dfNatOriginLong, 0, 1 ) == 0 01281 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01282 &dfNatOriginLong, 0, 1 ) == 0 ) 01283 dfNatOriginLong = 0.0; 01284 01285 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01286 &dfNatOriginLat, 0, 1 ) == 0 01287 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01288 &dfNatOriginLat, 0, 1 ) == 0 01289 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01290 &dfNatOriginLat, 0, 1 ) == 0 ) 01291 dfNatOriginLat = 0.0; 01292 01293 if( GTIFKeyGet(psGTIF, ProjScaleAtNatOriginGeoKey, 01294 &dfNatOriginScale, 0, 1 ) == 0 01295 && GTIFKeyGet(psGTIF, ProjScaleAtCenterGeoKey, 01296 &dfNatOriginScale, 0, 1 ) == 0 ) 01297 dfNatOriginScale = 1.0; 01298 01299 /* notdef: should transform to decimal degrees at this point */ 01300 01301 psDefn->ProjParm[0] = dfNatOriginLat; 01302 psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey; 01303 psDefn->ProjParm[1] = dfNatOriginLong; 01304 psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey; 01305 psDefn->ProjParm[4] = dfNatOriginScale; 01306 psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey; 01307 psDefn->ProjParm[5] = dfFalseEasting; 01308 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01309 psDefn->ProjParm[6] = dfFalseNorthing; 01310 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01311 01312 psDefn->nParms = 7; 01313 break; 01314 01315 /* -------------------------------------------------------------------- */ 01316 case CT_AzimuthalEquidistant: 01317 case CT_MillerCylindrical: 01318 case CT_Equirectangular: 01319 case CT_Gnomonic: 01320 case CT_LambertAzimEqualArea: 01321 case CT_Orthographic: 01322 /* -------------------------------------------------------------------- */ 01323 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01324 &dfNatOriginLong, 0, 1 ) == 0 01325 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01326 &dfNatOriginLong, 0, 1 ) == 0 01327 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01328 &dfNatOriginLong, 0, 1 ) == 0 ) 01329 dfNatOriginLong = 0.0; 01330 01331 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01332 &dfNatOriginLat, 0, 1 ) == 0 01333 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01334 &dfNatOriginLat, 0, 1 ) == 0 01335 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01336 &dfNatOriginLat, 0, 1 ) == 0 ) 01337 dfNatOriginLat = 0.0; 01338 01339 /* notdef: should transform to decimal degrees at this point */ 01340 01341 psDefn->ProjParm[0] = dfNatOriginLat; 01342 psDefn->ProjParmId[0] = ProjCenterLatGeoKey; 01343 psDefn->ProjParm[1] = dfNatOriginLong; 01344 psDefn->ProjParmId[1] = ProjCenterLongGeoKey; 01345 psDefn->ProjParm[5] = dfFalseEasting; 01346 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01347 psDefn->ProjParm[6] = dfFalseNorthing; 01348 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01349 01350 psDefn->nParms = 7; 01351 break; 01352 01353 /* -------------------------------------------------------------------- */ 01354 case CT_Robinson: 01355 case CT_Sinusoidal: 01356 case CT_VanDerGrinten: 01357 /* -------------------------------------------------------------------- */ 01358 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01359 &dfNatOriginLong, 0, 1 ) == 0 01360 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01361 &dfNatOriginLong, 0, 1 ) == 0 01362 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01363 &dfNatOriginLong, 0, 1 ) == 0 ) 01364 dfNatOriginLong = 0.0; 01365 01366 /* notdef: should transform to decimal degrees at this point */ 01367 01368 psDefn->ProjParm[1] = dfNatOriginLong; 01369 psDefn->ProjParmId[1] = ProjCenterLongGeoKey; 01370 psDefn->ProjParm[5] = dfFalseEasting; 01371 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01372 psDefn->ProjParm[6] = dfFalseNorthing; 01373 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01374 01375 psDefn->nParms = 7; 01376 break; 01377 01378 /* -------------------------------------------------------------------- */ 01379 case CT_PolarStereographic: 01380 /* -------------------------------------------------------------------- */ 01381 if( GTIFKeyGet(psGTIF, ProjStraightVertPoleLongGeoKey, 01382 &dfNatOriginLong, 0, 1 ) == 0 01383 && GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01384 &dfNatOriginLong, 0, 1 ) == 0 01385 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01386 &dfNatOriginLong, 0, 1 ) == 0 01387 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01388 &dfNatOriginLong, 0, 1 ) == 0 ) 01389 dfNatOriginLong = 0.0; 01390 01391 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01392 &dfNatOriginLat, 0, 1 ) == 0 01393 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01394 &dfNatOriginLat, 0, 1 ) == 0 01395 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01396 &dfNatOriginLat, 0, 1 ) == 0 ) 01397 dfNatOriginLat = 0.0; 01398 01399 if( GTIFKeyGet(psGTIF, ProjScaleAtNatOriginGeoKey, 01400 &dfNatOriginScale, 0, 1 ) == 0 01401 && GTIFKeyGet(psGTIF, ProjScaleAtCenterGeoKey, 01402 &dfNatOriginScale, 0, 1 ) == 0 ) 01403 dfNatOriginScale = 1.0; 01404 01405 /* notdef: should transform to decimal degrees at this point */ 01406 01407 psDefn->ProjParm[0] = dfNatOriginLat; 01408 psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;; 01409 psDefn->ProjParm[1] = dfNatOriginLong; 01410 psDefn->ProjParmId[1] = ProjStraightVertPoleLongGeoKey; 01411 psDefn->ProjParm[4] = dfNatOriginScale; 01412 psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey; 01413 psDefn->ProjParm[5] = dfFalseEasting; 01414 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01415 psDefn->ProjParm[6] = dfFalseNorthing; 01416 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01417 01418 psDefn->nParms = 7; 01419 break; 01420 01421 /* -------------------------------------------------------------------- */ 01422 case CT_LambertConfConic_2SP: 01423 /* -------------------------------------------------------------------- */ 01424 if( GTIFKeyGet(psGTIF, ProjStdParallel1GeoKey, 01425 &dfStdParallel1, 0, 1 ) == 0 ) 01426 dfStdParallel1 = 0.0; 01427 01428 if( GTIFKeyGet(psGTIF, ProjStdParallel2GeoKey, 01429 &dfStdParallel2, 0, 1 ) == 0 ) 01430 dfStdParallel1 = 0.0; 01431 01432 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01433 &dfNatOriginLong, 0, 1 ) == 0 01434 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01435 &dfNatOriginLong, 0, 1 ) == 0 01436 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01437 &dfNatOriginLong, 0, 1 ) == 0 ) 01438 dfNatOriginLong = 0.0; 01439 01440 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01441 &dfNatOriginLat, 0, 1 ) == 0 01442 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01443 &dfNatOriginLat, 0, 1 ) == 0 01444 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01445 &dfNatOriginLat, 0, 1 ) == 0 ) 01446 dfNatOriginLat = 0.0; 01447 01448 /* notdef: should transform to decimal degrees at this point */ 01449 01450 psDefn->ProjParm[0] = dfNatOriginLat; 01451 psDefn->ProjParmId[0] = ProjFalseOriginLatGeoKey; 01452 psDefn->ProjParm[1] = dfNatOriginLong; 01453 psDefn->ProjParmId[1] = ProjFalseOriginLongGeoKey; 01454 psDefn->ProjParm[2] = dfStdParallel1; 01455 psDefn->ProjParmId[2] = ProjStdParallel1GeoKey; 01456 psDefn->ProjParm[3] = dfStdParallel2; 01457 psDefn->ProjParmId[3] = ProjStdParallel2GeoKey; 01458 psDefn->ProjParm[5] = dfFalseEasting; 01459 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01460 psDefn->ProjParm[6] = dfFalseNorthing; 01461 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01462 01463 psDefn->nParms = 7; 01464 break; 01465 01466 /* -------------------------------------------------------------------- */ 01467 case CT_AlbersEqualArea: 01468 case CT_EquidistantConic: 01469 /* -------------------------------------------------------------------- */ 01470 if( GTIFKeyGet(psGTIF, ProjStdParallel1GeoKey, 01471 &dfStdParallel1, 0, 1 ) == 0 ) 01472 dfStdParallel1 = 0.0; 01473 01474 if( GTIFKeyGet(psGTIF, ProjStdParallel2GeoKey, 01475 &dfStdParallel2, 0, 1 ) == 0 ) 01476 dfStdParallel1 = 0.0; 01477 01478 if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, 01479 &dfNatOriginLong, 0, 1 ) == 0 01480 && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, 01481 &dfNatOriginLong, 0, 1 ) == 0 01482 && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, 01483 &dfNatOriginLong, 0, 1 ) == 0 ) 01484 dfNatOriginLong = 0.0; 01485 01486 if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, 01487 &dfNatOriginLat, 0, 1 ) == 0 01488 && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, 01489 &dfNatOriginLat, 0, 1 ) == 0 01490 && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, 01491 &dfNatOriginLat, 0, 1 ) == 0 ) 01492 dfNatOriginLat = 0.0; 01493 01494 /* notdef: should transform to decimal degrees at this point */ 01495 01496 psDefn->ProjParm[0] = dfStdParallel1; 01497 psDefn->ProjParmId[0] = ProjStdParallel1GeoKey; 01498 psDefn->ProjParm[1] = dfStdParallel2; 01499 psDefn->ProjParmId[1] = ProjStdParallel2GeoKey; 01500 psDefn->ProjParm[2] = dfNatOriginLat; 01501 psDefn->ProjParmId[2] = ProjNatOriginLatGeoKey; 01502 psDefn->ProjParm[3] = dfNatOriginLong; 01503 psDefn->ProjParmId[3] = ProjNatOriginLongGeoKey; 01504 psDefn->ProjParm[5] = dfFalseEasting; 01505 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01506 psDefn->ProjParm[6] = dfFalseNorthing; 01507 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01508 01509 psDefn->nParms = 7; 01510 break; 01511 } 01512 } 01513 01514 /************************************************************************/ 01515 /* GTIFGetDefn() */ 01516 /************************************************************************/ 01517 01634 int GTIFGetDefn( GTIF * psGTIF, GTIFDefn * psDefn ) 01635 01636 { 01637 int i; 01638 short nGeogUOMLinear; 01639 01640 /* -------------------------------------------------------------------- */ 01641 /* Initially we default all the information we can. */ 01642 /* -------------------------------------------------------------------- */ 01643 psDefn->Model = KvUserDefined; 01644 psDefn->PCS = KvUserDefined; 01645 psDefn->GCS = KvUserDefined; 01646 psDefn->UOMLength = KvUserDefined; 01647 psDefn->UOMLengthInMeters = 1.0; 01648 psDefn->UOMAngle = KvUserDefined; 01649 psDefn->UOMAngleInDegrees = 1.0; 01650 psDefn->Datum = KvUserDefined; 01651 psDefn->Ellipsoid = KvUserDefined; 01652 psDefn->SemiMajor = 0.0; 01653 psDefn->SemiMinor = 0.0; 01654 psDefn->PM = KvUserDefined; 01655 psDefn->PMLongToGreenwich = 0.0; 01656 01657 psDefn->ProjCode = KvUserDefined; 01658 psDefn->Projection = KvUserDefined; 01659 psDefn->CTProjection = KvUserDefined; 01660 01661 psDefn->nParms = 0; 01662 for( i = 0; i < MAX_GTIF_PROJPARMS; i++ ) 01663 { 01664 psDefn->ProjParm[i] = 0.0; 01665 psDefn->ProjParmId[i] = 0; 01666 } 01667 01668 psDefn->MapSys = KvUserDefined; 01669 psDefn->Zone = 0; 01670 01671 /* -------------------------------------------------------------------- */ 01672 /* Try to get the overall model type. */ 01673 /* -------------------------------------------------------------------- */ 01674 GTIFKeyGet(psGTIF,GTModelTypeGeoKey,&(psDefn->Model),0,1); 01675 01676 /* -------------------------------------------------------------------- */ 01677 /* Extract the Geog units. */ 01678 /* -------------------------------------------------------------------- */ 01679 nGeogUOMLinear = 9001; /* Linear_Meter */ 01680 GTIFKeyGet(psGTIF, GeogLinearUnitsGeoKey, &nGeogUOMLinear, 0, 1 ); 01681 01682 /* -------------------------------------------------------------------- */ 01683 /* Try to get a PCS. */ 01684 /* -------------------------------------------------------------------- */ 01685 if( GTIFKeyGet(psGTIF,ProjectedCSTypeGeoKey, &(psDefn->PCS),0,1) == 1 01686 && psDefn->PCS != KvUserDefined ) 01687 { 01688 /* 01689 * Translate this into useful information. 01690 */ 01691 GTIFGetPCSInfo( psDefn->PCS, NULL, 01692 &(psDefn->UOMLength), &(psDefn->UOMAngle), 01693 &(psDefn->GCS), &(psDefn->ProjCode) ); 01694 } 01695 01696 /* -------------------------------------------------------------------- */ 01697 /* If we have the PCS code, but didn't find it in the CSV files */ 01698 /* (likely because we can't find them) we will try some ``jiffy */ 01699 /* rules'' for UTM and state plane. */ 01700 /* -------------------------------------------------------------------- */ 01701 if( psDefn->PCS != KvUserDefined && psDefn->ProjCode == KvUserDefined ) 01702 { 01703 int nMapSys, nZone; 01704 int nGCS = psDefn->GCS; 01705 01706 nMapSys = GTIFPCSToMapSys( psDefn->PCS, &nGCS, &nZone ); 01707 if( nMapSys != KvUserDefined ) 01708 { 01709 psDefn->ProjCode = GTIFMapSysToProj( nMapSys, nZone ); 01710 psDefn->GCS = nGCS; 01711 } 01712 } 01713 01714 /* -------------------------------------------------------------------- */ 01715 /* If the Proj_ code is specified directly, use that. */ 01716 /* -------------------------------------------------------------------- */ 01717 if( psDefn->ProjCode == KvUserDefined ) 01718 GTIFKeyGet(psGTIF, ProjectionGeoKey, &(psDefn->ProjCode), 0, 1 ); 01719 01720 if( psDefn->ProjCode != KvUserDefined ) 01721 { 01722 /* 01723 * We have an underlying projection transformation value. Look 01724 * this up. For a PCS of ``WGS 84 / UTM 11'' the transformation 01725 * would be Transverse Mercator, with a particular set of options. 01726 * The nProjTRFCode itself would correspond to the name 01727 * ``UTM zone 11N'', and doesn't include datum info. 01728 */ 01729 GTIFGetProjTRFInfo( psDefn->ProjCode, NULL, 01730 &(psDefn->Projection), 01731 psDefn->ProjParm ); 01732 psDefn->CTProjection = 01733 EPSGProjMethodToCTProjMethod( psDefn->Projection ); 01734 01735 /* 01736 * Set the GeoTIFF identity of the parameters. 01737 */ 01738 SetGTParmIds( psDefn ); 01739 } 01740 01741 /* -------------------------------------------------------------------- */ 01742 /* Try to get a GCS. If found, it will override any implied by */ 01743 /* the PCS. */ 01744 /* -------------------------------------------------------------------- */ 01745 GTIFKeyGet(psGTIF, GeographicTypeGeoKey, &(psDefn->GCS), 0, 1 ); 01746 01747 /* -------------------------------------------------------------------- */ 01748 /* Derive the datum, and prime meridian from the GCS. */ 01749 /* -------------------------------------------------------------------- */ 01750 if( psDefn->GCS != KvUserDefined ) 01751 { 01752 GTIFGetGCSInfo( psDefn->GCS, NULL, &(psDefn->Datum), &(psDefn->PM), 01753 &(psDefn->UOMAngle) ); 01754 } 01755 01756 /* -------------------------------------------------------------------- */ 01757 /* Handle the GCS angular units. GeogAngularUnitsGeoKey */ 01758 /* overrides the GCS or PCS setting. */ 01759 /* -------------------------------------------------------------------- */ 01760 GTIFKeyGet(psGTIF, GeogAngularUnitsGeoKey, &(psDefn->UOMAngle), 0, 1 ); 01761 if( psDefn->UOMAngle != KvUserDefined ) 01762 { 01763 GTIFGetUOMAngleInfo( psDefn->UOMAngle, NULL, 01764 &(psDefn->UOMAngleInDegrees) ); 01765 } 01766 01767 /* -------------------------------------------------------------------- */ 01768 /* Check for a datum setting, and then use the datum to derive */ 01769 /* an ellipsoid. */ 01770 /* -------------------------------------------------------------------- */ 01771 GTIFKeyGet(psGTIF, GeogGeodeticDatumGeoKey, &(psDefn->Datum), 0, 1 ); 01772 01773 if( psDefn->Datum != KvUserDefined ) 01774 { 01775 GTIFGetDatumInfo( psDefn->Datum, NULL, &(psDefn->Ellipsoid) ); 01776 } 01777 01778 /* -------------------------------------------------------------------- */ 01779 /* Check for an explicit ellipsoid. Use the ellipsoid to */ 01780 /* derive the ellipsoid characteristics, if possible. */ 01781 /* -------------------------------------------------------------------- */ 01782 GTIFKeyGet(psGTIF, GeogEllipsoidGeoKey, &(psDefn->Ellipsoid), 0, 1 ); 01783 01784 if( psDefn->Ellipsoid != KvUserDefined ) 01785 { 01786 GTIFGetEllipsoidInfo( psDefn->Ellipsoid, NULL, 01787 &(psDefn->SemiMajor), &(psDefn->SemiMinor) ); 01788 } 01789 01790 /* -------------------------------------------------------------------- */ 01791 /* Get the prime meridian info. */ 01792 /* -------------------------------------------------------------------- */ 01793 GTIFKeyGet(psGTIF, GeogPrimeMeridianGeoKey, &(psDefn->PM), 0, 1 ); 01794 01795 if( psDefn->PM != KvUserDefined ) 01796 { 01797 GTIFGetPMInfo( psDefn->PM, NULL, &(psDefn->PMLongToGreenwich) ); 01798 } 01799 else 01800 { 01801 GTIFKeyGet(psGTIF, GeogPrimeMeridianLongGeoKey, 01802 &(psDefn->PMLongToGreenwich), 0, 1 ); 01803 01804 psDefn->PMLongToGreenwich = 01805 GTIFAngleToDD( psDefn->PMLongToGreenwich, 01806 psDefn->UOMAngle ); 01807 } 01808 01809 /* -------------------------------------------------------------------- */ 01810 /* Have the projection units of measure been overridden? We */ 01811 /* should likely be doing something about angular units too, */ 01812 /* but these are very rarely not decimal degrees for actual */ 01813 /* file coordinates. */ 01814 /* -------------------------------------------------------------------- */ 01815 GTIFKeyGet(psGTIF,ProjLinearUnitsGeoKey,&(psDefn->UOMLength),0,1); 01816 01817 if( psDefn->UOMLength != KvUserDefined ) 01818 { 01819 GTIFGetUOMLengthInfo( psDefn->UOMLength, NULL, 01820 &(psDefn->UOMLengthInMeters) ); 01821 } 01822 01823 /* -------------------------------------------------------------------- */ 01824 /* Handle a variety of user defined transform types. */ 01825 /* -------------------------------------------------------------------- */ 01826 if( GTIFKeyGet(psGTIF,ProjCoordTransGeoKey, 01827 &(psDefn->CTProjection),0,1) == 1) 01828 { 01829 GTIFFetchProjParms( psGTIF, psDefn ); 01830 } 01831 01832 /* -------------------------------------------------------------------- */ 01833 /* Try to set the zoned map system information. */ 01834 /* -------------------------------------------------------------------- */ 01835 psDefn->MapSys = GTIFProjToMapSys( psDefn->ProjCode, &(psDefn->Zone) ); 01836 01837 /* -------------------------------------------------------------------- */ 01838 /* If this is UTM, and we were unable to extract the projection */ 01839 /* parameters from the CSV file, just set them directly now, */ 01840 /* since it's pretty easy, and a common case. */ 01841 /* -------------------------------------------------------------------- */ 01842 if( (psDefn->MapSys == MapSys_UTM_North 01843 || psDefn->MapSys == MapSys_UTM_South) 01844 && psDefn->CTProjection == KvUserDefined ) 01845 { 01846 psDefn->CTProjection = CT_TransverseMercator; 01847 psDefn->nParms = 7; 01848 psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey; 01849 psDefn->ProjParm[0] = 0.0; 01850 01851 psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey; 01852 psDefn->ProjParm[1] = psDefn->Zone*6 - 183.0; 01853 01854 psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey; 01855 psDefn->ProjParm[4] = 0.9996; 01856 01857 psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; 01858 psDefn->ProjParm[5] = 500000.0; 01859 01860 psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; 01861 01862 if( psDefn->MapSys == MapSys_UTM_North ) 01863 psDefn->ProjParm[6] = 0.0; 01864 else 01865 psDefn->ProjParm[6] = 10000000.0; 01866 } 01867 01868 /* -------------------------------------------------------------------- */ 01869 /* For now we forceable deaccess all CSV files to reduce the */ 01870 /* chance of "leakage". Really, this should be application */ 01871 /* controlled. */ 01872 /* -------------------------------------------------------------------- */ 01873 CSVDeaccess( NULL ); 01874 01875 return TRUE; 01876 } 01877 01878 /************************************************************************/ 01879 /* GTIFDecToDMS() */ 01880 /* */ 01881 /* Convenient function to translate decimal degrees to DMS */ 01882 /* format for reporting to a user. */ 01883 /************************************************************************/ 01884 01885 const char *GTIFDecToDMS( double dfAngle, const char * pszAxis, 01886 int nPrecision ) 01887 01888 { 01889 int nDegrees, nMinutes; 01890 double dfSeconds; 01891 char szFormat[30]; 01892 static char szBuffer[50]; 01893 const char *pszHemisphere = NULL; 01894 double dfRound; 01895 int i; 01896 01897 dfRound = 0.5/60; 01898 for( i = 0; i < nPrecision; i++ ) 01899 dfRound = dfRound * 0.1; 01900 01901 nDegrees = (int) ABS(dfAngle); 01902 nMinutes = (int) ((ABS(dfAngle) - nDegrees) * 60 + dfRound); 01903 dfSeconds = ABS((ABS(dfAngle) * 3600 - nDegrees*3600 - nMinutes*60)); 01904 01905 if( EQUAL(pszAxis,"Long") && dfAngle < 0.0 ) 01906 pszHemisphere = "W"; 01907 else if( EQUAL(pszAxis,"Long") ) 01908 pszHemisphere = "E"; 01909 else if( dfAngle < 0.0 ) 01910 pszHemisphere = "S"; 01911 else 01912 pszHemisphere = "N"; 01913 01914 sprintf( szFormat, "%%3dd%%2d\'%%%d.%df\"%s", 01915 nPrecision+3, nPrecision, pszHemisphere ); 01916 sprintf( szBuffer, szFormat, nDegrees, nMinutes, dfSeconds ); 01917 01918 return( szBuffer ); 01919 } 01920 01921 /************************************************************************/ 01922 /* GTIFPrintDefn() */ 01923 /* */ 01924 /* Report the contents of a GTIFDefn structure ... mostly for */ 01925 /* debugging. */ 01926 /************************************************************************/ 01927 01928 void GTIFPrintDefn( GTIFDefn * psDefn, FILE * fp ) 01929 01930 { 01931 /* -------------------------------------------------------------------- */ 01932 /* Get the PCS name if possible. */ 01933 /* -------------------------------------------------------------------- */ 01934 if( psDefn->PCS != KvUserDefined ) 01935 { 01936 char *pszPCSName = CPLStrdup("name unknown"); 01937 01938 GTIFGetPCSInfo( psDefn->PCS, &pszPCSName, NULL, NULL, NULL, NULL ); 01939 fprintf( fp, "PCS = %d (%s)\n", psDefn->PCS, pszPCSName ); 01940 CPLFree( pszPCSName ); 01941 } 01942 01943 /* -------------------------------------------------------------------- */ 01944 /* Dump the projection code if possible. */ 01945 /* -------------------------------------------------------------------- */ 01946 if( psDefn->ProjCode != KvUserDefined ) 01947 { 01948 char *pszTRFName = CPLStrdup(""); 01949 01950 GTIFGetProjTRFInfo( psDefn->ProjCode, &pszTRFName, NULL, NULL ); 01951 fprintf( fp, "Projection = %d (%s)\n", 01952 psDefn->ProjCode, pszTRFName ); 01953 01954 CPLFree( pszTRFName ); 01955 } 01956 01957 /* -------------------------------------------------------------------- */ 01958 /* Try to dump the projection method name, and parameters if possible.*/ 01959 /* -------------------------------------------------------------------- */ 01960 if( psDefn->CTProjection != KvUserDefined ) 01961 { 01962 char *pszName = GTIFValueName(ProjCoordTransGeoKey, 01963 psDefn->CTProjection); 01964 int i; 01965 01966 if( pszName == NULL ) 01967 pszName = "(unknown)"; 01968 01969 fprintf( fp, "Projection Method: %s\n", pszName ); 01970 01971 for( i = 0; i < psDefn->nParms; i++ ) 01972 { 01973 if( psDefn->ProjParmId[i] == 0 ) 01974 continue; 01975 01976 pszName = GTIFKeyName((geokey_t) psDefn->ProjParmId[i]); 01977 if( pszName == NULL ) 01978 pszName = "(unknown)"; 01979 01980 if( i < 4 ) 01981 { 01982 char *pszAxisName; 01983 01984 if( strstr(pszName,"Long") != NULL ) 01985 pszAxisName = "Long"; 01986 else if( strstr(pszName,"Lat") != NULL ) 01987 pszAxisName = "Lat"; 01988 else 01989 pszAxisName = "?"; 01990 01991 fprintf( fp, " %s: %f (%s)\n", 01992 pszName, psDefn->ProjParm[i], 01993 GTIFDecToDMS( psDefn->ProjParm[i], pszAxisName, 2 ) ); 01994 } 01995 else if( i == 4 ) 01996 fprintf( fp, " %s: %f\n", pszName, psDefn->ProjParm[i] ); 01997 else 01998 fprintf( fp, " %s: %f m\n", pszName, psDefn->ProjParm[i] ); 01999 } 02000 } 02001 02002 /* -------------------------------------------------------------------- */ 02003 /* Report the GCS name, and number. */ 02004 /* -------------------------------------------------------------------- */ 02005 if( psDefn->GCS != KvUserDefined ) 02006 { 02007 char *pszName = CPLStrdup("(unknown)"); 02008 02009 GTIFGetGCSInfo( psDefn->GCS, &pszName, NULL, NULL, NULL ); 02010 fprintf( fp, "GCS: %d/%s\n", psDefn->GCS, pszName ); 02011 CPLFree( pszName ); 02012 } 02013 02014 /* -------------------------------------------------------------------- */ 02015 /* Report the datum name. */ 02016 /* -------------------------------------------------------------------- */ 02017 if( psDefn->Datum != KvUserDefined ) 02018 { 02019 char *pszName = CPLStrdup("(unknown)"); 02020 02021 GTIFGetDatumInfo( psDefn->Datum, &pszName, NULL ); 02022 fprintf( fp, "Datum: %d/%s\n", psDefn->Datum, pszName ); 02023 CPLFree( pszName ); 02024 } 02025 02026 /* -------------------------------------------------------------------- */ 02027 /* Report the ellipsoid. */ 02028 /* -------------------------------------------------------------------- */ 02029 if( psDefn->Ellipsoid != KvUserDefined ) 02030 { 02031 char *pszName = CPLStrdup("(unknown)"); 02032 02033 GTIFGetEllipsoidInfo( psDefn->Ellipsoid, &pszName, NULL, NULL ); 02034 fprintf( fp, "Ellipsoid: %d/%s (%.2f,%.2f)\n", 02035 psDefn->Ellipsoid, pszName, 02036 psDefn->SemiMajor, psDefn->SemiMinor ); 02037 CPLFree( pszName ); 02038 } 02039 02040 /* -------------------------------------------------------------------- */ 02041 /* Report the prime meridian. */ 02042 /* -------------------------------------------------------------------- */ 02043 if( psDefn->PM != KvUserDefined ) 02044 { 02045 char *pszName = CPLStrdup("(unknown)"); 02046 02047 GTIFGetPMInfo( psDefn->PM, &pszName, NULL ); 02048 fprintf( fp, "Prime Meridian: %d/%s (%f/%s)\n", 02049 psDefn->PM, pszName, 02050 psDefn->PMLongToGreenwich, 02051 GTIFDecToDMS( psDefn->PMLongToGreenwich, "Long", 2 ) ); 02052 CPLFree( pszName ); 02053 } 02054 02055 /* -------------------------------------------------------------------- */ 02056 /* Report the projection units of measure (currently just */ 02057 /* linear). */ 02058 /* -------------------------------------------------------------------- */ 02059 if( psDefn->UOMLength != KvUserDefined ) 02060 { 02061 char *pszName = CPLStrdup( "(unknown)" ); 02062 02063 GTIFGetUOMLengthInfo( psDefn->UOMLength, &pszName, NULL ); 02064 fprintf( fp, "Projection Linear Units: %d/%s (%fm)\n", 02065 psDefn->UOMLength, pszName, psDefn->UOMLengthInMeters ); 02066 CPLFree( pszName ); 02067 } 02068 } 02069
1.2.3-20001105 written by Dimitri van Heesch,
© 1997-2000