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