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