1 : /******************************************************************************
2 : * $Id: ogr_fromepsg.cpp 19236 2010-03-28 16:40:33Z rouault $
3 : *
4 : * Project: OpenGIS Simple Features Reference Implementation
5 : * Purpose: Generate an OGRSpatialReference object based on an EPSG
6 : * PROJCS, or GEOGCS code.
7 : * Author: Frank Warmerdam, warmerdam@pobox.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2000, Frank Warmerdam
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_csv.h"
34 :
35 : CPL_CVSID("$Id: ogr_fromepsg.cpp 19236 2010-03-28 16:40:33Z rouault $");
36 :
37 : #ifndef PI
38 : # define PI 3.14159265358979323846
39 : #endif
40 :
41 : static const char *papszDatumEquiv[] =
42 : {
43 : "Militar_Geographische_Institut",
44 : "Militar_Geographische_Institute",
45 : "World_Geodetic_System_1984",
46 : "WGS_1984",
47 : "WGS_72_Transit_Broadcast_Ephemeris",
48 : "WGS_1972_Transit_Broadcast_Ephemeris",
49 : "World_Geodetic_System_1972",
50 : "WGS_1972",
51 : "European_Terrestrial_Reference_System_89",
52 : "European_Reference_System_1989",
53 : NULL
54 : };
55 :
56 : /************************************************************************/
57 : /* OGREPSGDatumNameMassage() */
58 : /* */
59 : /* Massage an EPSG datum name into WMT format. Also transform */
60 : /* specific exception cases into WKT versions. */
61 : /************************************************************************/
62 :
63 40661 : void OGREPSGDatumNameMassage( char ** ppszDatum )
64 :
65 : {
66 : int i, j;
67 40661 : char *pszDatum = *ppszDatum;
68 :
69 40661 : if (pszDatum[0] == '\0')
70 0 : return;
71 :
72 : /* -------------------------------------------------------------------- */
73 : /* Translate non-alphanumeric values to underscores. */
74 : /* -------------------------------------------------------------------- */
75 898810 : for( i = 0; pszDatum[i] != '\0'; i++ )
76 : {
77 858149 : if( !(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z')
78 : && !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z')
79 : && !(pszDatum[i] >= '0' && pszDatum[i] <= '9') )
80 : {
81 104633 : pszDatum[i] = '_';
82 : }
83 : }
84 :
85 : /* -------------------------------------------------------------------- */
86 : /* Remove repeated and trailing underscores. */
87 : /* -------------------------------------------------------------------- */
88 858149 : for( i = 1, j = 0; pszDatum[i] != '\0'; i++ )
89 : {
90 817488 : if( pszDatum[j] == '_' && pszDatum[i] == '_' )
91 7919 : continue;
92 :
93 809569 : pszDatum[++j] = pszDatum[i];
94 : }
95 40661 : if( pszDatum[j] == '_' )
96 6221 : pszDatum[j] = '\0';
97 : else
98 34440 : pszDatum[j+1] = '\0';
99 :
100 : /* -------------------------------------------------------------------- */
101 : /* Search for datum equivelences. Specific massaged names get */
102 : /* mapped to OpenGIS specified names. */
103 : /* -------------------------------------------------------------------- */
104 239031 : for( i = 0; papszDatumEquiv[i] != NULL; i += 2 )
105 : {
106 199801 : if( EQUAL(*ppszDatum,papszDatumEquiv[i]) )
107 : {
108 1431 : CPLFree( *ppszDatum );
109 1431 : *ppszDatum = CPLStrdup( papszDatumEquiv[i+1] );
110 1431 : break;
111 : }
112 : }
113 : }
114 :
115 : /************************************************************************/
116 : /* EPSGAngleStringToDD() */
117 : /* */
118 : /* Convert an angle in the specified units to decimal degrees. */
119 : /************************************************************************/
120 :
121 : static double
122 25131 : EPSGAngleStringToDD( const char * pszAngle, int nUOMAngle )
123 :
124 : {
125 : double dfAngle;
126 :
127 25131 : if( nUOMAngle == 9110 ) /* DDD.MMSSsss */
128 : {
129 : char *pszDecimal;
130 :
131 7739 : dfAngle = ABS(atoi(pszAngle));
132 7739 : pszDecimal = (char *) strchr(pszAngle,'.');
133 7739 : if( pszDecimal != NULL && strlen(pszDecimal) > 1 )
134 : {
135 : char szMinutes[3];
136 : char szSeconds[64];
137 :
138 6143 : szMinutes[0] = pszDecimal[1];
139 9139 : if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' )
140 2996 : szMinutes[1] = pszDecimal[2];
141 : else
142 3147 : szMinutes[1] = '0';
143 :
144 6143 : szMinutes[2] = '\0';
145 6143 : dfAngle += atoi(szMinutes) / 60.0;
146 :
147 6143 : if( strlen(pszDecimal) > 3 )
148 : {
149 728 : szSeconds[0] = pszDecimal[3];
150 1399 : if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' )
151 : {
152 671 : szSeconds[1] = pszDecimal[4];
153 671 : szSeconds[2] = '.';
154 671 : strncpy( szSeconds+3, pszDecimal + 5, sizeof(szSeconds)-3 );
155 671 : szSeconds[sizeof(szSeconds)-1] = 0;
156 : }
157 : else
158 : {
159 57 : szSeconds[1] = '0';
160 57 : szSeconds[2] = '\0';
161 : }
162 728 : dfAngle += CPLAtof(szSeconds) / 3600.0;
163 : }
164 : }
165 :
166 7739 : if( pszAngle[0] == '-' )
167 2533 : dfAngle *= -1;
168 : }
169 17602 : else if( nUOMAngle == 9105 || nUOMAngle == 9106 ) /* grad */
170 : {
171 210 : dfAngle = 180 * (CPLAtof(pszAngle ) / 200);
172 : }
173 17182 : else if( nUOMAngle == 9101 ) /* radians */
174 : {
175 0 : dfAngle = 180 * (CPLAtof(pszAngle ) / PI);
176 : }
177 17182 : else if( nUOMAngle == 9103 ) /* arc-minute */
178 : {
179 0 : dfAngle = CPLAtof(pszAngle) / 60;
180 : }
181 17182 : else if( nUOMAngle == 9104 ) /* arc-second */
182 : {
183 0 : dfAngle = CPLAtof(pszAngle) / 3600;
184 : }
185 : else /* decimal degrees ... some cases missing but seeminly never used */
186 : {
187 17182 : CPLAssert( nUOMAngle == 9102 || nUOMAngle == 0 );
188 :
189 17182 : dfAngle = CPLAtof(pszAngle );
190 : }
191 :
192 25131 : return( dfAngle );
193 : }
194 :
195 : /************************************************************************/
196 : /* EPSGGetUOMAngleInfo() */
197 : /************************************************************************/
198 :
199 : int EPSGGetUOMAngleInfo( int nUOMAngleCode,
200 : char **ppszUOMName,
201 8125 : double * pdfInDegrees )
202 :
203 : {
204 8125 : const char *pszUOMName = NULL;
205 8125 : double dfInDegrees = 1.0;
206 8125 : const char *pszFilename = CSVFilename( "unit_of_measure.csv" );
207 : char szSearchKey[24];
208 :
209 8125 : sprintf( szSearchKey, "%d", nUOMAngleCode );
210 : pszUOMName = CSVGetField( pszFilename,
211 : "UOM_CODE", szSearchKey, CC_Integer,
212 8125 : "UNIT_OF_MEAS_NAME" );
213 :
214 : /* -------------------------------------------------------------------- */
215 : /* If the file is found, read from there. Note that FactorC is */
216 : /* an empty field for any of the DMS style formats, and in this */
217 : /* case we really want to return the default InDegrees value */
218 : /* (1.0) from above. */
219 : /* -------------------------------------------------------------------- */
220 8125 : if( pszUOMName != NULL )
221 : {
222 : double dfFactorB, dfFactorC;
223 :
224 : dfFactorB =
225 : CPLAtof(CSVGetField( pszFilename,
226 : "UOM_CODE", szSearchKey, CC_Integer,
227 8125 : "FACTOR_B" ));
228 :
229 : dfFactorC =
230 : CPLAtof(CSVGetField( pszFilename,
231 : "UOM_CODE", szSearchKey, CC_Integer,
232 8125 : "FACTOR_C" ));
233 :
234 8125 : if( dfFactorC != 0.0 )
235 8011 : dfInDegrees = (dfFactorB / dfFactorC) * (180.0 / PI);
236 :
237 : /* We do a special override of some of the DMS formats name */
238 8125 : if( nUOMAngleCode == 9102 || nUOMAngleCode == 9107
239 : || nUOMAngleCode == 9108 || nUOMAngleCode == 9110
240 : || nUOMAngleCode == 9122 )
241 8065 : pszUOMName = "degree";
242 :
243 : // For some reason, (FactorB) is not very precise in EPSG, use
244 : // a more exact form for grads.
245 8125 : if( nUOMAngleCode == 9105 )
246 60 : dfInDegrees = 180.0 / 200.0;
247 : }
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Otherwise handle a few well known units directly. */
251 : /* -------------------------------------------------------------------- */
252 : else
253 : {
254 0 : switch( nUOMAngleCode )
255 : {
256 : case 9101:
257 0 : pszUOMName = "radian";
258 0 : dfInDegrees = 180.0 / PI;
259 0 : break;
260 :
261 : case 9102:
262 : case 9107:
263 : case 9108:
264 : case 9110:
265 : case 9122:
266 0 : pszUOMName = "degree";
267 0 : dfInDegrees = 1.0;
268 0 : break;
269 :
270 : case 9103:
271 0 : pszUOMName = "arc-minute";
272 0 : dfInDegrees = 1 / 60.0;
273 0 : break;
274 :
275 : case 9104:
276 0 : pszUOMName = "arc-second";
277 0 : dfInDegrees = 1 / 3600.0;
278 0 : break;
279 :
280 : case 9105:
281 0 : pszUOMName = "grad";
282 0 : dfInDegrees = 180.0 / 200.0;
283 0 : break;
284 :
285 : case 9106:
286 0 : pszUOMName = "gon";
287 0 : dfInDegrees = 180.0 / 200.0;
288 0 : break;
289 :
290 : case 9109:
291 0 : pszUOMName = "microradian";
292 0 : dfInDegrees = 180.0 / (3.14159265358979 * 1000000.0);
293 0 : break;
294 :
295 : default:
296 0 : return FALSE;
297 : }
298 : }
299 :
300 : /* -------------------------------------------------------------------- */
301 : /* Return to caller. */
302 : /* -------------------------------------------------------------------- */
303 8125 : if( ppszUOMName != NULL )
304 : {
305 8125 : if( pszUOMName != NULL )
306 8125 : *ppszUOMName = CPLStrdup( pszUOMName );
307 : else
308 0 : *ppszUOMName = NULL;
309 : }
310 :
311 8125 : if( pdfInDegrees != NULL )
312 8125 : *pdfInDegrees = dfInDegrees;
313 :
314 8125 : return( TRUE );
315 : }
316 :
317 : /************************************************************************/
318 : /* EPSGGetUOMLengthInfo() */
319 : /* */
320 : /* Note: This function should eventually also know how to */
321 : /* lookup length aliases in the UOM_LE_ALIAS table. */
322 : /************************************************************************/
323 :
324 : static int
325 : EPSGGetUOMLengthInfo( int nUOMLengthCode,
326 : char **ppszUOMName,
327 29411 : double * pdfInMeters )
328 :
329 : {
330 : char **papszUnitsRecord;
331 : char szSearchKey[24];
332 : int iNameField;
333 :
334 : #define UOM_FILENAME CSVFilename( "unit_of_measure.csv" )
335 :
336 : /* -------------------------------------------------------------------- */
337 : /* We short cut meter to save work in the most common case. */
338 : /* -------------------------------------------------------------------- */
339 29411 : if( nUOMLengthCode == 9001 )
340 : {
341 26174 : if( ppszUOMName != NULL )
342 6063 : *ppszUOMName = CPLStrdup( "metre" );
343 26174 : if( pdfInMeters != NULL )
344 26174 : *pdfInMeters = 1.0;
345 :
346 26174 : return TRUE;
347 : }
348 :
349 : /* -------------------------------------------------------------------- */
350 : /* Search the units database for this unit. If we don't find */
351 : /* it return failure. */
352 : /* -------------------------------------------------------------------- */
353 3237 : sprintf( szSearchKey, "%d", nUOMLengthCode );
354 : papszUnitsRecord =
355 3237 : CSVScanFileByName( UOM_FILENAME, "UOM_CODE", szSearchKey, CC_Integer );
356 :
357 3237 : if( papszUnitsRecord == NULL )
358 0 : return FALSE;
359 :
360 : /* -------------------------------------------------------------------- */
361 : /* Get the name, if requested. */
362 : /* -------------------------------------------------------------------- */
363 3237 : if( ppszUOMName != NULL )
364 : {
365 1026 : iNameField = CSVGetFileFieldId( UOM_FILENAME, "UNIT_OF_MEAS_NAME" );
366 1026 : *ppszUOMName = CPLStrdup( CSLGetField(papszUnitsRecord, iNameField) );
367 : }
368 :
369 : /* -------------------------------------------------------------------- */
370 : /* Get the A and B factor fields, and create the multiplicative */
371 : /* factor. */
372 : /* -------------------------------------------------------------------- */
373 3237 : if( pdfInMeters != NULL )
374 : {
375 : int iBFactorField, iCFactorField;
376 :
377 3237 : iBFactorField = CSVGetFileFieldId( UOM_FILENAME, "FACTOR_B" );
378 3237 : iCFactorField = CSVGetFileFieldId( UOM_FILENAME, "FACTOR_C" );
379 :
380 3237 : if( CPLAtof(CSLGetField(papszUnitsRecord, iCFactorField)) > 0.0 )
381 : *pdfInMeters = CPLAtof(CSLGetField(papszUnitsRecord,iBFactorField))
382 3237 : / CPLAtof(CSLGetField(papszUnitsRecord, iCFactorField));
383 : else
384 0 : *pdfInMeters = 0.0;
385 : }
386 :
387 3237 : return( TRUE );
388 : }
389 :
390 : /************************************************************************/
391 : /* EPSGGetWGS84Transform() */
392 : /* */
393 : /* The following code attempts to find a bursa-wolf */
394 : /* transformation from this GeogCS to WGS84 (4326). */
395 : /* */
396 : /* Faults: */
397 : /* o I think there are codes other than 9603 and 9607 that */
398 : /* return compatible, or easily transformed parameters. */
399 : /* o Only the first path from the given GeogCS is checked due */
400 : /* to limitations in the CSV API. */
401 : /************************************************************************/
402 :
403 9033 : int EPSGGetWGS84Transform( int nGeogCS, double *padfTransform )
404 :
405 : {
406 : int nMethodCode, iDXField, iField;
407 : char szCode[32];
408 : const char *pszFilename;
409 : char **papszLine;
410 :
411 : /* -------------------------------------------------------------------- */
412 : /* Fetch the line from the GCS table. */
413 : /* -------------------------------------------------------------------- */
414 9033 : pszFilename = CSVFilename("gcs.override.csv");
415 9033 : sprintf( szCode, "%d", nGeogCS );
416 : papszLine = CSVScanFileByName( pszFilename,
417 : "COORD_REF_SYS_CODE",
418 9033 : szCode, CC_Integer );
419 9033 : if( papszLine == NULL )
420 : {
421 9033 : pszFilename = CSVFilename("gcs.csv");
422 9033 : sprintf( szCode, "%d", nGeogCS );
423 : papszLine = CSVScanFileByName( pszFilename,
424 : "COORD_REF_SYS_CODE",
425 9033 : szCode, CC_Integer );
426 : }
427 :
428 9033 : if( papszLine == NULL )
429 0 : return FALSE;
430 :
431 : /* -------------------------------------------------------------------- */
432 : /* Verify that the method code is one of our accepted ones. */
433 : /* -------------------------------------------------------------------- */
434 : nMethodCode =
435 : atoi(CSLGetField( papszLine,
436 : CSVGetFileFieldId(pszFilename,
437 9033 : "COORD_OP_METHOD_CODE")));
438 9033 : if( nMethodCode != 9603 && nMethodCode != 9607 && nMethodCode != 9606 )
439 2442 : return FALSE;
440 :
441 : /* -------------------------------------------------------------------- */
442 : /* Fetch the transformation parameters. */
443 : /* -------------------------------------------------------------------- */
444 6591 : iDXField = CSVGetFileFieldId(pszFilename, "DX");
445 6591 : if (iDXField < 0 || CSLCount(papszLine) < iDXField + 7)
446 0 : return FALSE;
447 :
448 52728 : for( iField = 0; iField < 7; iField++ )
449 46137 : padfTransform[iField] = CPLAtof(papszLine[iDXField+iField]);
450 :
451 : /* -------------------------------------------------------------------- */
452 : /* 9607 - coordinate frame rotation has reverse signs on the */
453 : /* rotational coefficients. Fix up now since we internal */
454 : /* operate according to method 9606 (position vector 7-parameter). */
455 : /* -------------------------------------------------------------------- */
456 6591 : if( nMethodCode == 9607 )
457 : {
458 1087 : padfTransform[3] *= -1;
459 1087 : padfTransform[4] *= -1;
460 1087 : padfTransform[5] *= -1;
461 : }
462 :
463 6591 : return TRUE;
464 : }
465 :
466 : /************************************************************************/
467 : /* EPSGGetPMInfo() */
468 : /* */
469 : /* Get the offset between a given prime meridian and Greenwich */
470 : /* in degrees. */
471 : /************************************************************************/
472 :
473 : static int
474 8125 : EPSGGetPMInfo( int nPMCode, char ** ppszName, double *pdfOffset )
475 :
476 : {
477 : char szSearchKey[24];
478 : int nUOMAngle;
479 :
480 : #define PM_FILENAME CSVFilename("prime_meridian.csv")
481 :
482 : /* -------------------------------------------------------------------- */
483 : /* Use a special short cut for Greenwich, since it is so common. */
484 : /* -------------------------------------------------------------------- */
485 8125 : if( nPMCode == 7022 /* PM_Greenwich */ )
486 : {
487 0 : if( pdfOffset != NULL )
488 0 : *pdfOffset = 0.0;
489 0 : if( ppszName != NULL )
490 0 : *ppszName = CPLStrdup( "Greenwich" );
491 0 : return TRUE;
492 : }
493 :
494 : /* -------------------------------------------------------------------- */
495 : /* Search the database for the corresponding datum code. */
496 : /* -------------------------------------------------------------------- */
497 8125 : sprintf( szSearchKey, "%d", nPMCode );
498 :
499 : nUOMAngle =
500 : atoi(CSVGetField( PM_FILENAME,
501 : "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer,
502 8125 : "UOM_CODE" ) );
503 8125 : if( nUOMAngle < 1 )
504 0 : return FALSE;
505 :
506 : /* -------------------------------------------------------------------- */
507 : /* Get the PM offset. */
508 : /* -------------------------------------------------------------------- */
509 8125 : if( pdfOffset != NULL )
510 : {
511 : *pdfOffset =
512 : EPSGAngleStringToDD(
513 : CSVGetField( PM_FILENAME,
514 : "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer,
515 : "GREENWICH_LONGITUDE" ),
516 8125 : nUOMAngle );
517 : }
518 :
519 : /* -------------------------------------------------------------------- */
520 : /* Get the name, if requested. */
521 : /* -------------------------------------------------------------------- */
522 8125 : if( ppszName != NULL )
523 : *ppszName =
524 : CPLStrdup(
525 : CSVGetField( PM_FILENAME,
526 : "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer,
527 8125 : "PRIME_MERIDIAN_NAME" ));
528 :
529 8125 : return( TRUE );
530 : }
531 :
532 : /************************************************************************/
533 : /* EPSGGetGCSInfo() */
534 : /* */
535 : /* Fetch the datum, and prime meridian related to a particular */
536 : /* GCS. */
537 : /************************************************************************/
538 :
539 : static int
540 : EPSGGetGCSInfo( int nGCSCode, char ** ppszName,
541 : int * pnDatum, char **ppszDatumName,
542 : int * pnPM, int *pnEllipsoid, int *pnUOMAngle,
543 16262 : int * pnCoordSysCode )
544 :
545 : {
546 : char szSearchKey[24];
547 : int nDatum, nPM, nUOMAngle, nEllipsoid;
548 : const char *pszFilename;
549 :
550 :
551 : /* -------------------------------------------------------------------- */
552 : /* Search the database for the corresponding datum code. */
553 : /* -------------------------------------------------------------------- */
554 16262 : pszFilename = CSVFilename("gcs.override.csv");
555 16262 : sprintf( szSearchKey, "%d", nGCSCode );
556 :
557 : nDatum = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
558 : szSearchKey, CC_Integer,
559 16262 : "DATUM_CODE" ) );
560 :
561 16262 : if( nDatum < 1 )
562 : {
563 16262 : pszFilename = CSVFilename("gcs.csv");
564 16262 : sprintf( szSearchKey, "%d", nGCSCode );
565 :
566 : nDatum = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
567 : szSearchKey, CC_Integer,
568 16262 : "DATUM_CODE" ) );
569 : }
570 :
571 16262 : if( nDatum < 1 )
572 7101 : return FALSE;
573 :
574 9161 : if( pnDatum != NULL )
575 8125 : *pnDatum = nDatum;
576 :
577 : /* -------------------------------------------------------------------- */
578 : /* Get the PM. */
579 : /* -------------------------------------------------------------------- */
580 : nPM = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
581 : szSearchKey, CC_Integer,
582 9161 : "PRIME_MERIDIAN_CODE" ) );
583 :
584 9161 : if( nPM < 1 )
585 0 : return FALSE;
586 :
587 9161 : if( pnPM != NULL )
588 8125 : *pnPM = nPM;
589 :
590 : /* -------------------------------------------------------------------- */
591 : /* Get the Ellipsoid. */
592 : /* -------------------------------------------------------------------- */
593 : nEllipsoid = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
594 : szSearchKey, CC_Integer,
595 9161 : "ELLIPSOID_CODE" ) );
596 :
597 9161 : if( nEllipsoid < 1 )
598 0 : return FALSE;
599 :
600 9161 : if( pnEllipsoid != NULL )
601 8125 : *pnEllipsoid = nEllipsoid;
602 :
603 : /* -------------------------------------------------------------------- */
604 : /* Get the angular units. */
605 : /* -------------------------------------------------------------------- */
606 : nUOMAngle = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
607 : szSearchKey, CC_Integer,
608 9161 : "UOM_CODE" ) );
609 :
610 9161 : if( nUOMAngle < 1 )
611 0 : return FALSE;
612 :
613 9161 : if( pnUOMAngle != NULL )
614 8125 : *pnUOMAngle = nUOMAngle;
615 :
616 : /* -------------------------------------------------------------------- */
617 : /* Get the name, if requested. */
618 : /* -------------------------------------------------------------------- */
619 9161 : if( ppszName != NULL )
620 : *ppszName =
621 : CPLStrdup(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
622 : szSearchKey, CC_Integer,
623 8125 : "COORD_REF_SYS_NAME" ));
624 :
625 : /* -------------------------------------------------------------------- */
626 : /* Get the datum name, if requested. */
627 : /* -------------------------------------------------------------------- */
628 9161 : if( ppszDatumName != NULL )
629 : *ppszDatumName =
630 : CPLStrdup(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
631 : szSearchKey, CC_Integer,
632 8125 : "DATUM_NAME" ));
633 :
634 : /* -------------------------------------------------------------------- */
635 : /* Get the CoordSysCode */
636 : /* -------------------------------------------------------------------- */
637 : int nCSC;
638 :
639 : nCSC = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
640 : szSearchKey, CC_Integer,
641 9161 : "COORD_SYS_CODE" ) );
642 :
643 9161 : if( pnCoordSysCode != NULL )
644 8125 : *pnCoordSysCode = nCSC;
645 :
646 9161 : return( TRUE );
647 : }
648 :
649 : /************************************************************************/
650 : /* OSRGetEllipsoidInfo() */
651 : /************************************************************************/
652 :
653 : /**
654 : * Fetch info about an ellipsoid.
655 : *
656 : * This helper function will return ellipsoid parameters corresponding to EPSG
657 : * code provided. Axes are always returned in meters. Semi major computed
658 : * based on inverse flattening where that is provided.
659 : *
660 : * @param nCode EPSG code of the requested ellipsoid
661 : *
662 : * @param ppszName pointer to string where ellipsoid name will be returned. It
663 : * is caller responsibility to free this string after using with CPLFree().
664 : *
665 : * @param pdfSemiMajor pointer to variable where semi major axis will be
666 : * returned.
667 : *
668 : * @param pdfInvFlattening pointer to variable where inverse flattening will
669 : * be returned.
670 : *
671 : * @return OGRERR_NONE on success or an error code in case of failure.
672 : **/
673 :
674 : OGRErr
675 : OSRGetEllipsoidInfo( int nCode, char ** ppszName,
676 8148 : double * pdfSemiMajor, double * pdfInvFlattening )
677 :
678 : {
679 : char szSearchKey[24];
680 8148 : double dfSemiMajor, dfToMeters = 1.0;
681 : int nUOMLength;
682 :
683 : /* -------------------------------------------------------------------- */
684 : /* Get the semi major axis. */
685 : /* -------------------------------------------------------------------- */
686 8148 : snprintf( szSearchKey, sizeof(szSearchKey), "%d", nCode );
687 8148 : szSearchKey[sizeof(szSearchKey) - 1] = '\n';
688 :
689 : dfSemiMajor =
690 : CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ),
691 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
692 8148 : "SEMI_MAJOR_AXIS" ) );
693 8148 : if( dfSemiMajor == 0.0 )
694 0 : return OGRERR_UNSUPPORTED_SRS;
695 :
696 : /* -------------------------------------------------------------------- */
697 : /* Get the translation factor into meters. */
698 : /* -------------------------------------------------------------------- */
699 : nUOMLength = atoi(CSVGetField( CSVFilename("ellipsoid.csv" ),
700 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
701 8148 : "UOM_CODE" ));
702 8148 : EPSGGetUOMLengthInfo( nUOMLength, NULL, &dfToMeters );
703 :
704 8148 : dfSemiMajor *= dfToMeters;
705 :
706 8148 : if( pdfSemiMajor != NULL )
707 8148 : *pdfSemiMajor = dfSemiMajor;
708 :
709 : /* -------------------------------------------------------------------- */
710 : /* Get the semi-minor if requested. If the Semi-minor axis */
711 : /* isn't available, compute it based on the inverse flattening. */
712 : /* -------------------------------------------------------------------- */
713 8148 : if( pdfInvFlattening != NULL )
714 : {
715 : *pdfInvFlattening =
716 : CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ),
717 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
718 8148 : "INV_FLATTENING" ));
719 :
720 8148 : if( *pdfInvFlattening == 0.0 )
721 : {
722 : double dfSemiMinor;
723 :
724 : dfSemiMinor =
725 : CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ),
726 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
727 837 : "SEMI_MINOR_AXIS" )) * dfToMeters;
728 :
729 1645 : if( dfSemiMajor != 0.0 && dfSemiMajor != dfSemiMinor )
730 : *pdfInvFlattening =
731 808 : -1.0 / (dfSemiMinor/dfSemiMajor - 1.0);
732 : else
733 29 : *pdfInvFlattening = 0.0;
734 : }
735 : }
736 :
737 : /* -------------------------------------------------------------------- */
738 : /* Get the name, if requested. */
739 : /* -------------------------------------------------------------------- */
740 8148 : if( ppszName != NULL )
741 : *ppszName =
742 : CPLStrdup(CSVGetField( CSVFilename("ellipsoid.csv" ),
743 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
744 8142 : "ELLIPSOID_NAME" ));
745 :
746 8148 : return OGRERR_NONE;
747 : }
748 :
749 : #define NatOriginLat 8801
750 : #define NatOriginLong 8802
751 : #define NatOriginScaleFactor 8805
752 : #define FalseEasting 8806
753 : #define FalseNorthing 8807
754 : #define ProjCenterLat 8811
755 : #define ProjCenterLong 8812
756 : #define Azimuth 8813
757 : #define AngleRectifiedToSkewedGrid 8814
758 : #define InitialLineScaleFactor 8815
759 : #define ProjCenterEasting 8816
760 : #define ProjCenterNorthing 8817
761 : #define PseudoStdParallelLat 8818
762 : #define PseudoStdParallelScaleFactor 8819
763 : #define FalseOriginLat 8821
764 : #define FalseOriginLong 8822
765 : #define StdParallel1Lat 8823
766 : #define StdParallel2Lat 8824
767 : #define FalseOriginEasting 8826
768 : #define FalseOriginNorthing 8827
769 : #define SphericalOriginLat 8828
770 : #define SphericalOriginLong 8829
771 : #define InitialLongitude 8830
772 : #define ZoneWidth 8831
773 : #define PolarLatStdParallel 8832
774 : #define PolarLongOrigin 8833
775 :
776 : /************************************************************************/
777 : /* EPSGGetProjTRFInfo() */
778 : /* */
779 : /* Transform a PROJECTION_TRF_CODE into a projection method, */
780 : /* and a set of parameters. The parameters identify will */
781 : /* depend on the returned method, but they will all have been */
782 : /* normalized into degrees and meters. */
783 : /************************************************************************/
784 :
785 : static int
786 : EPSGGetProjTRFInfo( int nPCS, int * pnProjMethod,
787 7089 : int *panParmIds, double * padfProjParms )
788 :
789 : {
790 : int nProjMethod, i;
791 : double adfProjParms[7];
792 : char szTRFCode[16];
793 7089 : CPLString osFilename;
794 :
795 : /* -------------------------------------------------------------------- */
796 : /* Get the proj method. If this fails to return a meaningful */
797 : /* number, then the whole function fails. */
798 : /* -------------------------------------------------------------------- */
799 7089 : osFilename = CSVFilename( "pcs.override.csv" );
800 7089 : sprintf( szTRFCode, "%d", nPCS );
801 : nProjMethod =
802 : atoi( CSVGetField( osFilename,
803 : "COORD_REF_SYS_CODE", szTRFCode, CC_Integer,
804 7089 : "COORD_OP_METHOD_CODE" ) );
805 7089 : if( nProjMethod == 0 )
806 : {
807 7083 : osFilename = CSVFilename( "pcs.csv" );
808 7083 : sprintf( szTRFCode, "%d", nPCS );
809 : nProjMethod =
810 : atoi( CSVGetField( osFilename,
811 : "COORD_REF_SYS_CODE", szTRFCode, CC_Integer,
812 7083 : "COORD_OP_METHOD_CODE" ) );
813 : }
814 :
815 7089 : if( nProjMethod == 0 )
816 0 : return FALSE;
817 :
818 : /* -------------------------------------------------------------------- */
819 : /* Get the parameters for this projection. */
820 : /* -------------------------------------------------------------------- */
821 :
822 56712 : for( i = 0; i < 7; i++ )
823 : {
824 : char szParamUOMID[32], szParamValueID[32], szParamCodeID[32];
825 : char *pszValue;
826 : int nUOM;
827 :
828 49623 : sprintf( szParamCodeID, "PARAMETER_CODE_%d", i+1 );
829 49623 : sprintf( szParamUOMID, "PARAMETER_UOM_%d", i+1 );
830 49623 : sprintf( szParamValueID, "PARAMETER_VALUE_%d", i+1 );
831 :
832 49623 : if( panParmIds != NULL )
833 : panParmIds[i] =
834 : atoi(CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode,
835 49623 : CC_Integer, szParamCodeID ));
836 :
837 : nUOM = atoi(CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode,
838 49623 : CC_Integer, szParamUOMID ));
839 : pszValue = CPLStrdup(
840 : CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode,
841 49623 : CC_Integer, szParamValueID ));
842 :
843 : // there is a bug in the EPSG 6.2.2 database for PCS 2935 and 2936
844 : // such that they have foot units for the scale factor. Avoid this.
845 49623 : if( (panParmIds[i] == NatOriginScaleFactor
846 : || panParmIds[i] == InitialLineScaleFactor
847 : || panParmIds[i] == PseudoStdParallelScaleFactor)
848 : && nUOM < 9200 )
849 8 : nUOM = 9201;
850 :
851 66629 : if( nUOM >= 9100 && nUOM < 9200 )
852 17006 : adfProjParms[i] = EPSGAngleStringToDD( pszValue, nUOM );
853 46791 : else if( nUOM > 9000 && nUOM < 9100 )
854 : {
855 : double dfInMeters;
856 :
857 14174 : if( !EPSGGetUOMLengthInfo( nUOM, NULL, &dfInMeters ) )
858 0 : dfInMeters = 1.0;
859 14174 : adfProjParms[i] = CPLAtof(pszValue) * dfInMeters;
860 : }
861 18443 : else if( EQUAL(pszValue,"") ) /* null field */
862 : {
863 12873 : adfProjParms[i] = 0.0;
864 : }
865 : else /* really we should consider looking up other scaling factors */
866 : {
867 5570 : if( nUOM != 9201 )
868 4 : CPLDebug( "OGR", "Non-unity scale factor units!" );
869 5570 : adfProjParms[i] = CPLAtof(pszValue);
870 : }
871 :
872 49623 : CPLFree( pszValue );
873 : }
874 :
875 : /* -------------------------------------------------------------------- */
876 : /* Transfer requested data into passed variables. */
877 : /* -------------------------------------------------------------------- */
878 7089 : if( pnProjMethod != NULL )
879 7089 : *pnProjMethod = nProjMethod;
880 :
881 7089 : if( padfProjParms != NULL )
882 : {
883 56712 : for( i = 0; i < 7; i++ )
884 49623 : padfProjParms[i] = adfProjParms[i];
885 : }
886 :
887 7089 : return TRUE;
888 : }
889 :
890 :
891 : /************************************************************************/
892 : /* EPSGGetPCSInfo() */
893 : /************************************************************************/
894 :
895 : static int
896 : EPSGGetPCSInfo( int nPCSCode, char **ppszEPSGName,
897 : int *pnUOMLengthCode, int *pnUOMAngleCode,
898 7095 : int *pnGeogCS, int *pnTRFCode, int *pnCoordSysCode )
899 :
900 : {
901 : char **papszRecord;
902 : char szSearchKey[24];
903 : const char *pszFilename;
904 :
905 : /* -------------------------------------------------------------------- */
906 : /* Search the units database for this unit. If we don't find */
907 : /* it return failure. */
908 : /* -------------------------------------------------------------------- */
909 7095 : pszFilename = CSVFilename( "pcs.csv" );
910 7095 : sprintf( szSearchKey, "%d", nPCSCode );
911 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
912 7095 : szSearchKey, CC_Integer );
913 :
914 7095 : if( papszRecord == NULL )
915 : {
916 0 : pszFilename = CSVFilename( "pcs.override.csv" );
917 0 : sprintf( szSearchKey, "%d", nPCSCode );
918 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
919 0 : szSearchKey, CC_Integer );
920 :
921 : }
922 :
923 7095 : if( papszRecord == NULL )
924 0 : return FALSE;
925 :
926 : /* -------------------------------------------------------------------- */
927 : /* Get the name, if requested. */
928 : /* -------------------------------------------------------------------- */
929 7095 : if( ppszEPSGName != NULL )
930 : {
931 : CPLString osPCSName =
932 : CSLGetField( papszRecord,
933 : CSVGetFileFieldId(pszFilename,
934 7095 : "COORD_REF_SYS_NAME"));
935 :
936 : const char *pszDeprecated =
937 : CSLGetField( papszRecord,
938 : CSVGetFileFieldId(pszFilename,
939 7095 : "DEPRECATED") );
940 :
941 7095 : if( pszDeprecated != NULL && *pszDeprecated == '1' )
942 521 : osPCSName += " (deprecated)";
943 :
944 7095 : *ppszEPSGName = CPLStrdup(osPCSName);
945 : }
946 :
947 : /* -------------------------------------------------------------------- */
948 : /* Get the UOM Length code, if requested. */
949 : /* -------------------------------------------------------------------- */
950 7095 : if( pnUOMLengthCode != NULL )
951 : {
952 : const char *pszValue;
953 :
954 : pszValue =
955 : CSLGetField( papszRecord,
956 7095 : CSVGetFileFieldId(pszFilename,"UOM_CODE"));
957 7095 : if( atoi(pszValue) > 0 )
958 7095 : *pnUOMLengthCode = atoi(pszValue);
959 : else
960 0 : *pnUOMLengthCode = 0;
961 : }
962 :
963 : /* -------------------------------------------------------------------- */
964 : /* Get the UOM Angle code, if requested. */
965 : /* -------------------------------------------------------------------- */
966 7095 : if( pnUOMAngleCode != NULL )
967 : {
968 : const char *pszValue;
969 :
970 : pszValue =
971 : CSLGetField( papszRecord,
972 7095 : CSVGetFileFieldId(pszFilename,"UOM_ANGLE_CODE") );
973 :
974 7095 : if( atoi(pszValue) > 0 )
975 0 : *pnUOMAngleCode = atoi(pszValue);
976 : else
977 7095 : *pnUOMAngleCode = 0;
978 : }
979 :
980 : /* -------------------------------------------------------------------- */
981 : /* Get the GeogCS (Datum with PM) code, if requested. */
982 : /* -------------------------------------------------------------------- */
983 7095 : if( pnGeogCS != NULL )
984 : {
985 : const char *pszValue;
986 :
987 : pszValue =
988 : CSLGetField( papszRecord,
989 7095 : CSVGetFileFieldId(pszFilename,"SOURCE_GEOGCRS_CODE"));
990 7095 : if( atoi(pszValue) > 0 )
991 7095 : *pnGeogCS = atoi(pszValue);
992 : else
993 0 : *pnGeogCS = 0;
994 : }
995 :
996 : /* -------------------------------------------------------------------- */
997 : /* Get the GeogCS (Datum with PM) code, if requested. */
998 : /* -------------------------------------------------------------------- */
999 7095 : if( pnTRFCode != NULL )
1000 : {
1001 : const char *pszValue;
1002 :
1003 : pszValue =
1004 : CSLGetField( papszRecord,
1005 7095 : CSVGetFileFieldId(pszFilename,"COORD_OP_CODE"));
1006 :
1007 :
1008 7095 : if( atoi(pszValue) > 0 )
1009 7095 : *pnTRFCode = atoi(pszValue);
1010 : else
1011 0 : *pnTRFCode = 0;
1012 : }
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* Get the CoordSysCode */
1016 : /* -------------------------------------------------------------------- */
1017 : int nCSC;
1018 :
1019 : nCSC = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
1020 : szSearchKey, CC_Integer,
1021 7095 : "COORD_SYS_CODE" ) );
1022 :
1023 7095 : if( pnCoordSysCode != NULL )
1024 7095 : *pnCoordSysCode = nCSC;
1025 :
1026 7095 : return TRUE;
1027 : }
1028 :
1029 : /************************************************************************/
1030 : /* SetEPSGAxisInfo() */
1031 : /************************************************************************/
1032 :
1033 : static OGRErr SetEPSGAxisInfo( OGRSpatialReference *poSRS,
1034 : const char *pszTargetKey,
1035 15144 : int nCoordSysCode )
1036 :
1037 : {
1038 : /* -------------------------------------------------------------------- */
1039 : /* Special cases for well known and common values. We short */
1040 : /* circuit these to save time doing file lookups. */
1041 : /* -------------------------------------------------------------------- */
1042 : // Conventional and common Easting/Northing values.
1043 15144 : if( nCoordSysCode >= 4400 && nCoordSysCode <= 4410 )
1044 : {
1045 : return
1046 : poSRS->SetAxes( pszTargetKey,
1047 : "Easting", OAO_East,
1048 2469 : "Northing", OAO_North );
1049 : }
1050 :
1051 : // Conventional and common Easting/Northing values.
1052 12675 : if( nCoordSysCode >= 6400 && nCoordSysCode <= 6423 )
1053 : {
1054 : return
1055 : poSRS->SetAxes( pszTargetKey,
1056 : "Latitude", OAO_North,
1057 8125 : "Longitude", OAO_East );
1058 : }
1059 :
1060 : /* -------------------------------------------------------------------- */
1061 : /* Get the definition from the coordinate_axis.csv file. */
1062 : /* -------------------------------------------------------------------- */
1063 : char **papszRecord;
1064 4550 : char **papszAxis1=NULL, **papszAxis2=NULL;
1065 : char szSearchKey[24];
1066 : const char *pszFilename;
1067 :
1068 4550 : pszFilename = CSVFilename( "coordinate_axis.csv" );
1069 4550 : sprintf( szSearchKey, "%d", nCoordSysCode );
1070 : papszRecord = CSVScanFileByName( pszFilename, "COORD_SYS_CODE",
1071 4550 : szSearchKey, CC_Integer );
1072 :
1073 4550 : if( papszRecord != NULL )
1074 : {
1075 4550 : papszAxis1 = CSLDuplicate( papszRecord );
1076 4550 : papszRecord = CSVGetNextLine( pszFilename );
1077 4550 : if( CSLCount(papszRecord) > 0
1078 : && EQUAL(papszRecord[0],papszAxis1[0]) )
1079 4550 : papszAxis2 = CSLDuplicate( papszRecord );
1080 : }
1081 :
1082 4550 : if( papszAxis2 == NULL )
1083 : {
1084 0 : CSLDestroy( papszAxis1 );
1085 : CPLError( CE_Failure, CPLE_AppDefined,
1086 : "Failed to find entries for COORD_SYS_CODE %d in coordinate_axis.csv",
1087 0 : nCoordSysCode );
1088 0 : return OGRERR_FAILURE;
1089 : }
1090 :
1091 : /* -------------------------------------------------------------------- */
1092 : /* Confirm the records are complete, and work out which columns */
1093 : /* are which. */
1094 : /* -------------------------------------------------------------------- */
1095 : int iAxisOrientationField, iAxisAbbrevField, iAxisOrderField;
1096 : int iAxisNameCodeField;
1097 :
1098 : iAxisOrientationField =
1099 4550 : CSVGetFileFieldId( pszFilename, "coord_axis_orientation" );
1100 : iAxisAbbrevField =
1101 4550 : CSVGetFileFieldId( pszFilename, "coord_axis_abbreviation" );
1102 : iAxisOrderField =
1103 4550 : CSVGetFileFieldId( pszFilename, "coord_axis_order" );
1104 : iAxisNameCodeField =
1105 4550 : CSVGetFileFieldId( pszFilename, "coord_axis_name_code" );
1106 :
1107 : /* Check that all fields are available and that the axis_order field */
1108 : /* is the one with highest index */
1109 4550 : if ( !( iAxisOrientationField >= 0 &&
1110 : iAxisOrientationField < iAxisOrderField &&
1111 : iAxisAbbrevField >= 0 &&
1112 : iAxisAbbrevField < iAxisOrderField &&
1113 : iAxisOrderField >= 0 &&
1114 : iAxisNameCodeField >= 0 &&
1115 : iAxisNameCodeField < iAxisOrderField ) )
1116 : {
1117 0 : CSLDestroy( papszAxis1 );
1118 0 : CSLDestroy( papszAxis2 );
1119 : CPLError( CE_Failure, CPLE_AppDefined,
1120 0 : "coordinate_axis.csv corrupted" );
1121 0 : return OGRERR_FAILURE;
1122 : }
1123 :
1124 4550 : if( CSLCount(papszAxis1) < iAxisOrderField+1
1125 : || CSLCount(papszAxis2) < iAxisOrderField+1 )
1126 : {
1127 0 : CSLDestroy( papszAxis1 );
1128 0 : CSLDestroy( papszAxis2 );
1129 : CPLError( CE_Failure, CPLE_AppDefined,
1130 : "Axis records appear incomplete for COORD_SYS_CODE %d in coordinate_axis.csv",
1131 0 : nCoordSysCode );
1132 0 : return OGRERR_FAILURE;
1133 : }
1134 :
1135 : /* -------------------------------------------------------------------- */
1136 : /* Do we need to switch the axes around? */
1137 : /* -------------------------------------------------------------------- */
1138 4550 : if( atoi(papszAxis2[iAxisOrderField]) < atoi(papszAxis1[iAxisOrderField]) )
1139 : {
1140 2104 : papszRecord = papszAxis1;
1141 2104 : papszAxis1 = papszAxis2;
1142 2104 : papszAxis2 = papszRecord;
1143 : }
1144 :
1145 : /* -------------------------------------------------------------------- */
1146 : /* Work out axis enumeration values. */
1147 : /* -------------------------------------------------------------------- */
1148 4550 : OGRAxisOrientation eOAxis1 = OAO_Other, eOAxis2 = OAO_Other;
1149 : int iAO;
1150 : static int anCodes[7] = { -1, 9907, 9909, 9906, 9908, -1, -1 };
1151 :
1152 36400 : for( iAO = 0; iAO < 7; iAO++ )
1153 : {
1154 31850 : if( EQUAL(papszAxis1[iAxisOrientationField],
1155 : OSRAxisEnumToName((OGRAxisOrientation) iAO)) )
1156 4480 : eOAxis1 = (OGRAxisOrientation) iAO;
1157 31850 : if( EQUAL(papszAxis2[iAxisOrientationField],
1158 : OSRAxisEnumToName((OGRAxisOrientation) iAO)) )
1159 4480 : eOAxis2 = (OGRAxisOrientation) iAO;
1160 :
1161 31850 : if( eOAxis1 == OAO_Other
1162 : && anCodes[iAO] == atoi(papszAxis1[iAxisNameCodeField]) )
1163 70 : eOAxis1 = (OGRAxisOrientation) iAO;
1164 31850 : if( eOAxis2 == OAO_Other
1165 : && anCodes[iAO] == atoi(papszAxis2[iAxisNameCodeField]) )
1166 70 : eOAxis2 = (OGRAxisOrientation) iAO;
1167 : }
1168 :
1169 : /* -------------------------------------------------------------------- */
1170 : /* Work out the axis name. We try to expand the abbreviation */
1171 : /* to a longer name. */
1172 : /* -------------------------------------------------------------------- */
1173 : const char *apszAxisName[2];
1174 4550 : apszAxisName[0] = papszAxis1[iAxisAbbrevField];
1175 4550 : apszAxisName[1] = papszAxis2[iAxisAbbrevField];
1176 :
1177 13650 : for( iAO = 0; iAO < 2; iAO++ )
1178 : {
1179 9100 : if( EQUAL(apszAxisName[iAO],"N") )
1180 260 : apszAxisName[iAO] = "Northing";
1181 8840 : else if( EQUAL(apszAxisName[iAO],"E") )
1182 260 : apszAxisName[iAO] = "Easting";
1183 8580 : else if( EQUAL(apszAxisName[iAO],"S") )
1184 0 : apszAxisName[iAO] = "Southing";
1185 8580 : else if( EQUAL(apszAxisName[iAO],"W") )
1186 0 : apszAxisName[iAO] = "Westing";
1187 : }
1188 :
1189 : /* -------------------------------------------------------------------- */
1190 : /* Set the axes. */
1191 : /* -------------------------------------------------------------------- */
1192 : OGRErr eResult;
1193 : eResult = poSRS->SetAxes( pszTargetKey,
1194 : apszAxisName[0], eOAxis1,
1195 4550 : apszAxisName[1], eOAxis2 );
1196 :
1197 4550 : CSLDestroy( papszAxis1 );
1198 4550 : CSLDestroy( papszAxis2 );
1199 :
1200 4550 : return eResult;
1201 : }
1202 :
1203 : /************************************************************************/
1204 : /* SetEPSGGeogCS() */
1205 : /* */
1206 : /* FLAWS: */
1207 : /* o Units are all hardcoded. */
1208 : /************************************************************************/
1209 :
1210 8131 : static OGRErr SetEPSGGeogCS( OGRSpatialReference * poSRS, int nGeogCS )
1211 :
1212 : {
1213 : int nDatumCode, nPMCode, nUOMAngle, nEllipsoidCode, nCSC;
1214 8131 : char *pszGeogCSName = NULL, *pszDatumName = NULL, *pszEllipsoidName = NULL;
1215 8131 : char *pszPMName = NULL, *pszAngleName = NULL;
1216 : double dfPMOffset, dfSemiMajor, dfInvFlattening, adfBursaTransform[7];
1217 : double dfAngleInDegrees, dfAngleInRadians;
1218 :
1219 8131 : if( !EPSGGetGCSInfo( nGeogCS, &pszGeogCSName,
1220 : &nDatumCode, &pszDatumName,
1221 : &nPMCode, &nEllipsoidCode, &nUOMAngle, &nCSC ) )
1222 6 : return OGRERR_UNSUPPORTED_SRS;
1223 :
1224 8125 : if( !EPSGGetPMInfo( nPMCode, &pszPMName, &dfPMOffset ) )
1225 : {
1226 0 : CPLFree( pszDatumName );
1227 0 : CPLFree( pszGeogCSName );
1228 0 : return OGRERR_UNSUPPORTED_SRS;
1229 : }
1230 :
1231 8125 : OGREPSGDatumNameMassage( &pszDatumName );
1232 :
1233 8125 : if( OSRGetEllipsoidInfo( nEllipsoidCode, &pszEllipsoidName,
1234 : &dfSemiMajor, &dfInvFlattening ) != OGRERR_NONE )
1235 : {
1236 0 : CPLFree( pszDatumName );
1237 0 : CPLFree( pszGeogCSName );
1238 0 : CPLFree( pszPMName );
1239 0 : return OGRERR_UNSUPPORTED_SRS;
1240 : }
1241 :
1242 8125 : if( !EPSGGetUOMAngleInfo( nUOMAngle, &pszAngleName, &dfAngleInDegrees ) )
1243 : {
1244 0 : pszAngleName = CPLStrdup("degree");
1245 0 : dfAngleInDegrees = 1.0;
1246 0 : nUOMAngle = -1;
1247 : }
1248 :
1249 8125 : if( dfAngleInDegrees == 1.0 )
1250 114 : dfAngleInRadians = CPLAtof(SRS_UA_DEGREE_CONV);
1251 : else
1252 8011 : dfAngleInRadians = CPLAtof(SRS_UA_DEGREE_CONV) * dfAngleInDegrees;
1253 :
1254 : poSRS->SetGeogCS( pszGeogCSName, pszDatumName,
1255 : pszEllipsoidName, dfSemiMajor, dfInvFlattening,
1256 : pszPMName, dfPMOffset,
1257 8125 : pszAngleName, dfAngleInRadians );
1258 :
1259 8125 : if( EPSGGetWGS84Transform( nGeogCS, adfBursaTransform ) )
1260 : {
1261 : OGR_SRSNode *poWGS84;
1262 : char szValue[100];
1263 :
1264 6535 : poWGS84 = new OGR_SRSNode( "TOWGS84" );
1265 :
1266 104560 : for( int iCoeff = 0; iCoeff < 7; iCoeff++ )
1267 : {
1268 45745 : sprintf( szValue, "%g", adfBursaTransform[iCoeff] );
1269 45745 : poWGS84->AddChild( new OGR_SRSNode( szValue ) );
1270 : }
1271 :
1272 6535 : poSRS->GetAttrNode( "DATUM" )->AddChild( poWGS84 );
1273 : }
1274 :
1275 8125 : poSRS->SetAuthority( "GEOGCS", "EPSG", nGeogCS );
1276 8125 : poSRS->SetAuthority( "DATUM", "EPSG", nDatumCode );
1277 8125 : poSRS->SetAuthority( "SPHEROID", "EPSG", nEllipsoidCode );
1278 8125 : poSRS->SetAuthority( "PRIMEM", "EPSG", nPMCode );
1279 :
1280 8125 : if( nUOMAngle > 0 )
1281 8125 : poSRS->SetAuthority( "GEOGCS|UNIT", "EPSG", nUOMAngle );
1282 :
1283 8125 : CPLFree( pszAngleName );
1284 8125 : CPLFree( pszDatumName );
1285 8125 : CPLFree( pszEllipsoidName );
1286 8125 : CPLFree( pszGeogCSName );
1287 8125 : CPLFree( pszPMName );
1288 :
1289 : /* -------------------------------------------------------------------- */
1290 : /* Set axes */
1291 : /* -------------------------------------------------------------------- */
1292 8125 : if( nCSC > 0 )
1293 : {
1294 8125 : SetEPSGAxisInfo( poSRS, "GEOGCS", nCSC );
1295 8125 : CPLErrorReset();
1296 : }
1297 :
1298 8125 : return OGRERR_NONE;
1299 : }
1300 :
1301 : /************************************************************************/
1302 : /* OGR_FetchParm() */
1303 : /* */
1304 : /* Fetch a parameter from the parm list, based on it's EPSG */
1305 : /* parameter code. */
1306 : /************************************************************************/
1307 :
1308 : static double OGR_FetchParm( double *padfProjParms, int *panParmIds,
1309 36427 : int nTargetId, double dfFromGreenwich )
1310 :
1311 : {
1312 : int i;
1313 : double dfResult;
1314 :
1315 : /* -------------------------------------------------------------------- */
1316 : /* Set default in meters/degrees. */
1317 : /* -------------------------------------------------------------------- */
1318 36427 : switch( nTargetId )
1319 : {
1320 : case NatOriginScaleFactor:
1321 : case InitialLineScaleFactor:
1322 : case PseudoStdParallelScaleFactor:
1323 5527 : dfResult = 1.0;
1324 5527 : break;
1325 :
1326 : case AngleRectifiedToSkewedGrid:
1327 47 : dfResult = 90.0;
1328 47 : break;
1329 :
1330 : default:
1331 30853 : dfResult = 0.0;
1332 : }
1333 :
1334 : /* -------------------------------------------------------------------- */
1335 : /* Try to find actual value in parameter list. */
1336 : /* -------------------------------------------------------------------- */
1337 113523 : for( i = 0; i < 7; i++ )
1338 : {
1339 113508 : if( panParmIds[i] == nTargetId )
1340 : {
1341 36412 : dfResult = padfProjParms[i];
1342 36412 : break;
1343 : }
1344 : }
1345 :
1346 : /* -------------------------------------------------------------------- */
1347 : /* EPSG longitudes are relative to greenwich. The follow code */
1348 : /* could be used to make them relative to the prime meridian of */
1349 : /* the associated GCS if that was appropriate. However, the */
1350 : /* SetNormProjParm() method expects longitudes relative to */
1351 : /* greenwich, so there is nothing for us to do. */
1352 : /* -------------------------------------------------------------------- */
1353 : #ifdef notdef
1354 : switch( nTargetId )
1355 : {
1356 : case NatOriginLong:
1357 : case ProjCenterLong:
1358 : case FalseOriginLong:
1359 : case SphericalOriginLong:
1360 : case InitialLongitude:
1361 : // Note that the EPSG values are already relative to greenwich.
1362 : // This shift is really making it relative to the provided prime
1363 : // meridian, so that when SetTM() and company the correction back
1364 : // ends up back relative to greenwich.
1365 : dfResult = dfResult + dfFromGreenwich;
1366 : break;
1367 :
1368 : default:
1369 : ;
1370 : }
1371 : #endif
1372 :
1373 36427 : return dfResult;
1374 : }
1375 :
1376 : #define OGR_FP(x) OGR_FetchParm( adfProjParms, anParmIds, (x), \
1377 : dfFromGreenwich )
1378 :
1379 : /************************************************************************/
1380 : /* SetEPSGProjCS() */
1381 : /************************************************************************/
1382 :
1383 7095 : static OGRErr SetEPSGProjCS( OGRSpatialReference * poSRS, int nPCSCode )
1384 :
1385 : {
1386 7095 : int nGCSCode, nUOMAngleCode, nUOMLength, nTRFCode, nProjMethod=0;
1387 7095 : int anParmIds[7], nCSC = 0;
1388 7095 : char *pszPCSName = NULL, *pszUOMLengthName = NULL;
1389 : double adfProjParms[7], dfInMeters, dfFromGreenwich;
1390 : OGRErr nErr;
1391 : OGR_SRSNode *poNode;
1392 :
1393 7095 : if( !EPSGGetPCSInfo( nPCSCode, &pszPCSName, &nUOMLength, &nUOMAngleCode,
1394 : &nGCSCode, &nTRFCode, &nCSC ) )
1395 : {
1396 0 : CPLFree(pszPCSName);
1397 0 : return OGRERR_UNSUPPORTED_SRS;
1398 : }
1399 :
1400 7095 : poSRS->SetNode( "PROJCS", pszPCSName );
1401 :
1402 : /* -------------------------------------------------------------------- */
1403 : /* Set GEOGCS. */
1404 : /* -------------------------------------------------------------------- */
1405 7095 : nErr = SetEPSGGeogCS( poSRS, nGCSCode );
1406 7095 : if( nErr != OGRERR_NONE )
1407 : {
1408 6 : CPLFree(pszPCSName);
1409 6 : return nErr;
1410 : }
1411 :
1412 7089 : dfFromGreenwich = poSRS->GetPrimeMeridian();
1413 :
1414 : /* -------------------------------------------------------------------- */
1415 : /* Set linear units. */
1416 : /* -------------------------------------------------------------------- */
1417 7089 : if( !EPSGGetUOMLengthInfo( nUOMLength, &pszUOMLengthName, &dfInMeters ) )
1418 : {
1419 0 : CPLFree(pszPCSName);
1420 0 : return OGRERR_UNSUPPORTED_SRS;
1421 : }
1422 :
1423 7089 : poSRS->SetLinearUnits( pszUOMLengthName, dfInMeters );
1424 7089 : poSRS->SetAuthority( "PROJCS|UNIT", "EPSG", nUOMLength );
1425 :
1426 7089 : CPLFree( pszUOMLengthName );
1427 7089 : CPLFree( pszPCSName );
1428 :
1429 : /* -------------------------------------------------------------------- */
1430 : /* Set projection and parameters. */
1431 : /* -------------------------------------------------------------------- */
1432 7089 : if( !EPSGGetProjTRFInfo( nPCSCode, &nProjMethod, anParmIds, adfProjParms ))
1433 0 : return OGRERR_UNSUPPORTED_SRS;
1434 :
1435 7089 : switch( nProjMethod )
1436 : {
1437 : case 9801:
1438 : case 9817: /* really LCC near conformal */
1439 : poSRS->SetLCC1SP( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1440 : OGR_FP( NatOriginScaleFactor ),
1441 134 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1442 134 : break;
1443 :
1444 : case 9802:
1445 : poSRS->SetLCC( OGR_FP( StdParallel1Lat ), OGR_FP( StdParallel2Lat ),
1446 : OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ),
1447 : OGR_FP( FalseOriginEasting ),
1448 1322 : OGR_FP( FalseOriginNorthing ));
1449 1322 : break;
1450 :
1451 : case 9803:
1452 : poSRS->SetLCCB( OGR_FP( StdParallel1Lat ), OGR_FP( StdParallel2Lat ),
1453 : OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ),
1454 : OGR_FP( FalseOriginEasting ),
1455 2 : OGR_FP( FalseOriginNorthing ));
1456 2 : break;
1457 :
1458 : case 9805:
1459 : poSRS->SetMercator2SP( OGR_FP( StdParallel1Lat ),
1460 : OGR_FP( NatOriginLat ), OGR_FP(NatOriginLong),
1461 4 : OGR_FP( FalseEasting ), OGR_FP(FalseNorthing) );
1462 :
1463 4 : break;
1464 :
1465 : case 9804:
1466 : case 9841: /* Mercator 1SP (Spherical) */
1467 : case 1024: /* Google Mercator */
1468 : poSRS->SetMercator( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1469 : OGR_FP( NatOriginScaleFactor ),
1470 30 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1471 :
1472 30 : if( nProjMethod == 1024 || nProjMethod == 9841 ) // override hack for google mercator.
1473 : {
1474 : poSRS->SetExtension( "PROJCS", "PROJ4",
1475 10 : "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs" );
1476 : }
1477 30 : break;
1478 :
1479 : case 9806:
1480 : poSRS->SetCS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1481 40 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1482 40 : break;
1483 :
1484 : case 9807:
1485 : poSRS->SetTM( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1486 : OGR_FP( NatOriginScaleFactor ),
1487 5211 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1488 5211 : break;
1489 :
1490 : case 9808:
1491 : poSRS->SetTMSO( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1492 : OGR_FP( NatOriginScaleFactor ),
1493 56 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1494 56 : break;
1495 :
1496 : case 9809:
1497 : poSRS->SetOS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1498 : OGR_FP( NatOriginScaleFactor ),
1499 43 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1500 43 : break;
1501 :
1502 : case 9810:
1503 : poSRS->SetPS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1504 : OGR_FP( NatOriginScaleFactor ),
1505 4 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1506 4 : break;
1507 :
1508 : case 9811:
1509 : poSRS->SetNZMG( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1510 2 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1511 2 : break;
1512 :
1513 : case 9812:
1514 : case 9813:
1515 : poSRS->SetHOM( OGR_FP( ProjCenterLat ), OGR_FP( ProjCenterLong ),
1516 : OGR_FP( Azimuth ),
1517 : OGR_FP( AngleRectifiedToSkewedGrid ),
1518 : OGR_FP( InitialLineScaleFactor ),
1519 24 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1520 :
1521 24 : poNode = poSRS->GetAttrNode( "PROJECTION" )->GetChild( 0 );
1522 24 : if( nProjMethod == 9813 )
1523 2 : poNode->SetValue( SRS_PT_LABORDE_OBLIQUE_MERCATOR );
1524 24 : break;
1525 :
1526 : case 9814:
1527 : /* NOTE: This is no longer used! Swiss Oblique Mercator gets
1528 : ** implemented using 9815 instead.
1529 : */
1530 : poSRS->SetSOC( OGR_FP( ProjCenterLat ), OGR_FP( ProjCenterLong ),
1531 0 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1532 0 : break;
1533 :
1534 : case 9815:
1535 : poSRS->SetHOM( OGR_FP( ProjCenterLat ), OGR_FP( ProjCenterLong ),
1536 : OGR_FP( Azimuth ),
1537 : OGR_FP( AngleRectifiedToSkewedGrid ),
1538 : OGR_FP( InitialLineScaleFactor ),
1539 : OGR_FP( ProjCenterEasting ),
1540 23 : OGR_FP( ProjCenterNorthing ) );
1541 23 : break;
1542 :
1543 : case 9816:
1544 : poSRS->SetTMG( OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ),
1545 : OGR_FP( FalseOriginEasting ),
1546 2 : OGR_FP( FalseOriginNorthing ) );
1547 2 : break;
1548 :
1549 : case 9818:
1550 : poSRS->SetPolyconic( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1551 4 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1552 4 : break;
1553 :
1554 : case 9819:
1555 : {
1556 2 : double dfCenterLong = OGR_FP( ProjCenterLong );
1557 :
1558 2 : if( dfCenterLong == 0.0 ) // See ticket #2559
1559 2 : dfCenterLong = OGR_FP( PolarLongOrigin );
1560 :
1561 : poSRS->SetKrovak( OGR_FP( ProjCenterLat ), dfCenterLong,
1562 : OGR_FP( Azimuth ),
1563 : OGR_FP( PseudoStdParallelLat ),
1564 : OGR_FP( PseudoStdParallelScaleFactor ),
1565 : OGR_FP( ProjCenterEasting ),
1566 2 : OGR_FP( ProjCenterNorthing ) );
1567 : }
1568 2 : break;
1569 :
1570 : case 9820:
1571 : poSRS->SetLAEA( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1572 14 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1573 14 : break;
1574 :
1575 : case 9821: /* this is the spherical form, and really needs different
1576 : equations which give different results but PROJ.4 doesn't
1577 : seem to support the spherical form. */
1578 : poSRS->SetLAEA( OGR_FP( SphericalOriginLat ),
1579 : OGR_FP( SphericalOriginLong ),
1580 0 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1581 0 : break;
1582 :
1583 : case 9822: /* Albers (Conic) Equal Area */
1584 : poSRS->SetACEA( OGR_FP( StdParallel1Lat ),
1585 : OGR_FP( StdParallel2Lat ),
1586 : OGR_FP( FalseOriginLat ),
1587 : OGR_FP( FalseOriginLong ),
1588 : OGR_FP( FalseOriginEasting ),
1589 40 : OGR_FP( FalseOriginNorthing ) );
1590 40 : break;
1591 :
1592 : case 9823: /* Equidistant Cylindrical / Plate Carre / Equirectangular */
1593 : poSRS->SetEquirectangular( OGR_FP( NatOriginLat ),
1594 : OGR_FP( NatOriginLong ),
1595 4 : 0.0, 0.0 );
1596 4 : break;
1597 :
1598 : case 9829: /* Polar Stereographic (Variant B) */
1599 : poSRS->SetPS( OGR_FP( PolarLatStdParallel ), OGR_FP(PolarLongOrigin),
1600 : 1.0,
1601 54 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1602 54 : break;
1603 :
1604 : case 9834: /* Lambert Cylindrical Equal Area (Spherical) bug #2659 */
1605 : poSRS->SetCEA( OGR_FP( StdParallel1Lat ), OGR_FP( NatOriginLong ),
1606 4 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1607 4 : break;
1608 :
1609 : default:
1610 : CPLDebug( "EPSG", "No WKT support for projection method %d.",
1611 70 : nProjMethod );
1612 70 : return OGRERR_UNSUPPORTED_SRS;
1613 : }
1614 :
1615 : /* -------------------------------------------------------------------- */
1616 : /* Set overall PCS authority code. */
1617 : /* -------------------------------------------------------------------- */
1618 7019 : poSRS->SetAuthority( "PROJCS", "EPSG", nPCSCode );
1619 :
1620 : /* -------------------------------------------------------------------- */
1621 : /* Set axes */
1622 : /* -------------------------------------------------------------------- */
1623 7019 : if( nCSC > 0 )
1624 : {
1625 7019 : SetEPSGAxisInfo( poSRS, "PROJCS", nCSC );
1626 7019 : CPLErrorReset();
1627 : }
1628 :
1629 7019 : return OGRERR_NONE;
1630 : }
1631 :
1632 : /************************************************************************/
1633 : /* importFromEPSG() */
1634 : /************************************************************************/
1635 :
1636 : /**
1637 : * \brief Initialize SRS based on EPSG GCS or PCS code.
1638 : *
1639 : * This method will initialize the spatial reference based on the
1640 : * passed in EPSG GCS or PCS code. The coordinate system definitions
1641 : * are normally read from the EPSG derived support files such as
1642 : * pcs.csv, gcs.csv, pcs.override.csv, gcs.override.csv and falling
1643 : * back to search for a PROJ.4 epsg init file or a definition in epsg.wkt.
1644 : *
1645 : * These support files are normally searched for in /usr/local/share/gdal
1646 : * or in the directory identified by the GDAL_DATA configuration option.
1647 : * See CPLFindFile() for details.
1648 : *
1649 : * This method is relatively expensive, and generally involves quite a bit
1650 : * of text file scanning. Reasonable efforts should be made to avoid calling
1651 : * it many times for the same coordinate system.
1652 : *
1653 : * This method is similar to importFromEPSGA() except that EPSG preferred
1654 : * axis ordering will *not* be applied for geographic coordinate systems.
1655 : * EPSG normally defines geographic coordinate systems to use lat/long
1656 : * contrary to typical GIS use).
1657 : *
1658 : * This method is the same as the C function OSRImportFromEPSG().
1659 : *
1660 : * @param nCode a GCS or PCS code from the horizontal coordinate system table.
1661 : *
1662 : * @return OGRERR_NONE on success, or an error code on failure.
1663 : */
1664 :
1665 8119 : OGRErr OGRSpatialReference::importFromEPSG( int nCode )
1666 :
1667 : {
1668 8119 : OGRErr eErr = importFromEPSGA( nCode );
1669 :
1670 : // Strip any GCS axis settings found.
1671 8119 : if( eErr == OGRERR_NONE )
1672 : {
1673 8049 : OGR_SRSNode *poGEOGCS = GetAttrNode( "GEOGCS" );
1674 :
1675 8049 : if( poGEOGCS != NULL )
1676 8049 : poGEOGCS->StripNodes( "AXIS" );
1677 : }
1678 :
1679 8119 : return eErr;
1680 : }
1681 :
1682 : /************************************************************************/
1683 : /* OSRImportFromEPSG() */
1684 : /************************************************************************/
1685 :
1686 : /**
1687 : * \brief Initialize SRS based on EPSG GCS or PCS code.
1688 : *
1689 : * This function is the same as OGRSpatialReference::importFromEPSG().
1690 : */
1691 :
1692 55 : OGRErr CPL_STDCALL OSRImportFromEPSG( OGRSpatialReferenceH hSRS, int nCode )
1693 :
1694 : {
1695 55 : VALIDATE_POINTER1( hSRS, "OSRImportFromEPSG", CE_Failure );
1696 :
1697 55 : return ((OGRSpatialReference *) hSRS)->importFromEPSG( nCode );
1698 : }
1699 :
1700 : /************************************************************************/
1701 : /* importFromEPSGA() */
1702 : /************************************************************************/
1703 :
1704 : /**
1705 : * \brief Initialize SRS based on EPSG GCS or PCS code.
1706 : *
1707 : * This method will initialize the spatial reference based on the
1708 : * passed in EPSG GCS or PCS code.
1709 : *
1710 : * This method is similar to importFromEPSG() except that EPSG preferred
1711 : * axis ordering *will* be applied for geographic coordinate systems.
1712 : * EPSG normally defines geographic coordinate systems to use lat/long
1713 : * contrary to typical GIS use). See OGRSpatialReference::importFromEPSG()
1714 : * for more details on operation of this method.
1715 : *
1716 : * This method is the same as the C function OSRImportFromEPSGA().
1717 : *
1718 : * @param nCode a GCS or PCS code from the horizontal coordinate system table.
1719 : *
1720 : * @return OGRERR_NONE on success, or an error code on failure.
1721 : */
1722 :
1723 8131 : OGRErr OGRSpatialReference::importFromEPSGA( int nCode )
1724 :
1725 : {
1726 : OGRErr eErr;
1727 :
1728 8131 : bNormInfoSet = FALSE;
1729 :
1730 : /* -------------------------------------------------------------------- */
1731 : /* Clear any existing definition. */
1732 : /* -------------------------------------------------------------------- */
1733 8131 : if( GetRoot() != NULL )
1734 : {
1735 7891 : delete poRoot;
1736 7891 : poRoot = NULL;
1737 : }
1738 :
1739 : /* -------------------------------------------------------------------- */
1740 : /* Verify that we can find the required filename(s). */
1741 : /* -------------------------------------------------------------------- */
1742 8131 : if( CSVScanFileByName( CSVFilename( "gcs.csv" ),
1743 : "COORD_REF_SYS_CODE",
1744 : "4269", CC_Integer ) == NULL )
1745 : {
1746 : CPLError( CE_Failure, CPLE_OpenFailed,
1747 : "Unable to open EPSG support file %s.\n"
1748 : "Try setting the GDAL_DATA environment variable to point to the\n"
1749 : "directory containing EPSG csv files.",
1750 0 : CSVFilename( "gcs.csv" ) );
1751 0 : return OGRERR_FAILURE;
1752 : }
1753 :
1754 : /* -------------------------------------------------------------------- */
1755 : /* Is this a GeogCS code? this is inadequate as a criteria */
1756 : /* -------------------------------------------------------------------- */
1757 8131 : if( EPSGGetGCSInfo( nCode, NULL, NULL, NULL, NULL, NULL, NULL, NULL ) )
1758 1036 : eErr = SetEPSGGeogCS( this, nCode );
1759 : else
1760 7095 : eErr = SetEPSGProjCS( this, nCode );
1761 :
1762 : /* -------------------------------------------------------------------- */
1763 : /* If we get it as an unsupported code, try looking it up in */
1764 : /* the epsg.wkt coordinate system dictionary. */
1765 : /* -------------------------------------------------------------------- */
1766 8131 : if( eErr == OGRERR_UNSUPPORTED_SRS )
1767 : {
1768 : char szCode[32];
1769 76 : sprintf( szCode, "%d", nCode );
1770 76 : eErr = importFromDict( "epsg.wkt", szCode );
1771 : }
1772 :
1773 : /* -------------------------------------------------------------------- */
1774 : /* If we get it as an unsupported code, try looking it up in */
1775 : /* the PROJ.4 support file(s). */
1776 : /* -------------------------------------------------------------------- */
1777 8131 : if( eErr == OGRERR_UNSUPPORTED_SRS )
1778 : {
1779 : char szWrkDefn[100];
1780 : char *pszNormalized;
1781 :
1782 76 : sprintf( szWrkDefn, "+init=epsg:%d", nCode );
1783 :
1784 76 : pszNormalized = OCTProj4Normalize( szWrkDefn );
1785 :
1786 76 : if( strstr(pszNormalized,"proj=") != NULL )
1787 6 : eErr = importFromProj4( pszNormalized );
1788 :
1789 76 : CPLFree( pszNormalized );
1790 : }
1791 :
1792 : /* -------------------------------------------------------------------- */
1793 : /* Push in authority information if we were successful, and it */
1794 : /* is not already present. */
1795 : /* -------------------------------------------------------------------- */
1796 : const char *pszAuthName;
1797 :
1798 8131 : if( IsProjected() )
1799 7095 : pszAuthName = GetAuthorityName( "PROJCS" );
1800 : else
1801 1036 : pszAuthName = GetAuthorityName( "GEOGCS" );
1802 :
1803 :
1804 8131 : if( eErr == OGRERR_NONE && pszAuthName == NULL )
1805 : {
1806 6 : if( IsProjected() )
1807 6 : SetAuthority( "PROJCS", "EPSG", nCode );
1808 0 : else if( IsGeographic() )
1809 0 : SetAuthority( "GEOGCS", "EPSG", nCode );
1810 :
1811 6 : eErr = FixupOrdering();
1812 : }
1813 :
1814 : /* -------------------------------------------------------------------- */
1815 : /* Otherwise officially issue an error message. */
1816 : /* -------------------------------------------------------------------- */
1817 8131 : if( eErr == OGRERR_UNSUPPORTED_SRS )
1818 : {
1819 : CPLError( CE_Failure, CPLE_NotSupported,
1820 : "EPSG PCS/GCS code %d not found in EPSG support files. Is this a valid\nEPSG coordinate system?",
1821 70 : nCode );
1822 : }
1823 :
1824 8131 : return eErr;
1825 : }
1826 :
1827 : /************************************************************************/
1828 : /* OSRImportFromEPSGA() */
1829 : /************************************************************************/
1830 :
1831 : /**
1832 : * \brief Initialize SRS based on EPSG GCS or PCS code.
1833 : *
1834 : * This function is the same as OGRSpatialReference::importFromEPSGA().
1835 : */
1836 :
1837 2 : OGRErr CPL_STDCALL OSRImportFromEPSGA( OGRSpatialReferenceH hSRS, int nCode )
1838 :
1839 : {
1840 2 : VALIDATE_POINTER1( hSRS, "OSRImportFromEPSGA", CE_Failure );
1841 :
1842 2 : return ((OGRSpatialReference *) hSRS)->importFromEPSGA( nCode );
1843 : }
1844 :
1845 : /************************************************************************/
1846 : /* SetStatePlane() */
1847 : /************************************************************************/
1848 :
1849 : /**
1850 : * \brief Set State Plane projection definition.
1851 : *
1852 : * This will attempt to generate a complete definition of a state plane
1853 : * zone based on generating the entire SRS from the EPSG tables. If the
1854 : * EPSG tables are unavailable, it will produce a stubbed LOCAL_CS definition
1855 : * and return OGRERR_FAILURE.
1856 : *
1857 : * This method is the same as the C function OSRSetStatePlaneWithUnits().
1858 : *
1859 : * @param nZone State plane zone number, in the USGS numbering scheme (as
1860 : * dinstinct from the Arc/Info and Erdas numbering scheme.
1861 : *
1862 : * @param bNAD83 TRUE if the NAD83 zone definition should be used or FALSE
1863 : * if the NAD27 zone definition should be used.
1864 : *
1865 : * @param pszOverrideUnitName Linear unit name to apply overriding the
1866 : * legal definition for this zone.
1867 : *
1868 : * @param dfOverrideUnit Linear unit conversion factor to apply overriding
1869 : * the legal definition for this zone.
1870 : *
1871 : * @return OGRERR_NONE on success, or OGRERR_FAILURE on failure, mostly likely
1872 : * due to the EPSG tables not being accessable.
1873 : */
1874 :
1875 : OGRErr OGRSpatialReference::SetStatePlane( int nZone, int bNAD83,
1876 : const char *pszOverrideUnitName,
1877 10 : double dfOverrideUnit )
1878 :
1879 : {
1880 : int nAdjustedId;
1881 : int nPCSCode;
1882 : char szID[32];
1883 :
1884 : /* -------------------------------------------------------------------- */
1885 : /* Get the index id from stateplane.csv. */
1886 : /* -------------------------------------------------------------------- */
1887 10 : if( bNAD83 )
1888 9 : nAdjustedId = nZone;
1889 : else
1890 1 : nAdjustedId = nZone + 10000;
1891 :
1892 : /* -------------------------------------------------------------------- */
1893 : /* Turn this into a PCS code. We assume there will only be one */
1894 : /* PCS corresponding to each Proj_ code since the proj code */
1895 : /* already effectively indicates NAD27 or NAD83. */
1896 : /* -------------------------------------------------------------------- */
1897 10 : sprintf( szID, "%d", nAdjustedId );
1898 : nPCSCode =
1899 : atoi( CSVGetField( CSVFilename( "stateplane.csv" ),
1900 : "ID", szID, CC_Integer,
1901 10 : "EPSG_PCS_CODE" ) );
1902 10 : if( nPCSCode < 1 )
1903 : {
1904 : char szName[128];
1905 : static int bFailureReported = FALSE;
1906 :
1907 0 : if( !bFailureReported )
1908 : {
1909 0 : bFailureReported = TRUE;
1910 : CPLError( CE_Warning, CPLE_OpenFailed,
1911 : "Unable to find state plane zone in stateplane.csv,\n"
1912 : "likely because the GDAL data files cannot be found. Using\n"
1913 0 : "incomplete definition of state plane zone.\n" );
1914 : }
1915 :
1916 0 : Clear();
1917 0 : if( bNAD83 )
1918 : {
1919 0 : sprintf( szName, "State Plane Zone %d / NAD83", nZone );
1920 0 : SetLocalCS( szName );
1921 0 : SetLinearUnits( SRS_UL_METER, 1.0 );
1922 : }
1923 : else
1924 : {
1925 0 : sprintf( szName, "State Plane Zone %d / NAD27", nZone );
1926 0 : SetLocalCS( szName );
1927 0 : SetLinearUnits( SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV) );
1928 : }
1929 :
1930 0 : return OGRERR_FAILURE;
1931 : }
1932 :
1933 : /* -------------------------------------------------------------------- */
1934 : /* Define based on a full EPSG definition of the zone. */
1935 : /* -------------------------------------------------------------------- */
1936 10 : OGRErr eErr = importFromEPSG( nPCSCode );
1937 :
1938 10 : if( eErr != OGRERR_NONE )
1939 0 : return eErr;
1940 :
1941 : /* -------------------------------------------------------------------- */
1942 : /* Apply units override if required. */
1943 : /* */
1944 : /* We will need to adjust the linear projection parameter to */
1945 : /* match the provided units, and clear the authority code. */
1946 : /* -------------------------------------------------------------------- */
1947 10 : if( dfOverrideUnit != 0.0
1948 : && fabs(dfOverrideUnit - GetLinearUnits()) > 0.0000000001 )
1949 : {
1950 4 : double dfFalseEasting = GetNormProjParm( SRS_PP_FALSE_EASTING );
1951 4 : double dfFalseNorthing= GetNormProjParm( SRS_PP_FALSE_NORTHING);
1952 : OGR_SRSNode *poPROJCS;
1953 :
1954 4 : SetLinearUnits( pszOverrideUnitName, dfOverrideUnit );
1955 :
1956 4 : SetNormProjParm( SRS_PP_FALSE_EASTING, dfFalseEasting );
1957 4 : SetNormProjParm( SRS_PP_FALSE_NORTHING, dfFalseNorthing );
1958 :
1959 4 : poPROJCS = GetAttrNode( "PROJCS" );
1960 4 : if( poPROJCS != NULL && poPROJCS->FindChild( "AUTHORITY" ) != -1 )
1961 : {
1962 4 : poPROJCS->DestroyChild( poPROJCS->FindChild( "AUTHORITY" ) );
1963 : }
1964 : }
1965 :
1966 10 : return OGRERR_NONE;
1967 : }
1968 :
1969 : /************************************************************************/
1970 : /* OSRSetStatePlane() */
1971 : /************************************************************************/
1972 :
1973 0 : OGRErr OSRSetStatePlane( OGRSpatialReferenceH hSRS, int nZone, int bNAD83 )
1974 :
1975 : {
1976 0 : VALIDATE_POINTER1( hSRS, "OSRSetStatePlane", CE_Failure );
1977 :
1978 0 : return ((OGRSpatialReference *) hSRS)->SetStatePlane( nZone, bNAD83 );
1979 : }
1980 :
1981 : /************************************************************************/
1982 : /* OSRSetStatePlaneWithUnits() */
1983 : /************************************************************************/
1984 :
1985 : OGRErr OSRSetStatePlaneWithUnits( OGRSpatialReferenceH hSRS,
1986 : int nZone, int bNAD83,
1987 : const char *pszOverrideUnitName,
1988 2 : double dfOverrideUnit )
1989 :
1990 : {
1991 2 : VALIDATE_POINTER1( hSRS, "OSRSetStatePlaneWithUnits", CE_Failure );
1992 :
1993 : return ((OGRSpatialReference *) hSRS)->SetStatePlane( nZone, bNAD83,
1994 : pszOverrideUnitName,
1995 2 : dfOverrideUnit );
1996 : }
1997 :
1998 : /************************************************************************/
1999 : /* GetEPSGGeogCS() */
2000 : /* */
2001 : /* Try to establish what the EPSG code for this coordinate */
2002 : /* systems GEOGCS might be. Returns -1 if no reasonable guess */
2003 : /* can be made. */
2004 : /* */
2005 : /* TODO: We really need to do some name lookups. */
2006 : /************************************************************************/
2007 :
2008 106 : int OGRSpatialReference::GetEPSGGeogCS()
2009 :
2010 : {
2011 106 : const char *pszAuthName = GetAuthorityName( "GEOGCS" );
2012 :
2013 : /* -------------------------------------------------------------------- */
2014 : /* Do we already have it? */
2015 : /* -------------------------------------------------------------------- */
2016 106 : if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") )
2017 45 : return atoi(GetAuthorityCode( "GEOGCS" ));
2018 :
2019 : /* -------------------------------------------------------------------- */
2020 : /* Get the datum and geogcs names. */
2021 : /* -------------------------------------------------------------------- */
2022 61 : const char *pszGEOGCS = GetAttrValue( "GEOGCS" );
2023 61 : const char *pszDatum = GetAttrValue( "DATUM" );
2024 :
2025 : // We can only operate on coordinate systems with a geogcs.
2026 61 : if( pszGEOGCS == NULL || pszDatum == NULL )
2027 0 : return -1;
2028 :
2029 : /* -------------------------------------------------------------------- */
2030 : /* Is this a "well known" geographic coordinate system? */
2031 : /* -------------------------------------------------------------------- */
2032 : int bWGS, bNAD;
2033 :
2034 : bWGS = strstr(pszGEOGCS,"WGS") != NULL
2035 : || strstr(pszDatum, "WGS")
2036 : || strstr(pszGEOGCS,"World Geodetic System")
2037 : || strstr(pszGEOGCS,"World_Geodetic_System")
2038 : || strstr(pszDatum, "World Geodetic System")
2039 61 : || strstr(pszDatum, "World_Geodetic_System");
2040 :
2041 : bNAD = strstr(pszGEOGCS,"NAD") != NULL
2042 : || strstr(pszDatum, "NAD")
2043 : || strstr(pszGEOGCS,"North American")
2044 : || strstr(pszGEOGCS,"North_American")
2045 : || strstr(pszDatum, "North American")
2046 61 : || strstr(pszDatum, "North_American");
2047 :
2048 61 : if( bWGS && (strstr(pszGEOGCS,"84") || strstr(pszDatum,"84")) )
2049 44 : return 4326;
2050 :
2051 17 : if( bWGS && (strstr(pszGEOGCS,"72") || strstr(pszDatum,"72")) )
2052 0 : return 4322;
2053 :
2054 17 : if( bNAD && (strstr(pszGEOGCS,"83") || strstr(pszDatum,"83")) )
2055 2 : return 4269;
2056 :
2057 15 : if( bNAD && (strstr(pszGEOGCS,"27") || strstr(pszDatum,"27")) )
2058 1 : return 4267;
2059 :
2060 : /* -------------------------------------------------------------------- */
2061 : /* If we know the datum, associate the most likely GCS with */
2062 : /* it. */
2063 : /* -------------------------------------------------------------------- */
2064 14 : pszAuthName = GetAuthorityName( "GEOGCS|DATUM" );
2065 :
2066 14 : if( pszAuthName != NULL
2067 : && EQUAL(pszAuthName,"epsg")
2068 : && GetPrimeMeridian() == 0.0 )
2069 : {
2070 0 : int nDatum = atoi(GetAuthorityCode("GEOGCS|DATUM"));
2071 :
2072 0 : if( nDatum >= 6000 && nDatum <= 6999 )
2073 0 : return nDatum - 2000;
2074 : }
2075 :
2076 14 : return -1;
2077 : }
2078 :
2079 : /************************************************************************/
2080 : /* AutoIdentifyEPSG() */
2081 : /************************************************************************/
2082 :
2083 : /**
2084 : * \brief Set EPSG authority info if possible.
2085 : *
2086 : * This method inspects a WKT definition, and adds EPSG authority nodes
2087 : * where an aspect of the coordinate system can be easily and safely
2088 : * corresponded with an EPSG identifier. In practice, this method will
2089 : * evolve over time. In theory it can add authority nodes for any object
2090 : * (ie. spheroid, datum, GEOGCS, units, and PROJCS) that could have an
2091 : * authority node. Mostly this is useful to inserting appropriate
2092 : * PROJCS codes for common formulations (like UTM n WGS84).
2093 : *
2094 : * If it success the OGRSpatialReference is updated in place, and the
2095 : * method return OGRERR_NONE. If the method fails to identify the
2096 : * general coordinate system OGRERR_UNSUPPORTED_SRS is returned but no
2097 : * error message is posted via CPLError().
2098 : *
2099 : * This method is the same as the C function OSRAutoIdentifyEPSG().
2100 : *
2101 : * @return OGRERR_NONE or OGRERR_UNSUPPORTED_SRS.
2102 : */
2103 :
2104 44 : OGRErr OGRSpatialReference::AutoIdentifyEPSG()
2105 :
2106 : {
2107 : /* -------------------------------------------------------------------- */
2108 : /* Do we have a GEOGCS node, but no authority? If so, try */
2109 : /* guessing it. */
2110 : /* -------------------------------------------------------------------- */
2111 44 : if( (IsProjected() || IsGeographic())
2112 : && GetAuthorityCode( "GEOGCS" ) == NULL )
2113 : {
2114 18 : int nGCS = GetEPSGGeogCS();
2115 18 : if( nGCS != -1 )
2116 4 : SetAuthority( "GEOGCS", "EPSG", nGCS );
2117 : }
2118 :
2119 : /* -------------------------------------------------------------------- */
2120 : /* Is this a UTM coordinate system with a common GEOGCS? */
2121 : /* -------------------------------------------------------------------- */
2122 : int nZone, bNorth;
2123 44 : if( (nZone = GetUTMZone( &bNorth )) != 0
2124 : && GetAuthorityCode( "PROJCS") == NULL )
2125 : {
2126 : const char *pszAuthName, *pszAuthCode;
2127 :
2128 34 : pszAuthName = GetAuthorityName( "PROJCS|GEOGCS" );
2129 34 : pszAuthCode = GetAuthorityCode( "PROJCS|GEOGCS" );
2130 :
2131 34 : if( pszAuthName == NULL || pszAuthCode == NULL )
2132 : {
2133 : /* don't exactly recognise datum */
2134 : }
2135 23 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4326 )
2136 : { // WGS84
2137 0 : if( bNorth )
2138 0 : SetAuthority( "PROJCS", "EPSG", 32600 + nZone );
2139 : else
2140 0 : SetAuthority( "PROJCS", "EPSG", 32700 + nZone );
2141 : }
2142 46 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4267
2143 : && nZone >= 3 && nZone <= 22 && bNorth )
2144 23 : SetAuthority( "PROJCS", "EPSG", 26700 + nZone ); // NAD27
2145 0 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4269
2146 : && nZone >= 3 && nZone <= 23 && bNorth )
2147 0 : SetAuthority( "PROJCS", "EPSG", 26900 + nZone ); // NAD83
2148 0 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4322 )
2149 : { // WGS72
2150 0 : if( bNorth )
2151 0 : SetAuthority( "PROJCS", "EPSG", 32200 + nZone );
2152 : else
2153 0 : SetAuthority( "PROJCS", "EPSG", 32300 + nZone );
2154 : }
2155 : }
2156 :
2157 : /* -------------------------------------------------------------------- */
2158 : /* Return. */
2159 : /* -------------------------------------------------------------------- */
2160 44 : if( IsProjected() && GetAuthorityCode("PROJCS") != NULL )
2161 23 : return OGRERR_NONE;
2162 21 : else if( IsGeographic() && GetAuthorityCode("GEOGCS") != NULL )
2163 2 : return OGRERR_NONE;
2164 : else
2165 19 : return OGRERR_UNSUPPORTED_SRS;
2166 : }
2167 :
2168 : /************************************************************************/
2169 : /* OSRAutoIdentifyEPSG() */
2170 : /************************************************************************/
2171 :
2172 0 : OGRErr OSRAutoIdentifyEPSG( OGRSpatialReferenceH hSRS )
2173 :
2174 : {
2175 0 : VALIDATE_POINTER1( hSRS, "OSRAutoIdentifyEPSG", CE_Failure );
2176 :
2177 0 : return ((OGRSpatialReference *) hSRS)->AutoIdentifyEPSG();
2178 : }
2179 :
2180 : /************************************************************************/
2181 : /* EPSGTreatsAsLatLong() */
2182 : /************************************************************************/
2183 :
2184 : /**
2185 : * \brief This method returns TRUE if EPSG feels this geographic coordinate
2186 : * system should be treated as having lat/long coordinate ordering.
2187 : *
2188 : * Currently this returns TRUE for all geographic coordinate systems
2189 : * with an EPSG code set, and AXIS values set defining it as lat, long.
2190 : * Note that coordinate systems with an EPSG code and no axis settings
2191 : * will be assumed to not be lat/long.
2192 : *
2193 : * FALSE will be returned for all coordinate systems that are not geographic,
2194 : * or that do not have an EPSG code set.
2195 : *
2196 : * @return TRUE or FALSE.
2197 : */
2198 :
2199 2 : int OGRSpatialReference::EPSGTreatsAsLatLong()
2200 :
2201 : {
2202 2 : if( !IsGeographic() )
2203 0 : return FALSE;
2204 :
2205 2 : const char *pszAuth = GetAuthorityName( "GEOGCS" );
2206 :
2207 2 : if( pszAuth == NULL || !EQUAL(pszAuth,"EPSG") )
2208 0 : return FALSE;
2209 :
2210 2 : OGR_SRSNode *poFirstAxis = GetAttrNode( "GEOGCS|AXIS" );
2211 :
2212 2 : if( poFirstAxis == NULL )
2213 0 : return TRUE;
2214 :
2215 2 : if( poFirstAxis->GetChildCount() >= 2
2216 : && EQUAL(poFirstAxis->GetChild(1)->GetValue(),"NORTH") )
2217 2 : return TRUE;
2218 :
2219 0 : return FALSE;
2220 : }
2221 :
2222 : /************************************************************************/
2223 : /* OSREPSGTreatsAsLatLong() */
2224 : /************************************************************************/
2225 :
2226 2 : int OSREPSGTreatsAsLatLong( OGRSpatialReferenceH hSRS )
2227 :
2228 : {
2229 2 : VALIDATE_POINTER1( hSRS, "OSREPSGTreatsAsLatLong", CE_Failure );
2230 :
2231 2 : return ((OGRSpatialReference *) hSRS)->EPSGTreatsAsLatLong();
2232 : }
2233 :
|