1 /******************************************************************************
2 * $Id: ogr_srs_pci.cpp 28565 2015-02-27 10:26:21Z rouault $
3 *
4 * Project: OpenGIS Simple Features Reference Implementation
5 * Purpose: OGRSpatialReference translation to/from PCI georeferencing
6 * information.
7 * Author: Andrey Kiselev, dron@ak4719.spb.edu
8 *
9 ******************************************************************************
10 * Copyright (c) 2003, Andrey Kiselev <dron@ak4719.spb.edu>
11 * Copyright (c) 2009-2011, Even Rouault <even dot rouault at mines-paris dot org>
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining a
14 * copy of this software and associated documentation files (the "Software"),
15 * to deal in the Software without restriction, including without limitation
16 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 * and/or sell copies of the Software, and to permit persons to whom the
18 * Software is furnished to do so, subject to the following conditions:
19 *
20 * The above copyright notice and this permission notice shall be included
21 * in all copies or substantial portions of the Software.
22 *
23 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 * DEALINGS IN THE SOFTWARE.
30 ****************************************************************************/
31
32 #include "ogr_spatialref.h"
33 #include "ogr_p.h"
34 #include "cpl_conv.h"
35 #include "cpl_csv.h"
36
37 CPL_CVSID("$Id: ogr_srs_pci.cpp 28565 2015-02-27 10:26:21Z rouault $");
38
39 typedef struct
40 {
41 const char *pszPCIDatum;
42 int nEPSGCode;
43 } PCIDatums;
44
45 static const PCIDatums asDatums[] =
46 {
47 { "D-01", 4267 }, // NAD27 (USA, NADCON)
48 { "D-03", 4267 }, // NAD27 (Canada, NTv1)
49 { "D-02", 4269 }, // NAD83 (USA, NADCON)
50 { "D-04", 4269 }, // NAD83 (Canada, NTv1)
51 { "D000", 4326 }, // WGS 1984
52 { "D001", 4322 }, // WGS 1972
53 { "D008", 4296 }, // Sudan
54 { "D013", 4601 }, // Antigua Island Astro 1943
55 { "D029", 4202 }, // Australian Geodetic 1966
56 { "D030", 4203 }, // Australian Geodetic 1984
57 { "D033", 4216 }, // Bermuda 1957
58 { "D034", 4165 }, // Bissau
59 { "D036", 4219 }, // Bukit Rimpah
60 { "D038", 4221 }, // Campo Inchauspe
61 { "D040", 4222 }, // Cape
62 { "D042", 4223 }, // Carthage
63 { "D044", 4224 }, // Chua Astro
64 { "D045", 4225 }, // Corrego Alegre
65 { "D046", 4155 }, // Dabola (Guinea)
66 { "D066", 4272 }, // Geodetic Datum 1949 (New Zealand)
67 { "D071", 4255 }, // Herat North (Afghanistan)
68 { "D077", 4239 }, // Indian 1954 (Thailand, Vietnam)
69 { "D078", 4240 }, // Indian 1975 (Thailand)
70 { "D083", 4244 }, // Kandawala (Sri Lanka)
71 { "D085", 4245 }, // Kertau 1948 (West Malaysia & Singapore)
72 { "D088", 4250 }, // Leigon (Ghana)
73 { "D089", 4251 }, // Liberia 1964 (Liberia)
74 { "D092", 4256 }, // Mahe 1971 (Mahe Island)
75 { "D093", 4262 }, // Massawa (Ethiopia (Eritrea))
76 { "D094", 4261 }, // Merchich (Morocco)
77 { "D098", 4604 }, // Montserrat Island Astro 1958 (Montserrat (Leeward Islands))
78 { "D110", 4267 }, // NAD27 / Alaska
79 { "D139", 4282 }, // Pointe Noire 1948 (Congo)
80 { "D140", 4615 }, // Porto Santo 1936 (Porto Santo, Madeira Islands)
81 { "D151", 4139 }, // Puerto Rico (Puerto Rico, Virgin Islands)
82 { "D153", 4287 }, // Qornoq (Greenland (South))
83 { "D158", 4292 }, // Sapper Hill 1943 (East Falkland Island)
84 { "D159", 4293 }, // Schwarzeck (Namibia)
85 { "D160", 4616 }, // Selvagem Grande 1938 (Salvage Islands)
86 { "D176", 4297 }, // Tananarive Observatory 1925 (Madagascar)
87 { "D177", 4298 }, // Timbalai 1948 (Brunei, East Malaysia (Sabah, Sarawak))
88 { "D187", 4309 }, // Yacare (Uruguay)
89 { "D188", 4311 }, // Zanderij (Suriname)
90 { "D401", 4124 }, // RT90 (Sweden)
91 { "D501", 4312 }, // MGI (Hermannskogel, Austria)
92 { NULL, 0 }
93 };
94
95 static const PCIDatums asEllips[] =
96 {
97 { "E000", 7008 }, // Clarke, 1866 (NAD1927)
98 { "E001", 7034 }, // Clarke, 1880
99 { "E002", 7004 }, // Bessel, 1841
100 { "E004", 7022 }, // International, 1924 (Hayford, 1909)
101 { "E005", 7043 }, // WGS, 1972
102 { "E006", 7042 }, // Everest, 1830
103 { "E008", 7019 }, // GRS, 1980 (NAD1983)
104 { "E009", 7001 }, // Airy, 1830
105 { "E010", 7018 }, // Modified Everest
106 { "E011", 7002 }, // Modified Airy
107 { "E012", 7030 }, // WGS, 1984 (GPS)
108 { "E014", 7003 }, // Australian National, 1965
109 { "E015", 7024 }, // Krassovsky, 1940
110 { "E016", 7053 }, // Hough
111 { "E019", 7052 }, // normal sphere
112 { "E333", 7046 }, // Bessel 1841 (Japan By Law)
113 { "E900", 7006 }, // Bessel, 1841 (Namibia)
114 { "E901", 7044 }, // Everest, 1956
115 { "E902", 7056 }, // Everest, 1969
116 { "E903", 7016 }, // Everest (Sabah & Sarawak)
117 { "E904", 7020 }, // Helmert, 1906
118 { "E907", 7036 }, // South American, 1969
119 { "E910", 7041 }, // ATS77
120 { NULL, 0 }
121 };
122
123 /************************************************************************/
124 /* OSRImportFromPCI() */
125 /************************************************************************/
126
127 /**
128 * \brief Import coordinate system from PCI projection definition.
129 *
130 * This function is the same as OGRSpatialReference::importFromPCI().
131 */
132
OSRImportFromPCI(OGRSpatialReferenceH hSRS,const char * pszProj,const char * pszUnits,double * padfPrjParams)133 OGRErr OSRImportFromPCI( OGRSpatialReferenceH hSRS, const char *pszProj,
134 const char *pszUnits, double *padfPrjParams )
135
136 {
137 VALIDATE_POINTER1( hSRS, "OSRImportFromPCI", CE_Failure );
138
139 return ((OGRSpatialReference *) hSRS)->importFromPCI( pszProj,
140 pszUnits,
141 padfPrjParams );
142 }
143
144 /************************************************************************/
145 /* importFromPCI() */
146 /************************************************************************/
147
148 /**
149 * \brief Import coordinate system from PCI projection definition.
150 *
151 * PCI software uses 16-character string to specify coordinate system
152 * and datum/ellipsoid. You should supply at least this string to the
153 * importFromPCI() function.
154 *
155 * This function is the equivalent of the C function OSRImportFromPCI().
156 *
157 * @param pszProj NULL terminated string containing the definition. Looks
158 * like "pppppppppppp Ennn" or "pppppppppppp Dnnn", where "pppppppppppp" is
159 * a projection code, "Ennn" is an ellipsoid code, "Dnnn" --- a datum code.
160 *
161 * @param pszUnits Grid units code ("DEGREE" or "METRE"). If NULL "METRE" will
162 * be used.
163 *
164 * @param padfPrjParams Array of 17 coordinate system parameters:
165 *
166 * [0] Spheroid semi major axis
167 * [1] Spheroid semi minor axis
168 * [2] Reference Longitude
169 * [3] Reference Latitude
170 * [4] First Standard Parallel
171 * [5] Second Standard Parallel
172 * [6] False Easting
173 * [7] False Northing
174 * [8] Scale Factor
175 * [9] Height above sphere surface
176 * [10] Longitude of 1st point on center line
177 * [11] Latitude of 1st point on center line
178 * [12] Longitude of 2nd point on center line
179 * [13] Latitude of 2nd point on center line
180 * [14] Azimuth east of north for center line
181 * [15] Landsat satellite number
182 * [16] Landsat path number
183 *
184 * Particular projection uses different parameters, unused ones may be set to
185 * zero. If NULL suppliet instead of array pointer default values will be
186 * used (i.e., zeroes).
187 *
188 * @return OGRERR_NONE on success or an error code in case of failure.
189 */
190
importFromPCI(const char * pszProj,const char * pszUnits,double * padfPrjParams)191 OGRErr OGRSpatialReference::importFromPCI( const char *pszProj,
192 const char *pszUnits,
193 double *padfPrjParams )
194
195 {
196 Clear();
197
198 if( pszProj == NULL || CPLStrnlen(pszProj, 16) < 16 )
199 return OGRERR_CORRUPT_DATA;
200
201 CPLDebug( "OSR_PCI", "Trying to import projection \"%s\"", pszProj );
202
203 /* -------------------------------------------------------------------- */
204 /* Use safe defaults if projection parameters are not supplied. */
205 /* -------------------------------------------------------------------- */
206 int bProjAllocated = FALSE;
207
208 if( padfPrjParams == NULL )
209 {
210 int i;
211
212 padfPrjParams = (double *)CPLMalloc( 17 * sizeof(double) );
213 if ( !padfPrjParams )
214 return OGRERR_NOT_ENOUGH_MEMORY;
215 for ( i = 0; i < 17; i++ )
216 padfPrjParams[i] = 0.0;
217 bProjAllocated = TRUE;
218 }
219
220 /* -------------------------------------------------------------------- */
221 /* Extract and "normalize" the earthmodel to look like E001, */
222 /* D-02 or D109. */
223 /* -------------------------------------------------------------------- */
224 char szEarthModel[5];
225 const char *pszEM;
226 int bIsNAD27 = FALSE;
227
228 strcpy( szEarthModel, "" );
229 pszEM = pszProj + strlen(pszProj) - 1;
230 while( pszEM != pszProj )
231 {
232 if( *pszEM == 'e' || *pszEM == 'E' || *pszEM == 'd' || *pszEM == 'D' )
233 {
234 int nCode = atoi(pszEM+1);
235
236 if( nCode >= -99 && nCode <= 999 )
237 sprintf( szEarthModel, "%c%03d", toupper(*pszEM), nCode );
238
239 break;
240 }
241
242 pszEM--;
243 }
244
245 if( EQUAL(pszEM,"E000")
246 || EQUAL(pszEM,"D-01")
247 || EQUAL(pszEM,"D-03")
248 || EQUAL(pszEM,"D-07")
249 || EQUAL(pszEM,"D-09")
250 || EQUAL(pszEM,"D-11")
251 || EQUAL(pszEM,"D-13")
252 || EQUAL(pszEM,"D-17") )
253 bIsNAD27 = TRUE;
254
255 /* -------------------------------------------------------------------- */
256 /* Operate on the basis of the projection name. */
257 /* -------------------------------------------------------------------- */
258 if( EQUALN( pszProj, "LONG/LAT", 8 ) )
259 {
260 }
261
262 else if( EQUALN( pszProj, "METER", 5 )
263 || EQUALN( pszProj, "METRE", 5 ) )
264 {
265 SetLocalCS( "METER" );
266 SetLinearUnits( "METER", 1.0 );
267 }
268
269 else if( EQUALN( pszProj, "FEET", 4 )
270 || EQUALN( pszProj, "FOOT", 4 ) )
271 {
272 SetLocalCS( "FEET" );
273 SetLinearUnits( "FEET", CPLAtof(SRS_UL_FOOT_CONV) );
274 }
275
276 else if( EQUALN( pszProj, "ACEA", 4 ) )
277 {
278 SetACEA( padfPrjParams[4], padfPrjParams[5],
279 padfPrjParams[3], padfPrjParams[2],
280 padfPrjParams[6], padfPrjParams[7] );
281 }
282
283 else if( EQUALN( pszProj, "AE", 2 ) )
284 {
285 SetAE( padfPrjParams[3], padfPrjParams[2],
286 padfPrjParams[6], padfPrjParams[7] );
287 }
288
289 else if( EQUALN( pszProj, "CASS ", 5 ) )
290 {
291 SetCS( padfPrjParams[3], padfPrjParams[2],
292 padfPrjParams[6], padfPrjParams[7] );
293 }
294
295 else if( EQUALN( pszProj, "EC", 2 ) )
296 {
297 SetEC( padfPrjParams[4], padfPrjParams[5],
298 padfPrjParams[3], padfPrjParams[2],
299 padfPrjParams[6], padfPrjParams[7] );
300 }
301
302 else if( EQUALN( pszProj, "ER", 2 ) )
303 {
304 // PCI and GCTP don't support natural origin lat.
305 SetEquirectangular2( 0.0, padfPrjParams[2],
306 padfPrjParams[3],
307 padfPrjParams[6], padfPrjParams[7] );
308 }
309
310 else if( EQUALN( pszProj, "GNO", 3 ) )
311 {
312 SetGnomonic( padfPrjParams[3], padfPrjParams[2],
313 padfPrjParams[6], padfPrjParams[7] );
314 }
315
316 // FIXME: GVNP --- General Vertical Near- Side Perspective skipped
317
318 // FIXME: GOOD -- our Goode's is not the interrupted version from pci
319
320 else if( EQUALN( pszProj, "LAEA", 4 ) )
321 {
322 SetLAEA( padfPrjParams[3], padfPrjParams[2],
323 padfPrjParams[6], padfPrjParams[7] );
324 }
325
326 else if( EQUALN( pszProj, "LCC ", 4 ) )
327 {
328 SetLCC( padfPrjParams[4], padfPrjParams[5],
329 padfPrjParams[3], padfPrjParams[2],
330 padfPrjParams[6], padfPrjParams[7] );
331 }
332
333 else if( EQUALN( pszProj, "LCC_1SP ", 7 ) )
334 {
335 SetLCC1SP( padfPrjParams[3], padfPrjParams[2],
336 padfPrjParams[8],
337 padfPrjParams[6], padfPrjParams[7] );
338 }
339
340 else if( EQUALN( pszProj, "MC", 2 ) )
341 {
342 SetMC( padfPrjParams[3], padfPrjParams[2],
343 padfPrjParams[6], padfPrjParams[7] );
344 }
345
346 else if( EQUALN( pszProj, "MER", 3 ) )
347 {
348 SetMercator( padfPrjParams[3], padfPrjParams[2],
349 (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
350 padfPrjParams[6], padfPrjParams[7] );
351 }
352
353 else if( EQUALN( pszProj, "OG", 2 ) )
354 {
355 SetOrthographic( padfPrjParams[3], padfPrjParams[2],
356 padfPrjParams[6], padfPrjParams[7] );
357 }
358
359 else if( EQUALN( pszProj, "OM ", 3 ) )
360 {
361 if( padfPrjParams[10] == 0.0
362 && padfPrjParams[11] == 0.0
363 && padfPrjParams[12] == 0.0
364 && padfPrjParams[13] == 0.0 )
365 {
366 SetHOM( padfPrjParams[3], padfPrjParams[2],
367 padfPrjParams[14],
368 padfPrjParams[14], // use azimuth for grid angle
369 padfPrjParams[8],
370 padfPrjParams[6], padfPrjParams[7] );
371 }
372 else
373 {
374 SetHOM2PNO( padfPrjParams[3],
375 padfPrjParams[11], padfPrjParams[10],
376 padfPrjParams[13], padfPrjParams[12],
377 padfPrjParams[8],
378 padfPrjParams[6], padfPrjParams[7] );
379 }
380 }
381
382 else if( EQUALN( pszProj, "PC", 2 ) )
383 {
384 SetPolyconic( padfPrjParams[3], padfPrjParams[2],
385 padfPrjParams[6], padfPrjParams[7] );
386 }
387
388 else if( EQUALN( pszProj, "PS", 2 ) )
389 {
390 SetPS( padfPrjParams[3], padfPrjParams[2],
391 (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
392 padfPrjParams[6], padfPrjParams[7] );
393 }
394
395 else if( EQUALN( pszProj, "ROB", 3 ) )
396 {
397 SetRobinson( padfPrjParams[2],
398 padfPrjParams[6], padfPrjParams[7] );
399 }
400
401 else if( EQUALN( pszProj, "SGDO", 4 ) )
402 {
403 SetOS( padfPrjParams[3], padfPrjParams[2],
404 (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
405 padfPrjParams[6], padfPrjParams[7] );
406 }
407
408 else if( EQUALN( pszProj, "SG", 2 ) )
409 {
410 SetStereographic( padfPrjParams[3], padfPrjParams[2],
411 (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
412 padfPrjParams[6], padfPrjParams[7] );
413 }
414
415 else if( EQUALN( pszProj, "SIN", 3 ) )
416 {
417 SetSinusoidal( padfPrjParams[2],
418 padfPrjParams[6], padfPrjParams[7] );
419 }
420
421 // FIXME: SOM --- Space Oblique Mercator skipped
422
423 else if( EQUALN( pszProj, "SPCS", 4 ) )
424 {
425 int iZone;
426
427 iZone = CPLScanLong( (char *)pszProj + 5, 4 );
428
429 SetStatePlane( iZone, !bIsNAD27 );
430 SetLinearUnitsAndUpdateParameters( SRS_UL_METER, 1.0 );
431 }
432
433 else if( EQUALN( pszProj, "SPIF", 4 ) )
434 {
435 int iZone;
436
437 iZone = CPLScanLong( (char *)pszProj + 5, 4 );
438
439 SetStatePlane( iZone, !bIsNAD27 );
440 SetLinearUnitsAndUpdateParameters( SRS_UL_FOOT,
441 CPLAtof(SRS_UL_FOOT_CONV) );
442 }
443
444 else if( EQUALN( pszProj, "SPAF", 4 ) )
445 {
446 int iZone;
447
448 iZone = CPLScanLong( (char *)pszProj + 5, 4 );
449
450 SetStatePlane( iZone, !bIsNAD27 );
451 SetLinearUnitsAndUpdateParameters( SRS_UL_US_FOOT,
452 CPLAtof(SRS_UL_US_FOOT_CONV) );
453 }
454
455 else if( EQUALN( pszProj, "TM", 2 ) )
456 {
457 SetTM( padfPrjParams[3], padfPrjParams[2],
458 (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
459 padfPrjParams[6], padfPrjParams[7] );
460 }
461
462 else if( EQUALN( pszProj, "UTM", 3 ) )
463 {
464 int iZone, bNorth = TRUE;
465
466 iZone = CPLScanLong( (char *)pszProj + 4, 5 );;
467 if ( iZone < 0 )
468 {
469 iZone = -iZone;
470 bNorth = FALSE;
471 }
472
473 // Check for a zone letter. PCI uses, accidentally, MGRS
474 // type row lettering in its UTM projection
475 char byZoneID = 0;
476
477 if( strlen(pszProj) > 10 && pszProj[10] != ' ' )
478 byZoneID = pszProj[10];
479
480 // Determine if the MGRS zone falls above or below the equator
481 if (byZoneID != 0 )
482 {
483 CPLDebug("OSR_PCI", "Found MGRS zone in UTM projection string: %c",
484 byZoneID);
485
486 if (byZoneID >= 'N' && byZoneID <= 'X')
487 {
488 bNorth = TRUE;
489 }
490 else if (byZoneID >= 'C' && byZoneID <= 'M')
491 {
492 bNorth = FALSE;
493 }
494 else
495 {
496 // yikes, most likely we got something that was not really
497 // an MGRS zone code so we ignore it.
498 }
499 }
500
501 SetUTM( iZone, bNorth );
502 }
503
504 else if( EQUALN( pszProj, "VDG", 3 ) )
505 {
506 SetVDG( padfPrjParams[2],
507 padfPrjParams[6], padfPrjParams[7] );
508 }
509
510 else
511 {
512 CPLDebug( "OSR_PCI", "Unsupported projection: %s", pszProj );
513 SetLocalCS( pszProj );
514 }
515
516 /* ==================================================================== */
517 /* Translate the datum/spheroid. */
518 /* ==================================================================== */
519
520 /* -------------------------------------------------------------------- */
521 /* We have an earthmodel string, look it up in the datum list. */
522 /* -------------------------------------------------------------------- */
523 if( strlen(szEarthModel) > 0
524 && (poRoot == NULL || IsProjected() || IsGeographic()) )
525 {
526 const PCIDatums *pasDatum = asDatums;
527
528 // Search for matching datum
529 while ( pasDatum->pszPCIDatum )
530 {
531 if( EQUALN( szEarthModel, pasDatum->pszPCIDatum, 4 ) )
532 {
533 OGRSpatialReference oGCS;
534 oGCS.importFromEPSG( pasDatum->nEPSGCode );
535 CopyGeogCSFrom( &oGCS );
536 break;
537 }
538 pasDatum++;
539 }
540
541 /* -------------------------------------------------------------------- */
542 /* If we did not find a datum definition in our incode epsg */
543 /* lookup table, then try fetching from the pci_datum.txt */
544 /* file. */
545 /* -------------------------------------------------------------------- */
546 char **papszDatumDefn = NULL;
547
548 if( !pasDatum->pszPCIDatum && szEarthModel[0] == 'D' )
549 {
550 const char *pszDatumCSV = CSVFilename( "pci_datum.txt" );
551 FILE *fp = NULL;
552
553 if( pszDatumCSV )
554 fp = VSIFOpen( pszDatumCSV, "r" );
555
556 if( fp != NULL )
557 {
558 char **papszLineItems = NULL;
559
560 while( (papszLineItems = CSVReadParseLine( fp )) != NULL )
561 {
562 if( CSLCount(papszLineItems) > 3
563 && EQUALN(papszLineItems[0],szEarthModel,4) )
564 {
565 papszDatumDefn = papszLineItems;
566 strncpy( szEarthModel, papszLineItems[2], 4 );
567 break;
568 }
569 CSLDestroy( papszLineItems );
570 }
571
572 VSIFClose( fp );
573 }
574 }
575
576 /* -------------------------------------------------------------------- */
577 /* If not, look in the ellipsoid/EPSG matching list. */
578 /* -------------------------------------------------------------------- */
579 if ( !pasDatum->pszPCIDatum ) // No matching; search for ellipsoids
580 {
581 char *pszName = NULL;
582 double dfSemiMajor = 0.0;
583 double dfInvFlattening = 0.0;
584 int nEPSGCode = 0;
585
586 pasDatum = asEllips;
587
588 while ( pasDatum->pszPCIDatum )
589 {
590 if( EQUALN( szEarthModel, pasDatum->pszPCIDatum, 4 ) )
591 {
592 nEPSGCode = pasDatum->nEPSGCode;
593 OSRGetEllipsoidInfo( pasDatum->nEPSGCode, &pszName,
594 &dfSemiMajor, &dfInvFlattening );
595 break;
596
597 }
598 pasDatum++;
599 }
600
601 /* -------------------------------------------------------------------- */
602 /* If we don't find it in that list, do a lookup in the */
603 /* pci_ellips.txt file. */
604 /* -------------------------------------------------------------------- */
605 if( !pasDatum->pszPCIDatum && szEarthModel[0] == 'E' )
606 {
607 const char *pszCSV = CSVFilename( "pci_ellips.txt" );
608 FILE *fp = NULL;
609
610 if( pszCSV )
611 fp = VSIFOpen( pszCSV, "r" );
612
613 if( fp != NULL )
614 {
615 char **papszLineItems = NULL;
616
617 while( (papszLineItems = CSVReadParseLine( fp )) != NULL )
618 {
619 if( CSLCount(papszLineItems) > 3
620 && EQUALN(papszLineItems[0],szEarthModel,4) )
621 {
622 dfSemiMajor = CPLAtof( papszLineItems[2] );
623 double dfSemiMinor = CPLAtof( papszLineItems[3] );
624 dfInvFlattening = OSRCalcInvFlattening(dfSemiMajor, dfSemiMinor);
625 break;
626 }
627 CSLDestroy( papszLineItems );
628 }
629 CSLDestroy( papszLineItems );
630
631 VSIFClose( fp );
632 }
633 }
634
635 /* -------------------------------------------------------------------- */
636 /* Custom spheroid? */
637 /* -------------------------------------------------------------------- */
638 if( dfSemiMajor == 0.0 && EQUALN(szEarthModel,"E999",4)
639 && padfPrjParams[0] != 0.0 )
640 {
641 dfSemiMajor = padfPrjParams[0];
642 double dfSemiMinor = padfPrjParams[1];
643 dfInvFlattening = OSRCalcInvFlattening(dfSemiMajor, dfSemiMinor);
644 }
645
646 /* -------------------------------------------------------------------- */
647 /* If nothing else, fall back to WGS84 parameters. */
648 /* -------------------------------------------------------------------- */
649 if( dfSemiMajor == 0.0 )
650 {
651 dfSemiMajor = SRS_WGS84_SEMIMAJOR;
652 dfInvFlattening = SRS_WGS84_INVFLATTENING;
653 }
654
655 /* -------------------------------------------------------------------- */
656 /* Now try to put this all together into a GEOGCS definition. */
657 /* -------------------------------------------------------------------- */
658 CPLString osGCSName, osDatumName, osEllipseName;
659
660 if( pszName )
661 osEllipseName = pszName;
662 else
663 osEllipseName.Printf( "Unknown - PCI %s", szEarthModel );
664 CPLFree( pszName );
665
666 if( papszDatumDefn )
667 osDatumName = papszDatumDefn[1];
668 else
669 osDatumName.Printf( "Unknown - PCI %s", szEarthModel );
670 osGCSName = osDatumName;
671
672 SetGeogCS( osGCSName, osDatumName, osEllipseName,
673 dfSemiMajor, dfInvFlattening );
674
675 // Do we have an ellipsoid EPSG code?
676 if( nEPSGCode != 0 )
677 SetAuthority( "SPHEROID", "EPSG", nEPSGCode );
678
679 // Do we have 7 datum shift parameters?
680 if( CSLCount(papszDatumDefn) >= 15
681 && CPLAtof(papszDatumDefn[14]) != 0.0 )
682 {
683 double dfScale = CPLAtof(papszDatumDefn[14]);
684
685 // we want scale in parts per million off 1.0
686 // but pci uses a mix of forms.
687 if( dfScale >= 0.999 && dfScale <= 1.001 )
688 dfScale = (dfScale-1.0) * 1000000.0;
689
690 SetTOWGS84( CPLAtof(papszDatumDefn[3]),
691 CPLAtof(papszDatumDefn[4]),
692 CPLAtof(papszDatumDefn[5]),
693 CPLAtof(papszDatumDefn[11]),
694 CPLAtof(papszDatumDefn[12]),
695 CPLAtof(papszDatumDefn[13]),
696 dfScale );
697 }
698
699 // Do we have 7 datum shift parameters?
700 else if( CSLCount(papszDatumDefn) == 11
701 && (CPLAtof(papszDatumDefn[3]) != 0.0
702 || CPLAtof(papszDatumDefn[4]) != 0.0
703 || CPLAtof(papszDatumDefn[5]) != 0.0 ) )
704 {
705 SetTOWGS84( CPLAtof(papszDatumDefn[3]),
706 CPLAtof(papszDatumDefn[4]),
707 CPLAtof(papszDatumDefn[5]) );
708 }
709 }
710
711 CSLDestroy(papszDatumDefn);
712 }
713
714 /* -------------------------------------------------------------------- */
715 /* Grid units translation */
716 /* -------------------------------------------------------------------- */
717 if( (IsLocal() || IsProjected()) && pszUnits )
718 {
719 if( EQUAL( pszUnits, "METRE" ) )
720 SetLinearUnits( SRS_UL_METER, 1.0 );
721 else if( EQUAL( pszUnits, "DEGREE" ) )
722 SetAngularUnits( SRS_UA_DEGREE, CPLAtof(SRS_UA_DEGREE_CONV) );
723 else
724 SetLinearUnits( SRS_UL_METER, 1.0 );
725 }
726
727 FixupOrdering();
728
729 if ( bProjAllocated && padfPrjParams )
730 CPLFree( padfPrjParams );
731
732 return OGRERR_NONE;
733 }
734
735 /************************************************************************/
736 /* OSRExportToPCI() */
737 /************************************************************************/
738 /**
739 * \brief Export coordinate system in PCI projection definition.
740 *
741 * This function is the same as OGRSpatialReference::exportToPCI().
742 */
OSRExportToPCI(OGRSpatialReferenceH hSRS,char ** ppszProj,char ** ppszUnits,double ** ppadfPrjParams)743 OGRErr OSRExportToPCI( OGRSpatialReferenceH hSRS,
744 char **ppszProj, char **ppszUnits,
745 double **ppadfPrjParams )
746
747 {
748 VALIDATE_POINTER1( hSRS, "OSRExportToPCI", CE_Failure );
749
750 *ppszProj = NULL;
751 *ppszUnits = NULL;
752 *ppadfPrjParams = NULL;
753
754 return ((OGRSpatialReference *) hSRS)->exportToPCI( ppszProj, ppszUnits,
755 ppadfPrjParams );
756 }
757
758 /************************************************************************/
759 /* exportToPCI() */
760 /************************************************************************/
761
762 /**
763 * \brief Export coordinate system in PCI projection definition.
764 *
765 * Converts the loaded coordinate reference system into PCI projection
766 * definition to the extent possible. The strings returned in ppszProj,
767 * ppszUnits and ppadfPrjParams array should be deallocated by the caller
768 * with CPLFree() when no longer needed.
769 *
770 * LOCAL_CS coordinate systems are not translatable. An empty string
771 * will be returned along with OGRERR_NONE.
772 *
773 * This method is the equivelent of the C function OSRExportToPCI().
774 *
775 * @param ppszProj pointer to which dynamically allocated PCI projection
776 * definition will be assigned.
777 *
778 * @param ppszUnits pointer to which dynamically allocated units definition
779 * will be assigned.
780 *
781 * @param ppadfPrjParams pointer to which dynamically allocated array of
782 * 17 projection parameters will be assigned. See importFromPCI() for the list
783 * of parameters.
784 *
785 * @return OGRERR_NONE on success or an error code on failure.
786 */
787
exportToPCI(char ** ppszProj,char ** ppszUnits,double ** ppadfPrjParams) const788 OGRErr OGRSpatialReference::exportToPCI( char **ppszProj, char **ppszUnits,
789 double **ppadfPrjParams ) const
790
791 {
792 const char *pszProjection = GetAttrValue("PROJECTION");
793
794 /* -------------------------------------------------------------------- */
795 /* Fill all projection parameters with zero. */
796 /* -------------------------------------------------------------------- */
797 int i;
798
799 *ppadfPrjParams = (double *)CPLMalloc( 17 * sizeof(double) );
800 for ( i = 0; i < 17; i++ )
801 (*ppadfPrjParams)[i] = 0.0;
802
803 /* -------------------------------------------------------------------- */
804 /* Get the prime meridian info. */
805 /* -------------------------------------------------------------------- */
806 #if 0
807 const OGR_SRSNode *poPRIMEM = GetAttrNode( "PRIMEM" );
808 double dfFromGreenwich = 0.0;
809
810 if( poPRIMEM != NULL && poPRIMEM->GetChildCount() >= 2
811 && CPLAtof(poPRIMEM->GetChild(1)->GetValue()) != 0.0 )
812 {
813 dfFromGreenwich = CPLAtof(poPRIMEM->GetChild(1)->GetValue());
814 }
815 #endif
816
817 /* ==================================================================== */
818 /* Handle the projection definition. */
819 /* ==================================================================== */
820 char szProj[17];
821
822 memset( szProj, 0, sizeof(szProj) );
823
824 if( IsLocal() )
825 {
826 if( GetLinearUnits() > 0.30479999 && GetLinearUnits() < 0.3048010 )
827 CPLPrintStringFill( szProj, "FEET", 17 );
828 else
829 CPLPrintStringFill( szProj, "METER", 17 );
830 }
831
832 else if( pszProjection == NULL )
833 {
834 CPLPrintStringFill( szProj, "LONG/LAT", 16 );
835 }
836
837 else if( EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
838 {
839 CPLPrintStringFill( szProj, "ACEA", 16 );
840 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
841 (*ppadfPrjParams)[3] =
842 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
843 (*ppadfPrjParams)[4] =
844 GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 );
845 (*ppadfPrjParams)[5] =
846 GetNormProjParm( SRS_PP_STANDARD_PARALLEL_2, 0.0 );
847 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
848 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
849 }
850
851 else if( EQUAL(pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT) )
852 {
853 CPLPrintStringFill( szProj, "AE", 16 );
854 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
855 (*ppadfPrjParams)[3] =
856 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
857 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
858 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
859 }
860
861 else if( EQUAL(pszProjection, SRS_PT_CASSINI_SOLDNER) )
862 {
863 CPLPrintStringFill( szProj, "CASS", 16 );
864 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
865 (*ppadfPrjParams)[3] =
866 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
867 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
868 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
869 }
870
871 else if( EQUAL(pszProjection, SRS_PT_EQUIDISTANT_CONIC) )
872 {
873 CPLPrintStringFill( szProj, "EC", 16 );
874 (*ppadfPrjParams)[2] =
875 GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 );
876 (*ppadfPrjParams)[3] =
877 GetNormProjParm( SRS_PP_LATITUDE_OF_CENTER, 0.0 );
878 (*ppadfPrjParams)[4] =
879 GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 );
880 (*ppadfPrjParams)[5] =
881 GetNormProjParm( SRS_PP_STANDARD_PARALLEL_2, 0.0 );
882 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
883 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
884 }
885
886 else if( EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR) )
887 {
888 CPLPrintStringFill( szProj, "ER", 16 );
889 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
890 (*ppadfPrjParams)[3] =
891 GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 );
892 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
893 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
894 }
895
896 else if( EQUAL(pszProjection, SRS_PT_GNOMONIC) )
897 {
898 CPLPrintStringFill( szProj, "GNO", 16 );
899 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
900 (*ppadfPrjParams)[3] =
901 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
902 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
903 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
904 }
905
906 else if( EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
907 {
908 CPLPrintStringFill( szProj, "LAEA", 16 );
909 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
910 (*ppadfPrjParams)[3] =
911 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
912 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
913 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
914 }
915
916 else if( EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) )
917 {
918 CPLPrintStringFill( szProj, "LCC", 16 );
919 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
920 (*ppadfPrjParams)[3] =
921 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
922 (*ppadfPrjParams)[4] =
923 GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 );
924 (*ppadfPrjParams)[5] =
925 GetNormProjParm( SRS_PP_STANDARD_PARALLEL_2, 0.0 );
926 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
927 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
928 }
929
930 else if( EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP) )
931 {
932 CPLPrintStringFill( szProj, "LCC_1SP", 16 );
933 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
934 (*ppadfPrjParams)[3] =
935 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
936 (*ppadfPrjParams)[8] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
937 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
938 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
939 }
940
941 else if( EQUAL(pszProjection, SRS_PT_MILLER_CYLINDRICAL) )
942 {
943 CPLPrintStringFill( szProj, "MC", 16 );
944 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
945 (*ppadfPrjParams)[3] =
946 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
947 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
948 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
949 }
950
951 else if( EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) )
952 {
953 CPLPrintStringFill( szProj, "MER", 16 );
954 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
955 (*ppadfPrjParams)[3] =
956 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
957 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
958 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
959 (*ppadfPrjParams)[8] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
960 }
961
962 else if( EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC) )
963 {
964 CPLPrintStringFill( szProj, "OG", 16 );
965 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
966 (*ppadfPrjParams)[3] =
967 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
968 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
969 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
970 }
971
972 else if( EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR) )
973 {
974 CPLPrintStringFill( szProj, "OM", 16 );
975 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER,0.0);
976 (*ppadfPrjParams)[3] = GetNormProjParm( SRS_PP_LATITUDE_OF_CENTER, 0.0);
977 (*ppadfPrjParams)[14] = GetNormProjParm( SRS_PP_AZIMUTH, 0.0);
978 // note we are ignoring rectified_grid_angle which has no pci analog.
979 (*ppadfPrjParams)[8] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 0.0);
980 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
981 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
982 }
983
984 else if( EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN) )
985 {
986 CPLPrintStringFill( szProj, "OM", 16 );
987 (*ppadfPrjParams)[3] = GetNormProjParm( SRS_PP_LATITUDE_OF_CENTER, 0.0);
988 (*ppadfPrjParams)[11] = GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_1,0.0);
989 (*ppadfPrjParams)[10] = GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_1,0.0);
990 (*ppadfPrjParams)[13] = GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_2,0.0);
991 (*ppadfPrjParams)[12] = GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_2,0.0);
992 (*ppadfPrjParams)[8] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 0.0);
993 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
994 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
995 }
996
997 else if( EQUAL(pszProjection, SRS_PT_POLYCONIC) )
998 {
999 CPLPrintStringFill( szProj, "PC", 16 );
1000 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
1001 (*ppadfPrjParams)[3] =
1002 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
1003 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1004 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1005 }
1006
1007 else if( EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC) )
1008 {
1009 CPLPrintStringFill( szProj, "PS", 16 );
1010 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
1011 (*ppadfPrjParams)[3] =
1012 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
1013 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1014 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1015 (*ppadfPrjParams)[8] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
1016 }
1017
1018 else if( EQUAL(pszProjection, SRS_PT_ROBINSON) )
1019 {
1020 CPLPrintStringFill( szProj, "ROB", 16 );
1021 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
1022 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1023 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1024 }
1025
1026 else if( EQUAL(pszProjection, SRS_PT_OBLIQUE_STEREOGRAPHIC) )
1027 {
1028 CPLPrintStringFill( szProj, "SGDO", 16 );
1029 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
1030 (*ppadfPrjParams)[3] =
1031 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
1032 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1033 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1034 (*ppadfPrjParams)[8] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
1035 }
1036
1037 else if( EQUAL(pszProjection, SRS_PT_STEREOGRAPHIC) )
1038 {
1039 CPLPrintStringFill( szProj, "SG", 16 );
1040 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
1041 (*ppadfPrjParams)[3] =
1042 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
1043 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1044 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1045 (*ppadfPrjParams)[8] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
1046 }
1047
1048 else if( EQUAL(pszProjection, SRS_PT_SINUSOIDAL) )
1049 {
1050 CPLPrintStringFill( szProj, "SIN", 16 );
1051 (*ppadfPrjParams)[2] =
1052 GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 );
1053 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1054 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1055 }
1056
1057 else if( EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR) )
1058 {
1059 int bNorth;
1060 int nZone = GetUTMZone( &bNorth );
1061
1062 if( nZone != 0 )
1063 {
1064 CPLPrintStringFill( szProj, "UTM", 16 );
1065 if( bNorth )
1066 CPLPrintInt32( szProj + 5, nZone, 4 );
1067 else
1068 CPLPrintInt32( szProj + 5, -nZone, 4 );
1069 }
1070 else
1071 {
1072 CPLPrintStringFill( szProj, "TM", 16 );
1073 (*ppadfPrjParams)[2] =
1074 GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
1075 (*ppadfPrjParams)[3] =
1076 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 );
1077 (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
1078 (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
1079 (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1080 }
1081 }
1082
1083 else if( EQUAL(pszProjection, SRS_PT_VANDERGRINTEN) )
1084 {
1085 CPLPrintStringFill( szProj, "VDG", 16 );
1086 (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 );
1087 (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1088 (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1089 }
1090
1091 // Projection unsupported by PCI
1092 else
1093 {
1094 CPLDebug( "OSR_PCI",
1095 "Projection \"%s\" unsupported by PCI. "
1096 "PIXEL value will be used.", pszProjection );
1097 CPLPrintStringFill( szProj, "PIXEL", 16 );
1098 }
1099
1100 /* ==================================================================== */
1101 /* Translate the earth model. */
1102 /* ==================================================================== */
1103
1104 /* -------------------------------------------------------------------- */
1105 /* Is this a well known datum? */
1106 /* -------------------------------------------------------------------- */
1107 const char *pszDatum = GetAttrValue( "DATUM" );
1108 char szEarthModel[5];
1109
1110 memset( szEarthModel, 0, sizeof(szEarthModel) );
1111
1112 if( pszDatum == NULL || strlen(pszDatum) == 0 )
1113 /* do nothing */;
1114 else if( EQUAL( pszDatum, SRS_DN_NAD27 ) )
1115 CPLPrintStringFill( szEarthModel, "D-01", 4 );
1116
1117 else if( EQUAL( pszDatum, SRS_DN_NAD83 ) )
1118 CPLPrintStringFill( szEarthModel, "D-02", 4 );
1119
1120 else if( EQUAL( pszDatum, SRS_DN_WGS84 ) )
1121 CPLPrintStringFill( szEarthModel, "D000", 4 );
1122
1123 /* -------------------------------------------------------------------- */
1124 /* If not a very well known datum, try for an EPSG based */
1125 /* translation. */
1126 /* -------------------------------------------------------------------- */
1127 if( szEarthModel[0] == '\0' )
1128 {
1129 const char *pszAuthority = GetAuthorityName("GEOGCS");
1130
1131 if( pszAuthority && EQUAL(pszAuthority,"EPSG") )
1132 {
1133 int nGCS_EPSG = atoi(GetAuthorityCode("GEOGCS"));
1134 int i;
1135
1136 for( i = 0; asDatums[i].nEPSGCode != 0; i++ )
1137 {
1138 if( asDatums[i].nEPSGCode == nGCS_EPSG )
1139 {
1140 strncpy( szEarthModel, asDatums[i].pszPCIDatum, 5 );
1141 break;
1142 }
1143 }
1144 }
1145 }
1146
1147 /* -------------------------------------------------------------------- */
1148 /* If we haven't found something yet, try translating the */
1149 /* ellipsoid. */
1150 /* -------------------------------------------------------------------- */
1151 if( szEarthModel[0] == '\0' )
1152 {
1153 double dfSemiMajor = GetSemiMajor();
1154 double dfInvFlattening = GetInvFlattening();
1155
1156 const PCIDatums *pasDatum = asEllips;
1157
1158 while ( pasDatum->pszPCIDatum )
1159 {
1160 double dfSM;
1161 double dfIF;
1162
1163 if ( OSRGetEllipsoidInfo( pasDatum->nEPSGCode, NULL,
1164 &dfSM, &dfIF ) == OGRERR_NONE
1165 && CPLIsEqual( dfSemiMajor, dfSM )
1166 && CPLIsEqual( dfInvFlattening, dfIF ) )
1167 {
1168 CPLPrintStringFill( szEarthModel, pasDatum->pszPCIDatum, 4 );
1169 break;
1170 }
1171
1172 pasDatum++;
1173 }
1174
1175 // Try to find in pci_ellips.txt
1176 if( szEarthModel[0] == '\0' )
1177 {
1178 const char *pszCSV = CSVFilename( "pci_ellips.txt" );
1179 FILE *fp = NULL;
1180 double dfSemiMinor = OSRCalcSemiMinorFromInvFlattening(dfSemiMajor, dfInvFlattening);
1181
1182 if( pszCSV )
1183 fp = VSIFOpen( pszCSV, "r" );
1184
1185 if( fp != NULL )
1186 {
1187 char **papszLineItems = NULL;
1188
1189 while( (papszLineItems = CSVReadParseLine( fp )) != NULL )
1190 {
1191 if( CSLCount(papszLineItems) >= 4
1192 && CPLIsEqual(dfSemiMajor,CPLAtof(papszLineItems[2]))
1193 && CPLIsEqual(dfSemiMinor,CPLAtof(papszLineItems[3])) )
1194 {
1195 strncpy( szEarthModel, papszLineItems[0], 5 );
1196 break;
1197 }
1198
1199 CSLDestroy( papszLineItems );
1200 }
1201
1202 CSLDestroy( papszLineItems );
1203 VSIFClose( fp );
1204 }
1205 }
1206
1207 // custom ellipsoid parameters
1208 if( szEarthModel[0] == '\0' )
1209 {
1210 CPLPrintStringFill( szEarthModel, "E999", 4 );
1211 (*ppadfPrjParams)[0] = dfSemiMajor;
1212 (*ppadfPrjParams)[1] = OSRCalcSemiMinorFromInvFlattening(dfSemiMajor, dfInvFlattening);
1213 }
1214 }
1215
1216 /* -------------------------------------------------------------------- */
1217 /* If we have a non-parameteric ellipsoid, scan the */
1218 /* pci_datum.txt for a match. */
1219 /* -------------------------------------------------------------------- */
1220 if( szEarthModel[0] == 'E'
1221 && !EQUAL(szEarthModel,"E999")
1222 && pszDatum != NULL )
1223 {
1224 const char *pszDatumCSV = CSVFilename( "pci_datum.txt" );
1225 FILE *fp = NULL;
1226 double adfTOWGS84[7];
1227 int bHaveTOWGS84;
1228
1229 bHaveTOWGS84 = (GetTOWGS84( adfTOWGS84, 7 ) == OGRERR_NONE);
1230
1231 if( pszDatumCSV )
1232 fp = VSIFOpen( pszDatumCSV, "r" );
1233
1234 if( fp != NULL )
1235 {
1236 char **papszLineItems = NULL;
1237
1238 while( (papszLineItems = CSVReadParseLine( fp )) != NULL )
1239 {
1240 // Compare based on datum name. This is mostly for
1241 // PCI round-tripping. We won't usually get exact matches
1242 // from other sources.
1243 if( CSLCount(papszLineItems) > 3
1244 && EQUAL(papszLineItems[1],pszDatum)
1245 && EQUAL(papszLineItems[2],szEarthModel) )
1246 {
1247 strncpy( szEarthModel, papszLineItems[0], 5 );
1248 break;
1249 }
1250
1251 int bTOWGS84Match = bHaveTOWGS84;
1252
1253 if( CSLCount(papszLineItems) < 11 )
1254 bTOWGS84Match = FALSE;
1255
1256 if( bTOWGS84Match
1257 && (!CPLIsEqual(adfTOWGS84[0],CPLAtof(papszLineItems[3]))
1258 || !CPLIsEqual(adfTOWGS84[1],CPLAtof(papszLineItems[4]))
1259 || !CPLIsEqual(adfTOWGS84[2],CPLAtof(papszLineItems[5]))))
1260 bTOWGS84Match = FALSE;
1261
1262 if( bTOWGS84Match && CSLCount(papszLineItems) >= 15
1263 && (!CPLIsEqual(adfTOWGS84[3],CPLAtof(papszLineItems[11]))
1264 || !CPLIsEqual(adfTOWGS84[4],CPLAtof(papszLineItems[12]))
1265 || !CPLIsEqual(adfTOWGS84[5],CPLAtof(papszLineItems[13]))))
1266 bTOWGS84Match = FALSE;
1267
1268 if( bTOWGS84Match && CSLCount(papszLineItems) >= 15 )
1269 {
1270 double dfScale = CPLAtof(papszLineItems[14]);
1271
1272 // convert to parts per million if is a 1 based scaling.
1273 if( dfScale >= 0.999 && dfScale <= 1.001 )
1274 dfScale = (dfScale-1.0) * 1000000.0;
1275
1276 if( !CPLIsEqual(adfTOWGS84[6],dfScale) )
1277 bTOWGS84Match = FALSE;
1278 }
1279
1280 if( bTOWGS84Match && CSLCount(papszLineItems) < 15
1281 && (!CPLIsEqual(adfTOWGS84[3],0.0)
1282 || !CPLIsEqual(adfTOWGS84[4],0.0)
1283 || !CPLIsEqual(adfTOWGS84[5],0.0)
1284 || !CPLIsEqual(adfTOWGS84[6],0.0)) )
1285 bTOWGS84Match = FALSE;
1286
1287 if( bTOWGS84Match )
1288 {
1289 strncpy( szEarthModel, papszLineItems[0], 5 );
1290 break;
1291 }
1292
1293 CSLDestroy( papszLineItems );
1294 }
1295
1296 CSLDestroy( papszLineItems );
1297 VSIFClose( fp );
1298 }
1299 }
1300
1301 CPLPrintStringFill( szProj + 12, szEarthModel, 4 );
1302
1303 CPLDebug( "OSR_PCI", "Translated as '%s'", szProj );
1304
1305 /* -------------------------------------------------------------------- */
1306 /* Translate the linear units. */
1307 /* -------------------------------------------------------------------- */
1308 const char *pszUnits;
1309
1310 if( EQUALN( szProj, "LONG/LAT", 8 ) )
1311 pszUnits = "DEGREE";
1312 else
1313 pszUnits = "METRE";
1314
1315 /* -------------------------------------------------------------------- */
1316 /* Report results. */
1317 /* -------------------------------------------------------------------- */
1318 szProj[16] = '\0';
1319 *ppszProj = CPLStrdup( szProj );
1320
1321 *ppszUnits = CPLStrdup( pszUnits );
1322
1323 return OGRERR_NONE;
1324 }
1325