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