1 : /******************************************************************************
2 : * $Id: ogr_fromepsg.cpp 24323 2012-04-27 05:45:18Z warmerdam $
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 24323 2012-04-27 05:45:18Z warmerdam $");
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 137800 : void OGREPSGDatumNameMassage( char ** ppszDatum )
66 :
67 : {
68 : int i, j;
69 137800 : char *pszDatum = *ppszDatum;
70 :
71 137800 : if (pszDatum[0] == '\0')
72 0 : return;
73 :
74 : /* -------------------------------------------------------------------- */
75 : /* Translate non-alphanumeric values to underscores. */
76 : /* -------------------------------------------------------------------- */
77 3130025 : for( i = 0; pszDatum[i] != '\0'; i++ )
78 : {
79 13875505 : if( pszDatum[i] != '+'
80 5223496 : && !(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z')
81 4497820 : && !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z')
82 1161964 : && !(pszDatum[i] >= '0' && pszDatum[i] <= '9') )
83 : {
84 359700 : pszDatum[i] = '_';
85 : }
86 : }
87 :
88 : /* -------------------------------------------------------------------- */
89 : /* Remove repeated and trailing underscores. */
90 : /* -------------------------------------------------------------------- */
91 2992225 : for( i = 1, j = 0; pszDatum[i] != '\0'; i++ )
92 : {
93 2854425 : if( pszDatum[j] == '_' && pszDatum[i] == '_' )
94 25072 : continue;
95 :
96 2829353 : pszDatum[++j] = pszDatum[i];
97 : }
98 137800 : if( pszDatum[j] == '_' )
99 20028 : pszDatum[j] = '\0';
100 : else
101 117772 : pszDatum[j+1] = '\0';
102 :
103 : /* -------------------------------------------------------------------- */
104 : /* Search for datum equivelences. Specific massaged names get */
105 : /* mapped to OpenGIS specified names. */
106 : /* -------------------------------------------------------------------- */
107 806246 : for( i = 0; papszDatumEquiv[i] != NULL; i += 2 )
108 : {
109 674339 : if( EQUAL(*ppszDatum,papszDatumEquiv[i]) )
110 : {
111 5893 : CPLFree( *ppszDatum );
112 5893 : *ppszDatum = CPLStrdup( papszDatumEquiv[i+1] );
113 5893 : 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 56640 : EPSGAngleStringToDD( const char * pszAngle, int nUOMAngle )
126 :
127 : {
128 : double dfAngle;
129 :
130 56640 : if( nUOMAngle == 9110 ) /* DDD.MMSSsss */
131 : {
132 : char *pszDecimal;
133 :
134 25546 : dfAngle = ABS(atoi(pszAngle));
135 25546 : pszDecimal = (char *) strchr(pszAngle,'.');
136 25546 : if( pszDecimal != NULL && strlen(pszDecimal) > 1 )
137 : {
138 : char szMinutes[3];
139 : char szSeconds[64];
140 :
141 20238 : szMinutes[0] = pszDecimal[1];
142 29854 : if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' )
143 9616 : szMinutes[1] = pszDecimal[2];
144 : else
145 10622 : szMinutes[1] = '0';
146 :
147 20238 : szMinutes[2] = '\0';
148 20238 : dfAngle += atoi(szMinutes) / 60.0;
149 :
150 20238 : if( strlen(pszDecimal) > 3 )
151 : {
152 2654 : szSeconds[0] = pszDecimal[3];
153 5142 : if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' )
154 : {
155 2488 : szSeconds[1] = pszDecimal[4];
156 2488 : szSeconds[2] = '.';
157 2488 : strncpy( szSeconds+3, pszDecimal + 5, sizeof(szSeconds)-3 );
158 2488 : szSeconds[sizeof(szSeconds)-1] = 0;
159 : }
160 : else
161 : {
162 166 : szSeconds[1] = '0';
163 166 : szSeconds[2] = '\0';
164 : }
165 2654 : dfAngle += CPLAtof(szSeconds) / 3600.0;
166 : }
167 : }
168 :
169 25546 : if( pszAngle[0] == '-' )
170 7984 : dfAngle *= -1;
171 : }
172 31826 : else if( nUOMAngle == 9105 || nUOMAngle == 9106 ) /* grad */
173 : {
174 732 : dfAngle = 180 * (CPLAtof(pszAngle ) / 200);
175 : }
176 30362 : else if( nUOMAngle == 9101 ) /* radians */
177 : {
178 0 : dfAngle = 180 * (CPLAtof(pszAngle ) / PI);
179 : }
180 30362 : else if( nUOMAngle == 9103 ) /* arc-minute */
181 : {
182 0 : dfAngle = CPLAtof(pszAngle) / 60;
183 : }
184 30362 : 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 30362 : CPLAssert( nUOMAngle == 9102 || nUOMAngle == 0 );
191 :
192 30362 : dfAngle = CPLAtof(pszAngle );
193 : }
194 :
195 56640 : return( dfAngle );
196 : }
197 :
198 : /************************************************************************/
199 : /* EPSGGetUOMAngleInfo() */
200 : /************************************************************************/
201 :
202 36204 : int EPSGGetUOMAngleInfo( int nUOMAngleCode,
203 : char **ppszUOMName,
204 : double * pdfInDegrees )
205 :
206 : {
207 36204 : const char *pszUOMName = NULL;
208 36204 : 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 36204 : if( nUOMAngleCode == 9102 || nUOMAngleCode == 9107
216 : || nUOMAngleCode == 9108 || nUOMAngleCode == 9110
217 : || nUOMAngleCode == 9122 )
218 : {
219 35926 : if( ppszUOMName != NULL )
220 35926 : *ppszUOMName = CPLStrdup("degree");
221 35926 : if( pdfInDegrees != NULL )
222 35926 : *pdfInDegrees = 1.0;
223 35926 : return TRUE;
224 : }
225 :
226 278 : pszFilename = CSVFilename( "unit_of_measure.csv" );
227 :
228 278 : sprintf( szSearchKey, "%d", nUOMAngleCode );
229 : pszUOMName = CSVGetField( pszFilename,
230 : "UOM_CODE", szSearchKey, CC_Integer,
231 278 : "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 278 : if( pszUOMName != NULL )
240 : {
241 : double dfFactorB, dfFactorC;
242 :
243 : dfFactorB =
244 : CPLAtof(CSVGetField( pszFilename,
245 : "UOM_CODE", szSearchKey, CC_Integer,
246 278 : "FACTOR_B" ));
247 :
248 : dfFactorC =
249 : CPLAtof(CSVGetField( pszFilename,
250 : "UOM_CODE", szSearchKey, CC_Integer,
251 278 : "FACTOR_C" ));
252 :
253 278 : if( dfFactorC != 0.0 )
254 278 : 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 278 : if( nUOMAngleCode == 9105 )
259 278 : 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 278 : if( ppszUOMName != NULL )
317 : {
318 278 : if( pszUOMName != NULL )
319 278 : *ppszUOMName = CPLStrdup( pszUOMName );
320 : else
321 0 : *ppszUOMName = NULL;
322 : }
323 :
324 278 : if( pdfInDegrees != NULL )
325 278 : *pdfInDegrees = dfInDegrees;
326 :
327 278 : 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 106530 : 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 106530 : if( nUOMLengthCode == 9001 )
353 : {
354 95996 : if( ppszUOMName != NULL )
355 20128 : *ppszUOMName = CPLStrdup( "metre" );
356 95996 : if( pdfInMeters != NULL )
357 95996 : *pdfInMeters = 1.0;
358 :
359 95996 : return TRUE;
360 : }
361 :
362 : /* -------------------------------------------------------------------- */
363 : /* Search the units database for this unit. If we don't find */
364 : /* it return failure. */
365 : /* -------------------------------------------------------------------- */
366 10534 : sprintf( szSearchKey, "%d", nUOMLengthCode );
367 : papszUnitsRecord =
368 10534 : CSVScanFileByName( UOM_FILENAME, "UOM_CODE", szSearchKey, CC_Integer );
369 :
370 10534 : if( papszUnitsRecord == NULL )
371 0 : return FALSE;
372 :
373 : /* -------------------------------------------------------------------- */
374 : /* Get the name, if requested. */
375 : /* -------------------------------------------------------------------- */
376 10534 : if( ppszUOMName != NULL )
377 : {
378 3302 : iNameField = CSVGetFileFieldId( UOM_FILENAME, "UNIT_OF_MEAS_NAME" );
379 3302 : *ppszUOMName = CPLStrdup( CSLGetField(papszUnitsRecord, iNameField) );
380 : }
381 :
382 : /* -------------------------------------------------------------------- */
383 : /* Get the A and B factor fields, and create the multiplicative */
384 : /* factor. */
385 : /* -------------------------------------------------------------------- */
386 10534 : if( pdfInMeters != NULL )
387 : {
388 : int iBFactorField, iCFactorField;
389 :
390 10534 : iBFactorField = CSVGetFileFieldId( UOM_FILENAME, "FACTOR_B" );
391 10534 : iCFactorField = CSVGetFileFieldId( UOM_FILENAME, "FACTOR_C" );
392 :
393 10534 : if( CPLAtof(CSLGetField(papszUnitsRecord, iCFactorField)) > 0.0 )
394 : *pdfInMeters = CPLAtof(CSLGetField(papszUnitsRecord,iBFactorField))
395 10534 : / CPLAtof(CSLGetField(papszUnitsRecord, iCFactorField));
396 : else
397 0 : *pdfInMeters = 0.0;
398 : }
399 :
400 10534 : 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 36218 : 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 36218 : pszFilename = CSVFilename("gcs.override.csv");
428 36218 : sprintf( szCode, "%d", nGeogCS );
429 : papszLine = CSVScanFileByName( pszFilename,
430 : "COORD_REF_SYS_CODE",
431 36218 : szCode, CC_Integer );
432 36218 : if( papszLine == NULL )
433 : {
434 36218 : pszFilename = CSVFilename("gcs.csv");
435 36218 : sprintf( szCode, "%d", nGeogCS );
436 : papszLine = CSVScanFileByName( pszFilename,
437 : "COORD_REF_SYS_CODE",
438 36218 : szCode, CC_Integer );
439 : }
440 :
441 36218 : 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 36218 : "COORD_OP_METHOD_CODE")));
451 36218 : if( nMethodCode != 9603 && nMethodCode != 9607 && nMethodCode != 9606 )
452 8742 : return FALSE;
453 :
454 : /* -------------------------------------------------------------------- */
455 : /* Fetch the transformation parameters. */
456 : /* -------------------------------------------------------------------- */
457 27476 : iDXField = CSVGetFileFieldId(pszFilename, "DX");
458 27476 : if (iDXField < 0 || CSLCount(papszLine) < iDXField + 7)
459 0 : return FALSE;
460 :
461 219808 : for( iField = 0; iField < 7; iField++ )
462 192332 : 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 27476 : if( nMethodCode == 9607 )
470 : {
471 4430 : padfTransform[3] *= -1;
472 4430 : padfTransform[4] *= -1;
473 4430 : padfTransform[5] *= -1;
474 : }
475 :
476 27476 : 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 36208 : 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 36208 : if( nPMCode == 7022 /* PM_Greenwich */ || nPMCode == 8901 )
501 : {
502 35294 : if( pdfOffset != NULL )
503 35294 : *pdfOffset = 0.0;
504 35294 : if( ppszName != NULL )
505 35294 : *ppszName = CPLStrdup( "Greenwich" );
506 35294 : return TRUE;
507 : }
508 :
509 : /* -------------------------------------------------------------------- */
510 : /* Search the database for the corresponding datum code. */
511 : /* -------------------------------------------------------------------- */
512 914 : sprintf( szSearchKey, "%d", nPMCode );
513 :
514 : nUOMAngle =
515 : atoi(CSVGetField( PM_FILENAME,
516 : "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer,
517 914 : "UOM_CODE" ) );
518 914 : if( nUOMAngle < 1 )
519 0 : return FALSE;
520 :
521 : /* -------------------------------------------------------------------- */
522 : /* Get the PM offset. */
523 : /* -------------------------------------------------------------------- */
524 914 : if( pdfOffset != NULL )
525 : {
526 : *pdfOffset =
527 : EPSGAngleStringToDD(
528 : CSVGetField( PM_FILENAME,
529 : "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer,
530 : "GREENWICH_LONGITUDE" ),
531 914 : nUOMAngle );
532 : }
533 :
534 : /* -------------------------------------------------------------------- */
535 : /* Get the name, if requested. */
536 : /* -------------------------------------------------------------------- */
537 914 : if( ppszName != NULL )
538 : *ppszName =
539 : CPLStrdup(
540 : CSVGetField( PM_FILENAME,
541 : "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer,
542 914 : "PRIME_MERIDIAN_NAME" ));
543 :
544 914 : 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 60382 : 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 60382 : pszFilename = CSVFilename("gcs.override.csv");
570 60382 : sprintf( szSearchKey, "%d", nGCSCode );
571 :
572 : nDatum = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
573 : szSearchKey, CC_Integer,
574 60382 : "DATUM_CODE" ) );
575 :
576 60382 : if( nDatum < 1 )
577 : {
578 60382 : pszFilename = CSVFilename("gcs.csv");
579 60382 : sprintf( szSearchKey, "%d", nGCSCode );
580 :
581 : nDatum = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
582 : szSearchKey, CC_Integer,
583 60382 : "DATUM_CODE" ) );
584 : }
585 :
586 60382 : if( nDatum < 1 )
587 24178 : return FALSE;
588 :
589 36204 : if( pnDatum != NULL )
590 36204 : *pnDatum = nDatum;
591 :
592 : /* -------------------------------------------------------------------- */
593 : /* Get the PM. */
594 : /* -------------------------------------------------------------------- */
595 : nPM = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
596 : szSearchKey, CC_Integer,
597 36204 : "PRIME_MERIDIAN_CODE" ) );
598 :
599 36204 : if( nPM < 1 )
600 0 : return FALSE;
601 :
602 36204 : if( pnPM != NULL )
603 36204 : *pnPM = nPM;
604 :
605 : /* -------------------------------------------------------------------- */
606 : /* Get the Ellipsoid. */
607 : /* -------------------------------------------------------------------- */
608 : nEllipsoid = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
609 : szSearchKey, CC_Integer,
610 36204 : "ELLIPSOID_CODE" ) );
611 :
612 36204 : if( nEllipsoid < 1 )
613 0 : return FALSE;
614 :
615 36204 : if( pnEllipsoid != NULL )
616 36204 : *pnEllipsoid = nEllipsoid;
617 :
618 : /* -------------------------------------------------------------------- */
619 : /* Get the angular units. */
620 : /* -------------------------------------------------------------------- */
621 : nUOMAngle = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
622 : szSearchKey, CC_Integer,
623 36204 : "UOM_CODE" ) );
624 :
625 36204 : if( nUOMAngle < 1 )
626 0 : return FALSE;
627 :
628 36204 : if( pnUOMAngle != NULL )
629 36204 : *pnUOMAngle = nUOMAngle;
630 :
631 : /* -------------------------------------------------------------------- */
632 : /* Get the name, if requested. */
633 : /* -------------------------------------------------------------------- */
634 36204 : if( ppszName != NULL )
635 : *ppszName =
636 : CPLStrdup(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
637 : szSearchKey, CC_Integer,
638 36204 : "COORD_REF_SYS_NAME" ));
639 :
640 : /* -------------------------------------------------------------------- */
641 : /* Get the datum name, if requested. */
642 : /* -------------------------------------------------------------------- */
643 36204 : if( ppszDatumName != NULL )
644 : *ppszDatumName =
645 : CPLStrdup(CSVGetField( pszFilename, "COORD_REF_SYS_CODE",
646 : szSearchKey, CC_Integer,
647 36204 : "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 36204 : "COORD_SYS_CODE" ) );
657 :
658 36204 : if( pnCoordSysCode != NULL )
659 36204 : *pnCoordSysCode = nCSC;
660 :
661 36204 : 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 36340 : OSRGetEllipsoidInfo( int nCode, char ** ppszName,
691 : double * pdfSemiMajor, double * pdfInvFlattening )
692 :
693 : {
694 : char szSearchKey[24];
695 36340 : double dfSemiMajor, dfToMeters = 1.0;
696 : int nUOMLength;
697 :
698 : /* -------------------------------------------------------------------- */
699 : /* Get the semi major axis. */
700 : /* -------------------------------------------------------------------- */
701 36340 : snprintf( szSearchKey, sizeof(szSearchKey), "%d", nCode );
702 36340 : szSearchKey[sizeof(szSearchKey) - 1] = '\n';
703 :
704 : dfSemiMajor =
705 : CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ),
706 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
707 36340 : "SEMI_MAJOR_AXIS" ) );
708 36340 : 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 36340 : "UOM_CODE" ));
717 36340 : EPSGGetUOMLengthInfo( nUOMLength, NULL, &dfToMeters );
718 :
719 36340 : dfSemiMajor *= dfToMeters;
720 :
721 36340 : if( pdfSemiMajor != NULL )
722 36340 : *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 36340 : if( pdfInvFlattening != NULL )
729 : {
730 : *pdfInvFlattening =
731 : CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ),
732 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
733 36340 : "INV_FLATTENING" ));
734 :
735 36340 : if( *pdfInvFlattening == 0.0 )
736 : {
737 : double dfSemiMinor;
738 :
739 : dfSemiMinor =
740 : CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ),
741 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
742 3955 : "SEMI_MINOR_AXIS" )) * dfToMeters;
743 :
744 7800 : if( dfSemiMajor != 0.0 && dfSemiMajor != dfSemiMinor )
745 : *pdfInvFlattening =
746 3845 : -1.0 / (dfSemiMinor/dfSemiMajor - 1.0);
747 : else
748 110 : *pdfInvFlattening = 0.0;
749 : }
750 : }
751 :
752 : /* -------------------------------------------------------------------- */
753 : /* Get the name, if requested. */
754 : /* -------------------------------------------------------------------- */
755 36340 : if( ppszName != NULL )
756 : *ppszName =
757 : CPLStrdup(CSVGetField( CSVFilename("ellipsoid.csv" ),
758 : "ELLIPSOID_CODE", szSearchKey, CC_Integer,
759 36244 : "ELLIPSOID_NAME" ));
760 :
761 36340 : 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 23410 : 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 23410 : 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 23410 : osFilename = CSVFilename( "pcs.override.csv" );
816 23410 : sprintf( szTRFCode, "%d", nPCS );
817 : nProjMethod =
818 : atoi( CSVGetField( osFilename,
819 : "COORD_REF_SYS_CODE", szTRFCode, CC_Integer,
820 23410 : "COORD_OP_METHOD_CODE" ) );
821 23410 : if( nProjMethod == 0 )
822 : {
823 23388 : osFilename = CSVFilename( "pcs.csv" );
824 23388 : sprintf( szTRFCode, "%d", nPCS );
825 : nProjMethod =
826 : atoi( CSVGetField( osFilename,
827 : "COORD_REF_SYS_CODE", szTRFCode, CC_Integer,
828 23388 : "COORD_OP_METHOD_CODE" ) );
829 : }
830 :
831 23410 : if( nProjMethod == 0 )
832 0 : return FALSE;
833 :
834 : /* -------------------------------------------------------------------- */
835 : /* Get the parameters for this projection. */
836 : /* -------------------------------------------------------------------- */
837 :
838 187280 : for( i = 0; i < 7; i++ )
839 : {
840 : char szParamUOMID[32], szParamValueID[32], szParamCodeID[32];
841 : char *pszValue;
842 : int nUOM;
843 :
844 163870 : sprintf( szParamCodeID, "PARAMETER_CODE_%d", i+1 );
845 163870 : sprintf( szParamUOMID, "PARAMETER_UOM_%d", i+1 );
846 163870 : sprintf( szParamValueID, "PARAMETER_VALUE_%d", i+1 );
847 :
848 163870 : if( panParmIds != NULL )
849 163870 : panParmIds[i] =
850 : atoi(CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode,
851 163870 : CC_Integer, szParamCodeID ));
852 :
853 : nUOM = atoi(CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode,
854 163870 : CC_Integer, szParamUOMID ));
855 : pszValue = CPLStrdup(
856 : CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode,
857 163870 : 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 454826 : if( (panParmIds[i] == NatOriginScaleFactor
862 145557 : || panParmIds[i] == InitialLineScaleFactor
863 145399 : || panParmIds[i] == PseudoStdParallelScaleFactor)
864 : && nUOM < 9200 )
865 24 : nUOM = 9201;
866 :
867 219596 : if( nUOM >= 9100 && nUOM < 9200 )
868 55726 : adfProjParms[i] = EPSGAngleStringToDD( pszValue, nUOM );
869 154904 : else if( nUOM > 9000 && nUOM < 9100 )
870 : {
871 : double dfInMeters;
872 :
873 46760 : if( !EPSGGetUOMLengthInfo( nUOM, NULL, &dfInMeters ) )
874 0 : dfInMeters = 1.0;
875 46760 : adfProjParms[i] = CPLAtof(pszValue) * dfInMeters;
876 : }
877 61384 : else if( EQUAL(pszValue,"") ) /* null field */
878 : {
879 42581 : adfProjParms[i] = 0.0;
880 : }
881 : else /* really we should consider looking up other scaling factors */
882 : {
883 18803 : if( nUOM != 9201 )
884 : CPLDebug( "OGR",
885 : "Non-unity scale factor units! (UOM=%d, PCS=%d)",
886 180 : nUOM, nPCS );
887 18803 : adfProjParms[i] = CPLAtof(pszValue);
888 : }
889 :
890 163870 : CPLFree( pszValue );
891 : }
892 :
893 : /* -------------------------------------------------------------------- */
894 : /* Transfer requested data into passed variables. */
895 : /* -------------------------------------------------------------------- */
896 23410 : if( pnProjMethod != NULL )
897 23410 : *pnProjMethod = nProjMethod;
898 :
899 23410 : if( padfProjParms != NULL )
900 : {
901 187280 : for( i = 0; i < 7; i++ )
902 163870 : padfProjParms[i] = adfProjParms[i];
903 : }
904 :
905 23410 : return TRUE;
906 : }
907 :
908 :
909 : /************************************************************************/
910 : /* EPSGGetPCSInfo() */
911 : /************************************************************************/
912 :
913 : static int
914 24168 : 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 24168 : pszFilename = CSVFilename( "pcs.csv" );
928 24168 : sprintf( szSearchKey, "%d", nPCSCode );
929 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
930 24168 : szSearchKey, CC_Integer );
931 :
932 24168 : if( papszRecord == NULL )
933 : {
934 740 : pszFilename = CSVFilename( "pcs.override.csv" );
935 740 : sprintf( szSearchKey, "%d", nPCSCode );
936 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
937 740 : szSearchKey, CC_Integer );
938 :
939 : }
940 :
941 24168 : if( papszRecord == NULL )
942 740 : return FALSE;
943 :
944 : /* -------------------------------------------------------------------- */
945 : /* Get the name, if requested. */
946 : /* -------------------------------------------------------------------- */
947 23428 : if( ppszEPSGName != NULL )
948 : {
949 : CPLString osPCSName =
950 : CSLGetField( papszRecord,
951 : CSVGetFileFieldId(pszFilename,
952 23428 : "COORD_REF_SYS_NAME"));
953 :
954 : const char *pszDeprecated =
955 : CSLGetField( papszRecord,
956 : CSVGetFileFieldId(pszFilename,
957 23428 : "DEPRECATED") );
958 :
959 23428 : if( pszDeprecated != NULL && *pszDeprecated == '1' )
960 1786 : osPCSName += " (deprecated)";
961 :
962 23428 : *ppszEPSGName = CPLStrdup(osPCSName);
963 : }
964 :
965 : /* -------------------------------------------------------------------- */
966 : /* Get the UOM Length code, if requested. */
967 : /* -------------------------------------------------------------------- */
968 23428 : if( pnUOMLengthCode != NULL )
969 : {
970 : const char *pszValue;
971 :
972 : pszValue =
973 : CSLGetField( papszRecord,
974 23428 : CSVGetFileFieldId(pszFilename,"UOM_CODE"));
975 23428 : if( atoi(pszValue) > 0 )
976 23428 : *pnUOMLengthCode = atoi(pszValue);
977 : else
978 0 : *pnUOMLengthCode = 0;
979 : }
980 :
981 : /* -------------------------------------------------------------------- */
982 : /* Get the UOM Angle code, if requested. */
983 : /* -------------------------------------------------------------------- */
984 23428 : if( pnUOMAngleCode != NULL )
985 : {
986 : const char *pszValue;
987 :
988 : pszValue =
989 : CSLGetField( papszRecord,
990 23428 : CSVGetFileFieldId(pszFilename,"UOM_ANGLE_CODE") );
991 :
992 23428 : if( atoi(pszValue) > 0 )
993 0 : *pnUOMAngleCode = atoi(pszValue);
994 : else
995 23428 : *pnUOMAngleCode = 0;
996 : }
997 :
998 : /* -------------------------------------------------------------------- */
999 : /* Get the GeogCS (Datum with PM) code, if requested. */
1000 : /* -------------------------------------------------------------------- */
1001 23428 : if( pnGeogCS != NULL )
1002 : {
1003 : const char *pszValue;
1004 :
1005 : pszValue =
1006 : CSLGetField( papszRecord,
1007 23428 : CSVGetFileFieldId(pszFilename,"SOURCE_GEOGCRS_CODE"));
1008 23428 : if( atoi(pszValue) > 0 )
1009 23428 : *pnGeogCS = atoi(pszValue);
1010 : else
1011 0 : *pnGeogCS = 0;
1012 : }
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* Get the GeogCS (Datum with PM) code, if requested. */
1016 : /* -------------------------------------------------------------------- */
1017 23428 : if( pnTRFCode != NULL )
1018 : {
1019 : const char *pszValue;
1020 :
1021 : pszValue =
1022 : CSLGetField( papszRecord,
1023 23428 : CSVGetFileFieldId(pszFilename,"COORD_OP_CODE"));
1024 :
1025 :
1026 23428 : if( atoi(pszValue) > 0 )
1027 23428 : *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 23428 : "COORD_SYS_CODE" ) );
1040 :
1041 23428 : if( pnCoordSysCode != NULL )
1042 23428 : *pnCoordSysCode = nCSC;
1043 :
1044 23428 : return TRUE;
1045 : }
1046 :
1047 : /************************************************************************/
1048 : /* SetEPSGAxisInfo() */
1049 : /************************************************************************/
1050 :
1051 59404 : 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 59404 : if( nCoordSysCode >= 4400 && nCoordSysCode <= 4410 )
1062 : {
1063 : return
1064 : poSRS->SetAxes( pszTargetKey,
1065 : "Easting", OAO_East,
1066 8129 : "Northing", OAO_North );
1067 : }
1068 :
1069 : // Conventional and common Easting/Northing values.
1070 51275 : if( nCoordSysCode >= 6400 && nCoordSysCode <= 6423 )
1071 : {
1072 : return
1073 : poSRS->SetAxes( pszTargetKey,
1074 : "Latitude", OAO_North,
1075 36204 : "Longitude", OAO_East );
1076 : }
1077 :
1078 : /* -------------------------------------------------------------------- */
1079 : /* Get the definition from the coordinate_axis.csv file. */
1080 : /* -------------------------------------------------------------------- */
1081 : char **papszRecord;
1082 15071 : char **papszAxis1=NULL, **papszAxis2=NULL;
1083 : char szSearchKey[24];
1084 : const char *pszFilename;
1085 :
1086 15071 : pszFilename = CSVFilename( "coordinate_axis.csv" );
1087 15071 : sprintf( szSearchKey, "%d", nCoordSysCode );
1088 : papszRecord = CSVScanFileByName( pszFilename, "COORD_SYS_CODE",
1089 15071 : szSearchKey, CC_Integer );
1090 :
1091 15071 : if( papszRecord != NULL )
1092 : {
1093 15071 : papszAxis1 = CSLDuplicate( papszRecord );
1094 15071 : papszRecord = CSVGetNextLine( pszFilename );
1095 30142 : if( CSLCount(papszRecord) > 0
1096 15071 : && EQUAL(papszRecord[0],papszAxis1[0]) )
1097 15071 : papszAxis2 = CSLDuplicate( papszRecord );
1098 : }
1099 :
1100 15071 : 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 15071 : CSVGetFileFieldId( pszFilename, "coord_axis_orientation" );
1118 : iAxisAbbrevField =
1119 15071 : CSVGetFileFieldId( pszFilename, "coord_axis_abbreviation" );
1120 : iAxisOrderField =
1121 15071 : CSVGetFileFieldId( pszFilename, "coord_axis_order" );
1122 : iAxisNameCodeField =
1123 15071 : 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 15071 : 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 15071 : 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 15071 : if( atoi(papszAxis2[iAxisOrderField]) < atoi(papszAxis1[iAxisOrderField]) )
1157 : {
1158 7168 : papszRecord = papszAxis1;
1159 7168 : papszAxis1 = papszAxis2;
1160 7168 : papszAxis2 = papszRecord;
1161 : }
1162 :
1163 : /* -------------------------------------------------------------------- */
1164 : /* Work out axis enumeration values. */
1165 : /* -------------------------------------------------------------------- */
1166 15071 : OGRAxisOrientation eOAxis1 = OAO_Other, eOAxis2 = OAO_Other;
1167 : int iAO;
1168 : static int anCodes[7] = { -1, 9907, 9909, 9906, 9908, -1, -1 };
1169 :
1170 120568 : for( iAO = 0; iAO < 7; iAO++ )
1171 : {
1172 105497 : if( EQUAL(papszAxis1[iAxisOrientationField],
1173 : OSRAxisEnumToName((OGRAxisOrientation) iAO)) )
1174 14813 : eOAxis1 = (OGRAxisOrientation) iAO;
1175 105497 : if( EQUAL(papszAxis2[iAxisOrientationField],
1176 : OSRAxisEnumToName((OGRAxisOrientation) iAO)) )
1177 14813 : eOAxis2 = (OGRAxisOrientation) iAO;
1178 :
1179 168395 : if( eOAxis1 == OAO_Other
1180 62898 : && anCodes[iAO] == atoi(papszAxis1[iAxisNameCodeField]) )
1181 258 : eOAxis1 = (OGRAxisOrientation) iAO;
1182 165455 : if( eOAxis2 == OAO_Other
1183 59958 : && anCodes[iAO] == atoi(papszAxis2[iAxisNameCodeField]) )
1184 258 : 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 15071 : apszAxisName[0] = papszAxis1[iAxisAbbrevField];
1193 15071 : apszAxisName[1] = papszAxis2[iAxisAbbrevField];
1194 :
1195 45213 : for( iAO = 0; iAO < 2; iAO++ )
1196 : {
1197 30142 : if( EQUAL(apszAxisName[iAO],"N") )
1198 1338 : apszAxisName[iAO] = "Northing";
1199 28804 : else if( EQUAL(apszAxisName[iAO],"E") )
1200 1338 : apszAxisName[iAO] = "Easting";
1201 27466 : else if( EQUAL(apszAxisName[iAO],"S") )
1202 0 : apszAxisName[iAO] = "Southing";
1203 27466 : 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 15071 : apszAxisName[1], eOAxis2 );
1214 :
1215 15071 : CSLDestroy( papszAxis1 );
1216 15071 : CSLDestroy( papszAxis2 );
1217 :
1218 15071 : return eResult;
1219 : }
1220 :
1221 : /************************************************************************/
1222 : /* SetEPSGGeogCS() */
1223 : /* */
1224 : /* FLAWS: */
1225 : /* o Units are all hardcoded. */
1226 : /************************************************************************/
1227 :
1228 60382 : static OGRErr SetEPSGGeogCS( OGRSpatialReference * poSRS, int nGeogCS )
1229 :
1230 : {
1231 : int nDatumCode, nPMCode, nUOMAngle, nEllipsoidCode, nCSC;
1232 60382 : char *pszGeogCSName = NULL, *pszDatumName = NULL, *pszEllipsoidName = NULL;
1233 60382 : char *pszPMName = NULL, *pszAngleName = NULL;
1234 : double dfPMOffset, dfSemiMajor, dfInvFlattening, adfBursaTransform[7];
1235 : double dfAngleInDegrees, dfAngleInRadians;
1236 :
1237 60382 : if( !EPSGGetGCSInfo( nGeogCS, &pszGeogCSName,
1238 : &nDatumCode, &pszDatumName,
1239 : &nPMCode, &nEllipsoidCode, &nUOMAngle, &nCSC ) )
1240 24178 : return OGRERR_UNSUPPORTED_SRS;
1241 :
1242 36204 : if( !EPSGGetPMInfo( nPMCode, &pszPMName, &dfPMOffset ) )
1243 : {
1244 0 : CPLFree( pszDatumName );
1245 0 : CPLFree( pszGeogCSName );
1246 0 : return OGRERR_UNSUPPORTED_SRS;
1247 : }
1248 :
1249 36204 : OGREPSGDatumNameMassage( &pszDatumName );
1250 :
1251 36204 : 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 36204 : if( !EPSGGetUOMAngleInfo( nUOMAngle, &pszAngleName, &dfAngleInDegrees ) )
1261 : {
1262 0 : pszAngleName = CPLStrdup("degree");
1263 0 : dfAngleInDegrees = 1.0;
1264 0 : nUOMAngle = -1;
1265 : }
1266 :
1267 36204 : if( dfAngleInDegrees == 1.0 )
1268 35926 : dfAngleInRadians = CPLAtof(SRS_UA_DEGREE_CONV);
1269 : else
1270 278 : dfAngleInRadians = CPLAtof(SRS_UA_DEGREE_CONV) * dfAngleInDegrees;
1271 :
1272 : poSRS->SetGeogCS( pszGeogCSName, pszDatumName,
1273 : pszEllipsoidName, dfSemiMajor, dfInvFlattening,
1274 : pszPMName, dfPMOffset,
1275 36204 : pszAngleName, dfAngleInRadians );
1276 :
1277 36204 : if( EPSGGetWGS84Transform( nGeogCS, adfBursaTransform ) )
1278 : {
1279 : OGR_SRSNode *poWGS84;
1280 : char szValue[100];
1281 :
1282 27476 : poWGS84 = new OGR_SRSNode( "TOWGS84" );
1283 :
1284 439616 : for( int iCoeff = 0; iCoeff < 7; iCoeff++ )
1285 : {
1286 192332 : sprintf( szValue, "%g", adfBursaTransform[iCoeff] );
1287 192332 : poWGS84->AddChild( new OGR_SRSNode( szValue ) );
1288 : }
1289 :
1290 27476 : poSRS->GetAttrNode( "DATUM" )->AddChild( poWGS84 );
1291 : }
1292 :
1293 36204 : poSRS->SetAuthority( "GEOGCS", "EPSG", nGeogCS );
1294 36204 : poSRS->SetAuthority( "DATUM", "EPSG", nDatumCode );
1295 36204 : poSRS->SetAuthority( "SPHEROID", "EPSG", nEllipsoidCode );
1296 36204 : poSRS->SetAuthority( "PRIMEM", "EPSG", nPMCode );
1297 :
1298 36204 : if( nUOMAngle > 0 )
1299 36204 : poSRS->SetAuthority( "GEOGCS|UNIT", "EPSG", nUOMAngle );
1300 :
1301 36204 : CPLFree( pszAngleName );
1302 36204 : CPLFree( pszDatumName );
1303 36204 : CPLFree( pszEllipsoidName );
1304 36204 : CPLFree( pszGeogCSName );
1305 36204 : CPLFree( pszPMName );
1306 :
1307 : /* -------------------------------------------------------------------- */
1308 : /* Set axes */
1309 : /* -------------------------------------------------------------------- */
1310 36204 : if( nCSC > 0 )
1311 : {
1312 36204 : SetEPSGAxisInfo( poSRS, "GEOGCS", nCSC );
1313 36204 : CPLErrorReset();
1314 : }
1315 :
1316 36204 : 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 120372 : 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 120372 : switch( nTargetId )
1337 : {
1338 : case NatOriginScaleFactor:
1339 : case InitialLineScaleFactor:
1340 : case PseudoStdParallelScaleFactor:
1341 18556 : dfResult = 1.0;
1342 18556 : break;
1343 :
1344 : case AngleRectifiedToSkewedGrid:
1345 158 : dfResult = 90.0;
1346 158 : break;
1347 :
1348 : default:
1349 101658 : dfResult = 0.0;
1350 : }
1351 :
1352 : /* -------------------------------------------------------------------- */
1353 : /* Try to find actual value in parameter list. */
1354 : /* -------------------------------------------------------------------- */
1355 377295 : for( i = 0; i < 7; i++ )
1356 : {
1357 376774 : if( panParmIds[i] == nTargetId )
1358 : {
1359 119851 : dfResult = padfProjParms[i];
1360 119851 : 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 120372 : return dfResult;
1392 : }
1393 :
1394 : #define OGR_FP(x) OGR_FetchParm( adfProjParms, anParmIds, (x), \
1395 : dfFromGreenwich )
1396 :
1397 : /************************************************************************/
1398 : /* SetEPSGProjCS() */
1399 : /************************************************************************/
1400 :
1401 24168 : static OGRErr SetEPSGProjCS( OGRSpatialReference * poSRS, int nPCSCode )
1402 :
1403 : {
1404 24168 : int nGCSCode, nUOMAngleCode, nUOMLength, nTRFCode, nProjMethod=0;
1405 24168 : int anParmIds[7], nCSC = 0;
1406 24168 : char *pszPCSName = NULL, *pszUOMLengthName = NULL;
1407 : double adfProjParms[7], dfInMeters, dfFromGreenwich;
1408 : OGRErr nErr;
1409 : OGR_SRSNode *poNode;
1410 :
1411 24168 : if( !EPSGGetPCSInfo( nPCSCode, &pszPCSName, &nUOMLength, &nUOMAngleCode,
1412 : &nGCSCode, &nTRFCode, &nCSC ) )
1413 : {
1414 740 : CPLFree(pszPCSName);
1415 740 : return OGRERR_UNSUPPORTED_SRS;
1416 : }
1417 :
1418 23428 : poSRS->SetNode( "PROJCS", pszPCSName );
1419 :
1420 : /* -------------------------------------------------------------------- */
1421 : /* Set GEOGCS. */
1422 : /* -------------------------------------------------------------------- */
1423 23428 : nErr = SetEPSGGeogCS( poSRS, nGCSCode );
1424 23428 : if( nErr != OGRERR_NONE )
1425 : {
1426 18 : CPLFree(pszPCSName);
1427 18 : return nErr;
1428 : }
1429 :
1430 23410 : dfFromGreenwich = poSRS->GetPrimeMeridian();
1431 :
1432 : /* -------------------------------------------------------------------- */
1433 : /* Set linear units. */
1434 : /* -------------------------------------------------------------------- */
1435 23410 : if( !EPSGGetUOMLengthInfo( nUOMLength, &pszUOMLengthName, &dfInMeters ) )
1436 : {
1437 0 : CPLFree(pszPCSName);
1438 0 : return OGRERR_UNSUPPORTED_SRS;
1439 : }
1440 :
1441 23410 : poSRS->SetLinearUnits( pszUOMLengthName, dfInMeters );
1442 23410 : poSRS->SetAuthority( "PROJCS|UNIT", "EPSG", nUOMLength );
1443 :
1444 23410 : CPLFree( pszUOMLengthName );
1445 23410 : CPLFree( pszPCSName );
1446 :
1447 : /* -------------------------------------------------------------------- */
1448 : /* Set projection and parameters. */
1449 : /* -------------------------------------------------------------------- */
1450 23410 : if( !EPSGGetProjTRFInfo( nPCSCode, &nProjMethod, anParmIds, adfProjParms ))
1451 0 : return OGRERR_UNSUPPORTED_SRS;
1452 :
1453 23410 : 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 468 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1460 468 : 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 4006 : OGR_FP( FalseOriginNorthing ));
1467 4006 : 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 6 : OGR_FP( FalseOriginNorthing ));
1474 6 : 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 175 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1489 :
1490 175 : if( nProjMethod == 1024 || nProjMethod == 9841 ) // override hack for google mercator.
1491 : {
1492 : poSRS->SetExtension( "PROJCS", "PROJ4",
1493 87 : "+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 175 : break;
1496 :
1497 : case 9806:
1498 : poSRS->SetCS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1499 126 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1500 126 : break;
1501 :
1502 : case 9807:
1503 : poSRS->SetTM( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1504 : OGR_FP( NatOriginScaleFactor ),
1505 17265 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1506 17265 : break;
1507 :
1508 : case 9808:
1509 : poSRS->SetTMSO( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1510 : OGR_FP( NatOriginScaleFactor ),
1511 174 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1512 174 : break;
1513 :
1514 : case 9809:
1515 : poSRS->SetOS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1516 : OGR_FP( NatOriginScaleFactor ),
1517 140 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1518 140 : break;
1519 :
1520 : case 9810:
1521 : poSRS->SetPS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1522 : OGR_FP( NatOriginScaleFactor ),
1523 36 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1524 36 : break;
1525 :
1526 : case 9811:
1527 : poSRS->SetNZMG( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1528 12 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1529 12 : 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 78 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1538 :
1539 78 : poNode = poSRS->GetAttrNode( "PROJECTION" )->GetChild( 0 );
1540 78 : if( nProjMethod == 9813 )
1541 6 : poNode->SetValue( SRS_PT_LABORDE_OBLIQUE_MERCATOR );
1542 78 : 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 80 : OGR_FP( ProjCenterNorthing ) );
1559 80 : break;
1560 :
1561 : case 9816:
1562 : poSRS->SetTMG( OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ),
1563 : OGR_FP( FalseOriginEasting ),
1564 6 : OGR_FP( FalseOriginNorthing ) );
1565 6 : break;
1566 :
1567 : case 9818:
1568 : poSRS->SetPolyconic( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1569 30 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1570 30 : break;
1571 :
1572 : case 9819:
1573 : {
1574 140 : double dfCenterLong = OGR_FP( ProjCenterLong );
1575 :
1576 140 : if( dfCenterLong == 0.0 ) // See ticket #2559
1577 140 : dfCenterLong = OGR_FP( PolarLongOrigin );
1578 :
1579 140 : double dfAzimuth = OGR_FP( CoLatConeAxis ); // See ticket #4223
1580 140 : if( dfAzimuth == 0.0 )
1581 0 : dfAzimuth = OGR_FP( Azimuth );
1582 :
1583 : poSRS->SetKrovak( OGR_FP( ProjCenterLat ), dfCenterLong,
1584 : dfAzimuth,
1585 : OGR_FP( PseudoStdParallelLat ),
1586 : OGR_FP( PseudoStdParallelScaleFactor ),
1587 : OGR_FP( ProjCenterEasting ),
1588 140 : OGR_FP( ProjCenterNorthing ) );
1589 : }
1590 140 : break;
1591 :
1592 : case 9820:
1593 : case 1027: /* used by EPSG:2163, 3408, 3409, 3973 and 3974 */
1594 : poSRS->SetLAEA( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ),
1595 78 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1596 78 : break;
1597 :
1598 : case 9821: /* DEPREACTED : this is the spherical form, and really needs different
1599 : equations which give different results but PROJ.4 doesn't
1600 : seem to support the spherical form. */
1601 : poSRS->SetLAEA( OGR_FP( SphericalOriginLat ),
1602 : OGR_FP( SphericalOriginLong ),
1603 0 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1604 0 : break;
1605 :
1606 : case 9822: /* Albers (Conic) Equal Area */
1607 : poSRS->SetACEA( OGR_FP( StdParallel1Lat ),
1608 : OGR_FP( StdParallel2Lat ),
1609 : OGR_FP( FalseOriginLat ),
1610 : OGR_FP( FalseOriginLong ),
1611 : OGR_FP( FalseOriginEasting ),
1612 152 : OGR_FP( FalseOriginNorthing ) );
1613 152 : break;
1614 :
1615 : case 9823: /* Equidistant Cylindrical / Plate Carre / Equirectangular */
1616 : case 1028:
1617 : case 1029:
1618 : poSRS->SetEquirectangular( OGR_FP( NatOriginLat ),
1619 : OGR_FP( NatOriginLong ),
1620 30 : 0.0, 0.0 );
1621 30 : break;
1622 :
1623 : case 9829: /* Polar Stereographic (Variant B) */
1624 : poSRS->SetPS( OGR_FP( PolarLatStdParallel ), OGR_FP(PolarLongOrigin),
1625 : 1.0,
1626 168 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1627 168 : break;
1628 :
1629 : case 9834: /* Lambert Cylindrical Equal Area (Spherical) bug #2659 */
1630 : poSRS->SetCEA( OGR_FP( StdParallel1Lat ), OGR_FP( NatOriginLong ),
1631 18 : OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) );
1632 18 : break;
1633 :
1634 : default:
1635 : CPLDebug( "EPSG", "No WKT support for projection method %d.",
1636 210 : nProjMethod );
1637 210 : return OGRERR_UNSUPPORTED_SRS;
1638 : }
1639 :
1640 : /* -------------------------------------------------------------------- */
1641 : /* Set overall PCS authority code. */
1642 : /* -------------------------------------------------------------------- */
1643 23200 : poSRS->SetAuthority( "PROJCS", "EPSG", nPCSCode );
1644 :
1645 : /* -------------------------------------------------------------------- */
1646 : /* Set axes */
1647 : /* -------------------------------------------------------------------- */
1648 23200 : if( nCSC > 0 )
1649 : {
1650 23200 : SetEPSGAxisInfo( poSRS, "PROJCS", nCSC );
1651 23200 : CPLErrorReset();
1652 : }
1653 :
1654 23200 : return OGRERR_NONE;
1655 : }
1656 :
1657 : /************************************************************************/
1658 : /* SetEPSGVertCS() */
1659 : /************************************************************************/
1660 :
1661 974 : static OGRErr SetEPSGVertCS( OGRSpatialReference * poSRS, int nVertCSCode )
1662 :
1663 : {
1664 : /* -------------------------------------------------------------------- */
1665 : /* Fetch record from the vertcs.csv or override file. */
1666 : /* -------------------------------------------------------------------- */
1667 : char **papszRecord;
1668 : char szSearchKey[24];
1669 : const char *pszFilename;
1670 :
1671 974 : pszFilename = CSVFilename( "vertcs.override.csv" );
1672 974 : sprintf( szSearchKey, "%d", nVertCSCode );
1673 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
1674 974 : szSearchKey, CC_Integer );
1675 :
1676 974 : if( papszRecord == NULL )
1677 : {
1678 970 : pszFilename = CSVFilename( "vertcs.csv" );
1679 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
1680 970 : szSearchKey, CC_Integer );
1681 :
1682 : }
1683 :
1684 974 : if( papszRecord == NULL )
1685 958 : return OGRERR_UNSUPPORTED_SRS;
1686 :
1687 :
1688 : /* -------------------------------------------------------------------- */
1689 : /* Setup the basic VERT_CS. */
1690 : /* -------------------------------------------------------------------- */
1691 : poSRS->SetVertCS(
1692 : CSLGetField( papszRecord,
1693 : CSVGetFileFieldId(pszFilename,
1694 : "COORD_REF_SYS_NAME")),
1695 : CSLGetField( papszRecord,
1696 : CSVGetFileFieldId(pszFilename,
1697 16 : "DATUM_NAME")) );
1698 : /* -------------------------------------------------------------------- */
1699 : /* Setup the VERT_DATUM node. */
1700 : /* -------------------------------------------------------------------- */
1701 : poSRS->SetAuthority( "VERT_CS|VERT_DATUM", "EPSG",
1702 : atoi(CSLGetField( papszRecord,
1703 : CSVGetFileFieldId(pszFilename,
1704 16 : "DATUM_CODE"))) );
1705 :
1706 : /* -------------------------------------------------------------------- */
1707 : /* Should we add a geoidgrids extension node? */
1708 : /* -------------------------------------------------------------------- */
1709 : const char *pszMethod =
1710 : CSLGetField( papszRecord,
1711 16 : CSVGetFileFieldId(pszFilename,"COORD_OP_METHOD_CODE_1"));
1712 16 : if( pszMethod && EQUAL(pszMethod,"9665") )
1713 : {
1714 : const char *pszParm11 =
1715 : CSLGetField( papszRecord,
1716 4 : CSVGetFileFieldId(pszFilename,"PARM_1_1"));
1717 :
1718 4 : poSRS->SetExtension( "VERT_CS|VERT_DATUM", "PROJ4_GRIDS", pszParm11 );
1719 : }
1720 :
1721 : /* -------------------------------------------------------------------- */
1722 : /* Set linear units. */
1723 : /* -------------------------------------------------------------------- */
1724 16 : char *pszUOMLengthName = NULL;
1725 : double dfInMeters;
1726 : int nUOM_CODE = atoi(CSLGetField( papszRecord,
1727 : CSVGetFileFieldId(pszFilename,
1728 16 : "UOM_CODE")));
1729 :
1730 16 : if( !EPSGGetUOMLengthInfo( nUOM_CODE, &pszUOMLengthName, &dfInMeters ) )
1731 : {
1732 : CPLError( CE_Failure, CPLE_AppDefined,
1733 0 : "Failed to lookup UOM CODE %d", nUOM_CODE );
1734 : }
1735 : else
1736 : {
1737 16 : poSRS->SetTargetLinearUnits( "VERT_CS", pszUOMLengthName, dfInMeters );
1738 16 : poSRS->SetAuthority( "VERT_CS|UNIT", "EPSG", nUOM_CODE );
1739 :
1740 16 : CPLFree( pszUOMLengthName );
1741 : }
1742 :
1743 : /* -------------------------------------------------------------------- */
1744 : /* Set overall authority code. */
1745 : /* -------------------------------------------------------------------- */
1746 16 : poSRS->SetAuthority( "VERT_CS", "EPSG", nVertCSCode );
1747 :
1748 16 : return OGRERR_NONE;
1749 : }
1750 :
1751 : /************************************************************************/
1752 : /* SetEPSGCompdCS() */
1753 : /************************************************************************/
1754 :
1755 958 : static OGRErr SetEPSGCompdCS( OGRSpatialReference * poSRS, int nCCSCode )
1756 :
1757 : {
1758 : /* -------------------------------------------------------------------- */
1759 : /* Fetch record from the compdcs.csv or override file. */
1760 : /* -------------------------------------------------------------------- */
1761 958 : char **papszRecord = NULL;
1762 : char szSearchKey[24];
1763 : const char *pszFilename;
1764 :
1765 958 : sprintf( szSearchKey, "%d", nCCSCode );
1766 :
1767 : // So far no override file needed.
1768 : // pszFilename = CSVFilename( "compdcs.override.csv" );
1769 : // papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
1770 : // szSearchKey, CC_Integer );
1771 :
1772 : //if( papszRecord == NULL )
1773 : {
1774 958 : pszFilename = CSVFilename( "compdcs.csv" );
1775 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
1776 958 : szSearchKey, CC_Integer );
1777 :
1778 : }
1779 :
1780 958 : if( papszRecord == NULL )
1781 950 : return OGRERR_UNSUPPORTED_SRS;
1782 :
1783 : /* -------------------------------------------------------------------- */
1784 : /* Fetch subinformation now before anything messes with the */
1785 : /* last loaded record. */
1786 : /* -------------------------------------------------------------------- */
1787 : int nPCSCode = atoi(CSLGetField( papszRecord,
1788 : CSVGetFileFieldId(pszFilename,
1789 8 : "CMPD_HORIZCRS_CODE")));
1790 : int nVertCSCode = atoi(CSLGetField( papszRecord,
1791 : CSVGetFileFieldId(pszFilename,
1792 8 : "CMPD_VERTCRS_CODE")));
1793 :
1794 : /* -------------------------------------------------------------------- */
1795 : /* Set the COMPD_CS node with a name. */
1796 : /* -------------------------------------------------------------------- */
1797 : poSRS->SetNode( "COMPD_CS",
1798 : CSLGetField( papszRecord,
1799 : CSVGetFileFieldId(pszFilename,
1800 8 : "COORD_REF_SYS_NAME")) );
1801 :
1802 : /* -------------------------------------------------------------------- */
1803 : /* Lookup the the projected coordinate system. Can the */
1804 : /* horizontal CRS be a GCS? */
1805 : /* -------------------------------------------------------------------- */
1806 8 : OGRSpatialReference oPCS;
1807 : OGRErr eErr;
1808 :
1809 8 : eErr = SetEPSGProjCS( &oPCS, nPCSCode );
1810 8 : if( eErr != OGRERR_NONE )
1811 : {
1812 : // perhaps it is a GCS?
1813 2 : eErr = SetEPSGGeogCS( &oPCS, nPCSCode );
1814 : }
1815 :
1816 8 : if( eErr != OGRERR_NONE )
1817 : {
1818 0 : return eErr;
1819 : }
1820 :
1821 : poSRS->GetRoot()->AddChild(
1822 8 : oPCS.GetRoot()->Clone() );
1823 :
1824 : /* -------------------------------------------------------------------- */
1825 : /* Lookup the VertCS. */
1826 : /* -------------------------------------------------------------------- */
1827 8 : OGRSpatialReference oVertCS;
1828 8 : eErr = SetEPSGVertCS( &oVertCS, nVertCSCode );
1829 8 : if( eErr != OGRERR_NONE )
1830 0 : return eErr;
1831 :
1832 : poSRS->GetRoot()->AddChild(
1833 8 : oVertCS.GetRoot()->Clone() );
1834 :
1835 : /* -------------------------------------------------------------------- */
1836 : /* Set overall authority code. */
1837 : /* -------------------------------------------------------------------- */
1838 8 : poSRS->SetAuthority( "COMPD_CS", "EPSG", nCCSCode );
1839 :
1840 8 : return OGRERR_NONE;
1841 : }
1842 :
1843 : /************************************************************************/
1844 : /* SetEPSGGeocCS() */
1845 : /************************************************************************/
1846 :
1847 950 : static OGRErr SetEPSGGeocCS( OGRSpatialReference * poSRS, int nGCSCode )
1848 :
1849 : {
1850 : /* -------------------------------------------------------------------- */
1851 : /* Fetch record from the geoccs.csv or override file. */
1852 : /* -------------------------------------------------------------------- */
1853 950 : char **papszRecord = NULL;
1854 : char szSearchKey[24];
1855 : const char *pszFilename;
1856 :
1857 950 : sprintf( szSearchKey, "%d", nGCSCode );
1858 :
1859 : // So far no override file needed.
1860 : // pszFilename = CSVFilename( "compdcs.override.csv" );
1861 : // papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
1862 : // szSearchKey, CC_Integer );
1863 :
1864 : //if( papszRecord == NULL )
1865 : {
1866 950 : pszFilename = CSVFilename( "geoccs.csv" );
1867 : papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
1868 950 : szSearchKey, CC_Integer );
1869 :
1870 : }
1871 :
1872 950 : if( papszRecord == NULL )
1873 946 : return OGRERR_UNSUPPORTED_SRS;
1874 :
1875 : /* -------------------------------------------------------------------- */
1876 : /* Set the GEOCCS node with a name. */
1877 : /* -------------------------------------------------------------------- */
1878 4 : poSRS->Clear();
1879 : poSRS->SetGeocCS( CSLGetField( papszRecord,
1880 : CSVGetFileFieldId(pszFilename,
1881 4 : "COORD_REF_SYS_NAME")) );
1882 :
1883 : /* -------------------------------------------------------------------- */
1884 : /* Get datum related information. */
1885 : /* -------------------------------------------------------------------- */
1886 : int nDatumCode, nEllipsoidCode, nPMCode;
1887 : char *pszDatumName;
1888 :
1889 : nDatumCode = atoi(CSLGetField( papszRecord,
1890 : CSVGetFileFieldId(pszFilename,
1891 4 : "DATUM_CODE")));
1892 :
1893 : pszDatumName =
1894 : CPLStrdup( CSLGetField( papszRecord,
1895 4 : CSVGetFileFieldId(pszFilename,"DATUM_NAME") ) );
1896 4 : OGREPSGDatumNameMassage( &pszDatumName );
1897 :
1898 :
1899 : nEllipsoidCode = atoi(CSLGetField( papszRecord,
1900 : CSVGetFileFieldId(pszFilename,
1901 4 : "ELLIPSOID_CODE")));
1902 :
1903 : nPMCode = atoi(CSLGetField( papszRecord,
1904 : CSVGetFileFieldId(pszFilename,
1905 4 : "PRIME_MERIDIAN_CODE")));
1906 :
1907 : /* -------------------------------------------------------------------- */
1908 : /* Get prime meridian information. */
1909 : /* -------------------------------------------------------------------- */
1910 4 : char *pszPMName = NULL;
1911 4 : double dfPMOffset = 0.0;
1912 :
1913 4 : if( !EPSGGetPMInfo( nPMCode, &pszPMName, &dfPMOffset ) )
1914 : {
1915 0 : CPLFree( pszDatumName );
1916 0 : return OGRERR_UNSUPPORTED_SRS;
1917 : }
1918 :
1919 : /* -------------------------------------------------------------------- */
1920 : /* Get the ellipsoid information. */
1921 : /* -------------------------------------------------------------------- */
1922 4 : char *pszEllipsoidName = NULL;
1923 : double dfSemiMajor, dfInvFlattening;
1924 :
1925 4 : if( OSRGetEllipsoidInfo( nEllipsoidCode, &pszEllipsoidName,
1926 : &dfSemiMajor, &dfInvFlattening ) != OGRERR_NONE )
1927 : {
1928 0 : CPLFree( pszDatumName );
1929 0 : CPLFree( pszPMName );
1930 0 : return OGRERR_UNSUPPORTED_SRS;
1931 : }
1932 :
1933 : /* -------------------------------------------------------------------- */
1934 : /* Setup the spheroid. */
1935 : /* -------------------------------------------------------------------- */
1936 : char szValue[128];
1937 :
1938 4 : OGR_SRSNode *poSpheroid = new OGR_SRSNode( "SPHEROID" );
1939 8 : poSpheroid->AddChild( new OGR_SRSNode( pszEllipsoidName ) );
1940 :
1941 4 : OGRPrintDouble( szValue, dfSemiMajor );
1942 8 : poSpheroid->AddChild( new OGR_SRSNode(szValue) );
1943 :
1944 4 : OGRPrintDouble( szValue, dfInvFlattening );
1945 8 : poSpheroid->AddChild( new OGR_SRSNode(szValue) );
1946 :
1947 4 : CPLFree( pszEllipsoidName );
1948 :
1949 : /* -------------------------------------------------------------------- */
1950 : /* Setup the Datum. */
1951 : /* -------------------------------------------------------------------- */
1952 8 : OGR_SRSNode *poDatum = new OGR_SRSNode( "DATUM" );
1953 8 : poDatum->AddChild( new OGR_SRSNode(pszDatumName) );
1954 4 : poDatum->AddChild( poSpheroid );
1955 :
1956 4 : poSRS->GetRoot()->AddChild( poDatum );
1957 :
1958 4 : CPLFree( pszDatumName );
1959 :
1960 : /* -------------------------------------------------------------------- */
1961 : /* Setup the prime meridian. */
1962 : /* -------------------------------------------------------------------- */
1963 4 : if( dfPMOffset == 0.0 )
1964 4 : strcpy( szValue, "0" );
1965 : else
1966 0 : OGRPrintDouble( szValue, dfPMOffset );
1967 :
1968 4 : OGR_SRSNode *poPM = new OGR_SRSNode( "PRIMEM" );
1969 8 : poPM->AddChild( new OGR_SRSNode( pszPMName ) );
1970 8 : poPM->AddChild( new OGR_SRSNode( szValue ) );
1971 :
1972 4 : poSRS->GetRoot()->AddChild( poPM );
1973 :
1974 4 : CPLFree( pszPMName );
1975 :
1976 : /* -------------------------------------------------------------------- */
1977 : /* Should we try to lookup a datum transform? */
1978 : /* -------------------------------------------------------------------- */
1979 : #ifdef notdef
1980 : if( EPSGGetWGS84Transform( nGeogCS, adfBursaTransform ) )
1981 : {
1982 : OGR_SRSNode *poWGS84;
1983 : char szValue[100];
1984 :
1985 : poWGS84 = new OGR_SRSNode( "TOWGS84" );
1986 :
1987 : for( int iCoeff = 0; iCoeff < 7; iCoeff++ )
1988 : {
1989 : sprintf( szValue, "%g", adfBursaTransform[iCoeff] );
1990 : poWGS84->AddChild( new OGR_SRSNode( szValue ) );
1991 : }
1992 :
1993 : poSRS->GetAttrNode( "DATUM" )->AddChild( poWGS84 );
1994 : }
1995 : #endif
1996 :
1997 : /* -------------------------------------------------------------------- */
1998 : /* Set linear units. */
1999 : /* -------------------------------------------------------------------- */
2000 4 : char *pszUOMLengthName = NULL;
2001 4 : double dfInMeters = 1.0;
2002 : int nUOMLength = atoi(CSLGetField( papszRecord,
2003 : CSVGetFileFieldId(pszFilename,
2004 4 : "UOM_CODE")));
2005 :
2006 4 : if( !EPSGGetUOMLengthInfo( nUOMLength, &pszUOMLengthName, &dfInMeters ) )
2007 : {
2008 0 : return OGRERR_UNSUPPORTED_SRS;
2009 : }
2010 :
2011 4 : poSRS->SetLinearUnits( pszUOMLengthName, dfInMeters );
2012 4 : poSRS->SetAuthority( "GEOCCS|UNIT", "EPSG", nUOMLength );
2013 :
2014 4 : CPLFree( pszUOMLengthName );
2015 :
2016 : /* -------------------------------------------------------------------- */
2017 : /* Set axes */
2018 : /* -------------------------------------------------------------------- */
2019 4 : OGR_SRSNode *poAxis = new OGR_SRSNode( "AXIS" );
2020 :
2021 8 : poAxis->AddChild( new OGR_SRSNode( "Geocentric X" ) );
2022 8 : poAxis->AddChild( new OGR_SRSNode( OSRAxisEnumToName(OAO_Other) ) );
2023 :
2024 4 : poSRS->GetRoot()->AddChild( poAxis );
2025 :
2026 8 : poAxis = new OGR_SRSNode( "AXIS" );
2027 :
2028 8 : poAxis->AddChild( new OGR_SRSNode( "Geocentric Y" ) );
2029 8 : poAxis->AddChild( new OGR_SRSNode( OSRAxisEnumToName(OAO_Other) ) );
2030 :
2031 4 : poSRS->GetRoot()->AddChild( poAxis );
2032 :
2033 8 : poAxis = new OGR_SRSNode( "AXIS" );
2034 :
2035 8 : poAxis->AddChild( new OGR_SRSNode( "Geocentric Z" ) );
2036 8 : poAxis->AddChild( new OGR_SRSNode( OSRAxisEnumToName(OAO_North) ) );
2037 :
2038 4 : poSRS->GetRoot()->AddChild( poAxis );
2039 :
2040 : /* -------------------------------------------------------------------- */
2041 : /* Set the authority codes. */
2042 : /* -------------------------------------------------------------------- */
2043 4 : poSRS->SetAuthority( "DATUM", "EPSG", nDatumCode );
2044 4 : poSRS->SetAuthority( "SPHEROID", "EPSG", nEllipsoidCode );
2045 4 : poSRS->SetAuthority( "PRIMEM", "EPSG", nPMCode );
2046 :
2047 : // if( nUOMAngle > 0 )
2048 : // poSRS->SetAuthority( "GEOGCS|UNIT", "EPSG", nUOMAngle );
2049 :
2050 4 : poSRS->SetAuthority( "GEOCCS", "EPSG", nGCSCode );
2051 :
2052 4 : return OGRERR_NONE;
2053 : }
2054 :
2055 : /************************************************************************/
2056 : /* importFromEPSG() */
2057 : /************************************************************************/
2058 :
2059 : /**
2060 : * \brief Initialize SRS based on EPSG GCS or PCS code.
2061 : *
2062 : * This method will initialize the spatial reference based on the
2063 : * passed in EPSG GCS or PCS code. The coordinate system definitions
2064 : * are normally read from the EPSG derived support files such as
2065 : * pcs.csv, gcs.csv, pcs.override.csv, gcs.override.csv and falling
2066 : * back to search for a PROJ.4 epsg init file or a definition in epsg.wkt.
2067 : *
2068 : * These support files are normally searched for in /usr/local/share/gdal
2069 : * or in the directory identified by the GDAL_DATA configuration option.
2070 : * See CPLFindFile() for details.
2071 : *
2072 : * This method is relatively expensive, and generally involves quite a bit
2073 : * of text file scanning. Reasonable efforts should be made to avoid calling
2074 : * it many times for the same coordinate system.
2075 : *
2076 : * This method is similar to importFromEPSGA() except that EPSG preferred
2077 : * axis ordering will *not* be applied for geographic coordinate systems.
2078 : * EPSG normally defines geographic coordinate systems to use lat/long
2079 : * contrary to typical GIS use).
2080 : *
2081 : * This method is the same as the C function OSRImportFromEPSG().
2082 : *
2083 : * @param nCode a GCS or PCS code from the horizontal coordinate system table.
2084 : *
2085 : * @return OGRERR_NONE on success, or an error code on failure.
2086 : */
2087 :
2088 36194 : OGRErr OGRSpatialReference::importFromEPSG( int nCode )
2089 :
2090 : {
2091 36194 : OGRErr eErr = importFromEPSGA( nCode );
2092 :
2093 : // Strip any GCS axis settings found.
2094 36194 : if( eErr == OGRERR_NONE )
2095 : {
2096 35402 : OGR_SRSNode *poGEOGCS = GetAttrNode( "GEOGCS" );
2097 :
2098 35402 : if( poGEOGCS != NULL )
2099 35390 : poGEOGCS->StripNodes( "AXIS" );
2100 : }
2101 :
2102 36194 : return eErr;
2103 : }
2104 :
2105 : /************************************************************************/
2106 : /* OSRImportFromEPSG() */
2107 : /************************************************************************/
2108 :
2109 : /**
2110 : * \brief Initialize SRS based on EPSG GCS or PCS code.
2111 : *
2112 : * This function is the same as OGRSpatialReference::importFromEPSG().
2113 : */
2114 :
2115 17976 : OGRErr CPL_STDCALL OSRImportFromEPSG( OGRSpatialReferenceH hSRS, int nCode )
2116 :
2117 : {
2118 17976 : VALIDATE_POINTER1( hSRS, "OSRImportFromEPSG", CE_Failure );
2119 :
2120 17976 : return ((OGRSpatialReference *) hSRS)->importFromEPSG( nCode );
2121 : }
2122 :
2123 : /************************************************************************/
2124 : /* importFromEPSGA() */
2125 : /************************************************************************/
2126 :
2127 : /**
2128 : * \brief Initialize SRS based on EPSG GCS or PCS code.
2129 : *
2130 : * This method will initialize the spatial reference based on the
2131 : * passed in EPSG GCS or PCS code.
2132 : *
2133 : * This method is similar to importFromEPSG() except that EPSG preferred
2134 : * axis ordering *will* be applied for geographic coordinate systems.
2135 : * EPSG normally defines geographic coordinate systems to use lat/long
2136 : * contrary to typical GIS use). See OGRSpatialReference::importFromEPSG()
2137 : * for more details on operation of this method.
2138 : *
2139 : * This method is the same as the C function OSRImportFromEPSGA().
2140 : *
2141 : * @param nCode a GCS or PCS code from the horizontal coordinate system table.
2142 : *
2143 : * @return OGRERR_NONE on success, or an error code on failure.
2144 : */
2145 :
2146 36952 : OGRErr OGRSpatialReference::importFromEPSGA( int nCode )
2147 :
2148 : {
2149 : OGRErr eErr;
2150 36952 : CPLLocaleC oLocaleForcer;
2151 :
2152 36952 : bNormInfoSet = FALSE;
2153 :
2154 : /* -------------------------------------------------------------------- */
2155 : /* Clear any existing definition. */
2156 : /* -------------------------------------------------------------------- */
2157 36952 : if( GetRoot() != NULL )
2158 : {
2159 8422 : delete poRoot;
2160 8422 : poRoot = NULL;
2161 : }
2162 :
2163 : /* -------------------------------------------------------------------- */
2164 : /* Verify that we can find the required filename(s). */
2165 : /* -------------------------------------------------------------------- */
2166 36952 : if( CSVScanFileByName( CSVFilename( "gcs.csv" ),
2167 : "COORD_REF_SYS_CODE",
2168 : "4269", CC_Integer ) == NULL )
2169 : {
2170 : CPLError( CE_Failure, CPLE_OpenFailed,
2171 : "Unable to open EPSG support file %s.\n"
2172 : "Try setting the GDAL_DATA environment variable to point to the\n"
2173 : "directory containing EPSG csv files.",
2174 0 : CSVFilename( "gcs.csv" ) );
2175 0 : return OGRERR_FAILURE;
2176 : }
2177 :
2178 : /* -------------------------------------------------------------------- */
2179 : /* Try this as various sorts of objects till one works. */
2180 : /* -------------------------------------------------------------------- */
2181 36952 : eErr = SetEPSGGeogCS( this, nCode );
2182 36952 : if( eErr == OGRERR_UNSUPPORTED_SRS )
2183 24160 : eErr = SetEPSGProjCS( this, nCode );
2184 36952 : if( eErr == OGRERR_UNSUPPORTED_SRS )
2185 966 : eErr = SetEPSGVertCS( this, nCode );
2186 36952 : if( eErr == OGRERR_UNSUPPORTED_SRS )
2187 958 : eErr = SetEPSGCompdCS( this, nCode );
2188 36952 : if( eErr == OGRERR_UNSUPPORTED_SRS )
2189 950 : eErr = SetEPSGGeocCS( this, nCode );
2190 :
2191 : /* -------------------------------------------------------------------- */
2192 : /* If we get it as an unsupported code, try looking it up in */
2193 : /* the epsg.wkt coordinate system dictionary. */
2194 : /* -------------------------------------------------------------------- */
2195 36952 : if( eErr == OGRERR_UNSUPPORTED_SRS )
2196 : {
2197 : char szCode[32];
2198 946 : sprintf( szCode, "%d", nCode );
2199 946 : eErr = importFromDict( "epsg.wkt", szCode );
2200 : }
2201 :
2202 : /* -------------------------------------------------------------------- */
2203 : /* If we get it as an unsupported code, try looking it up in */
2204 : /* the PROJ.4 support file(s). */
2205 : /* -------------------------------------------------------------------- */
2206 36952 : if( eErr == OGRERR_UNSUPPORTED_SRS )
2207 : {
2208 : char szWrkDefn[100];
2209 : char *pszNormalized;
2210 :
2211 792 : sprintf( szWrkDefn, "+init=epsg:%d", nCode );
2212 :
2213 792 : pszNormalized = OCTProj4Normalize( szWrkDefn );
2214 :
2215 792 : if( strstr(pszNormalized,"proj=") != NULL )
2216 0 : eErr = importFromProj4( pszNormalized );
2217 :
2218 792 : CPLFree( pszNormalized );
2219 : }
2220 :
2221 : /* -------------------------------------------------------------------- */
2222 : /* Push in authority information if we were successful, and it */
2223 : /* is not already present. */
2224 : /* -------------------------------------------------------------------- */
2225 : const char *pszAuthName;
2226 :
2227 36952 : if( IsProjected() )
2228 23438 : pszAuthName = GetAuthorityName( "PROJCS" );
2229 : else
2230 13514 : pszAuthName = GetAuthorityName( "GEOGCS" );
2231 :
2232 :
2233 36952 : if( eErr == OGRERR_NONE && pszAuthName == NULL )
2234 : {
2235 14 : if( IsProjected() )
2236 2 : SetAuthority( "PROJCS", "EPSG", nCode );
2237 12 : else if( IsGeographic() )
2238 0 : SetAuthority( "GEOGCS", "EPSG", nCode );
2239 : }
2240 :
2241 : /* -------------------------------------------------------------------- */
2242 : /* Otherwise officially issue an error message. */
2243 : /* -------------------------------------------------------------------- */
2244 36952 : if( eErr == OGRERR_UNSUPPORTED_SRS )
2245 : {
2246 : CPLError( CE_Failure, CPLE_NotSupported,
2247 : "EPSG PCS/GCS code %d not found in EPSG support files. Is this a valid\nEPSG coordinate system?",
2248 792 : nCode );
2249 : }
2250 :
2251 : /* -------------------------------------------------------------------- */
2252 : /* To the extent possible, we want to return the results in as */
2253 : /* close to standard OGC format as possible, so we fixup the */
2254 : /* ordering. */
2255 : /* -------------------------------------------------------------------- */
2256 36952 : if( eErr == OGRERR_NONE )
2257 : {
2258 36160 : eErr = FixupOrdering();
2259 : }
2260 :
2261 36952 : return eErr;
2262 : }
2263 :
2264 : /************************************************************************/
2265 : /* OSRImportFromEPSGA() */
2266 : /************************************************************************/
2267 :
2268 : /**
2269 : * \brief Initialize SRS based on EPSG GCS or PCS code.
2270 : *
2271 : * This function is the same as OGRSpatialReference::importFromEPSGA().
2272 : */
2273 :
2274 6 : OGRErr CPL_STDCALL OSRImportFromEPSGA( OGRSpatialReferenceH hSRS, int nCode )
2275 :
2276 : {
2277 6 : VALIDATE_POINTER1( hSRS, "OSRImportFromEPSGA", CE_Failure );
2278 :
2279 6 : return ((OGRSpatialReference *) hSRS)->importFromEPSGA( nCode );
2280 : }
2281 :
2282 : /************************************************************************/
2283 : /* SetStatePlane() */
2284 : /************************************************************************/
2285 :
2286 : /**
2287 : * \brief Set State Plane projection definition.
2288 : *
2289 : * This will attempt to generate a complete definition of a state plane
2290 : * zone based on generating the entire SRS from the EPSG tables. If the
2291 : * EPSG tables are unavailable, it will produce a stubbed LOCAL_CS definition
2292 : * and return OGRERR_FAILURE.
2293 : *
2294 : * This method is the same as the C function OSRSetStatePlaneWithUnits().
2295 : *
2296 : * @param nZone State plane zone number, in the USGS numbering scheme (as
2297 : * dinstinct from the Arc/Info and Erdas numbering scheme.
2298 : *
2299 : * @param bNAD83 TRUE if the NAD83 zone definition should be used or FALSE
2300 : * if the NAD27 zone definition should be used.
2301 : *
2302 : * @param pszOverrideUnitName Linear unit name to apply overriding the
2303 : * legal definition for this zone.
2304 : *
2305 : * @param dfOverrideUnit Linear unit conversion factor to apply overriding
2306 : * the legal definition for this zone.
2307 : *
2308 : * @return OGRERR_NONE on success, or OGRERR_FAILURE on failure, mostly likely
2309 : * due to the EPSG tables not being accessable.
2310 : */
2311 :
2312 14 : OGRErr OGRSpatialReference::SetStatePlane( int nZone, int bNAD83,
2313 : const char *pszOverrideUnitName,
2314 : double dfOverrideUnit )
2315 :
2316 : {
2317 : int nAdjustedId;
2318 : int nPCSCode;
2319 : char szID[32];
2320 :
2321 : /* -------------------------------------------------------------------- */
2322 : /* Get the index id from stateplane.csv. */
2323 : /* -------------------------------------------------------------------- */
2324 14 : if( bNAD83 )
2325 12 : nAdjustedId = nZone;
2326 : else
2327 2 : nAdjustedId = nZone + 10000;
2328 :
2329 : /* -------------------------------------------------------------------- */
2330 : /* Turn this into a PCS code. We assume there will only be one */
2331 : /* PCS corresponding to each Proj_ code since the proj code */
2332 : /* already effectively indicates NAD27 or NAD83. */
2333 : /* -------------------------------------------------------------------- */
2334 14 : sprintf( szID, "%d", nAdjustedId );
2335 : nPCSCode =
2336 : atoi( CSVGetField( CSVFilename( "stateplane.csv" ),
2337 : "ID", szID, CC_Integer,
2338 14 : "EPSG_PCS_CODE" ) );
2339 14 : if( nPCSCode < 1 )
2340 : {
2341 : char szName[128];
2342 : static int bFailureReported = FALSE;
2343 :
2344 0 : if( !bFailureReported )
2345 : {
2346 0 : bFailureReported = TRUE;
2347 : CPLError( CE_Warning, CPLE_OpenFailed,
2348 : "Unable to find state plane zone in stateplane.csv,\n"
2349 : "likely because the GDAL data files cannot be found. Using\n"
2350 0 : "incomplete definition of state plane zone.\n" );
2351 : }
2352 :
2353 0 : Clear();
2354 0 : if( bNAD83 )
2355 : {
2356 0 : sprintf( szName, "State Plane Zone %d / NAD83", nZone );
2357 0 : SetLocalCS( szName );
2358 0 : SetLinearUnits( SRS_UL_METER, 1.0 );
2359 : }
2360 : else
2361 : {
2362 0 : sprintf( szName, "State Plane Zone %d / NAD27", nZone );
2363 0 : SetLocalCS( szName );
2364 0 : SetLinearUnits( SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV) );
2365 : }
2366 :
2367 0 : return OGRERR_FAILURE;
2368 : }
2369 :
2370 : /* -------------------------------------------------------------------- */
2371 : /* Define based on a full EPSG definition of the zone. */
2372 : /* -------------------------------------------------------------------- */
2373 14 : OGRErr eErr = importFromEPSG( nPCSCode );
2374 :
2375 14 : if( eErr != OGRERR_NONE )
2376 0 : return eErr;
2377 :
2378 : /* -------------------------------------------------------------------- */
2379 : /* Apply units override if required. */
2380 : /* */
2381 : /* We will need to adjust the linear projection parameter to */
2382 : /* match the provided units, and clear the authority code. */
2383 : /* -------------------------------------------------------------------- */
2384 14 : if( dfOverrideUnit != 0.0
2385 : && fabs(dfOverrideUnit - GetLinearUnits()) > 0.0000000001 )
2386 : {
2387 2 : double dfFalseEasting = GetNormProjParm( SRS_PP_FALSE_EASTING );
2388 2 : double dfFalseNorthing= GetNormProjParm( SRS_PP_FALSE_NORTHING);
2389 : OGR_SRSNode *poPROJCS;
2390 :
2391 2 : SetLinearUnits( pszOverrideUnitName, dfOverrideUnit );
2392 :
2393 2 : SetNormProjParm( SRS_PP_FALSE_EASTING, dfFalseEasting );
2394 2 : SetNormProjParm( SRS_PP_FALSE_NORTHING, dfFalseNorthing );
2395 :
2396 2 : poPROJCS = GetAttrNode( "PROJCS" );
2397 2 : if( poPROJCS != NULL && poPROJCS->FindChild( "AUTHORITY" ) != -1 )
2398 : {
2399 2 : poPROJCS->DestroyChild( poPROJCS->FindChild( "AUTHORITY" ) );
2400 : }
2401 : }
2402 :
2403 14 : return OGRERR_NONE;
2404 : }
2405 :
2406 : /************************************************************************/
2407 : /* OSRSetStatePlane() */
2408 : /************************************************************************/
2409 :
2410 : /**
2411 : * \brief Set State Plane projection definition.
2412 : *
2413 : * This function is the same as OGRSpatialReference::SetStatePlane().
2414 : */
2415 :
2416 0 : OGRErr OSRSetStatePlane( OGRSpatialReferenceH hSRS, int nZone, int bNAD83 )
2417 :
2418 : {
2419 0 : VALIDATE_POINTER1( hSRS, "OSRSetStatePlane", CE_Failure );
2420 :
2421 0 : return ((OGRSpatialReference *) hSRS)->SetStatePlane( nZone, bNAD83 );
2422 : }
2423 :
2424 : /************************************************************************/
2425 : /* OSRSetStatePlaneWithUnits() */
2426 : /************************************************************************/
2427 :
2428 : /**
2429 : * \brief Set State Plane projection definition.
2430 : *
2431 : * This function is the same as OGRSpatialReference::SetStatePlane().
2432 : */
2433 :
2434 4 : OGRErr OSRSetStatePlaneWithUnits( OGRSpatialReferenceH hSRS,
2435 : int nZone, int bNAD83,
2436 : const char *pszOverrideUnitName,
2437 : double dfOverrideUnit )
2438 :
2439 : {
2440 4 : VALIDATE_POINTER1( hSRS, "OSRSetStatePlaneWithUnits", CE_Failure );
2441 :
2442 : return ((OGRSpatialReference *) hSRS)->SetStatePlane( nZone, bNAD83,
2443 : pszOverrideUnitName,
2444 4 : dfOverrideUnit );
2445 : }
2446 :
2447 : /************************************************************************/
2448 : /* GetEPSGGeogCS() */
2449 : /* */
2450 : /* Try to establish what the EPSG code for this coordinate */
2451 : /* systems GEOGCS might be. Returns -1 if no reasonable guess */
2452 : /* can be made. */
2453 : /* */
2454 : /* TODO: We really need to do some name lookups. */
2455 : /************************************************************************/
2456 :
2457 290 : int OGRSpatialReference::GetEPSGGeogCS()
2458 :
2459 : {
2460 290 : const char *pszAuthName = GetAuthorityName( "GEOGCS" );
2461 :
2462 : /* -------------------------------------------------------------------- */
2463 : /* Do we already have it? */
2464 : /* -------------------------------------------------------------------- */
2465 290 : if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") )
2466 142 : return atoi(GetAuthorityCode( "GEOGCS" ));
2467 :
2468 : /* -------------------------------------------------------------------- */
2469 : /* Get the datum and geogcs names. */
2470 : /* -------------------------------------------------------------------- */
2471 148 : const char *pszGEOGCS = GetAttrValue( "GEOGCS" );
2472 148 : const char *pszDatum = GetAttrValue( "DATUM" );
2473 :
2474 : // We can only operate on coordinate systems with a geogcs.
2475 148 : if( pszGEOGCS == NULL || pszDatum == NULL )
2476 0 : return -1;
2477 :
2478 : /* -------------------------------------------------------------------- */
2479 : /* Is this a "well known" geographic coordinate system? */
2480 : /* -------------------------------------------------------------------- */
2481 : int bWGS, bNAD;
2482 :
2483 : bWGS = strstr(pszGEOGCS,"WGS") != NULL
2484 : || strstr(pszDatum, "WGS")
2485 : || strstr(pszGEOGCS,"World Geodetic System")
2486 : || strstr(pszGEOGCS,"World_Geodetic_System")
2487 : || strstr(pszDatum, "World Geodetic System")
2488 148 : || strstr(pszDatum, "World_Geodetic_System");
2489 :
2490 : bNAD = strstr(pszGEOGCS,"NAD") != NULL
2491 : || strstr(pszDatum, "NAD")
2492 : || strstr(pszGEOGCS,"North American")
2493 : || strstr(pszGEOGCS,"North_American")
2494 : || strstr(pszDatum, "North American")
2495 148 : || strstr(pszDatum, "North_American");
2496 :
2497 148 : if( bWGS && (strstr(pszGEOGCS,"84") || strstr(pszDatum,"84")) )
2498 96 : return 4326;
2499 :
2500 52 : if( bWGS && (strstr(pszGEOGCS,"72") || strstr(pszDatum,"72")) )
2501 0 : return 4322;
2502 :
2503 52 : if( bNAD && (strstr(pszGEOGCS,"83") || strstr(pszDatum,"83")) )
2504 10 : return 4269;
2505 :
2506 42 : if( bNAD && (strstr(pszGEOGCS,"27") || strstr(pszDatum,"27")) )
2507 2 : return 4267;
2508 :
2509 : /* -------------------------------------------------------------------- */
2510 : /* If we know the datum, associate the most likely GCS with */
2511 : /* it. */
2512 : /* -------------------------------------------------------------------- */
2513 40 : pszAuthName = GetAuthorityName( "GEOGCS|DATUM" );
2514 :
2515 40 : if( pszAuthName != NULL
2516 : && EQUAL(pszAuthName,"epsg")
2517 : && GetPrimeMeridian() == 0.0 )
2518 : {
2519 6 : int nDatum = atoi(GetAuthorityCode("GEOGCS|DATUM"));
2520 :
2521 6 : if( nDatum >= 6000 && nDatum <= 6999 )
2522 6 : return nDatum - 2000;
2523 : }
2524 :
2525 34 : return -1;
2526 : }
2527 :
2528 : /************************************************************************/
2529 : /* AutoIdentifyEPSG() */
2530 : /************************************************************************/
2531 :
2532 : /**
2533 : * \brief Set EPSG authority info if possible.
2534 : *
2535 : * This method inspects a WKT definition, and adds EPSG authority nodes
2536 : * where an aspect of the coordinate system can be easily and safely
2537 : * corresponded with an EPSG identifier. In practice, this method will
2538 : * evolve over time. In theory it can add authority nodes for any object
2539 : * (ie. spheroid, datum, GEOGCS, units, and PROJCS) that could have an
2540 : * authority node. Mostly this is useful to inserting appropriate
2541 : * PROJCS codes for common formulations (like UTM n WGS84).
2542 : *
2543 : * If it success the OGRSpatialReference is updated in place, and the
2544 : * method return OGRERR_NONE. If the method fails to identify the
2545 : * general coordinate system OGRERR_UNSUPPORTED_SRS is returned but no
2546 : * error message is posted via CPLError().
2547 : *
2548 : * This method is the same as the C function OSRAutoIdentifyEPSG().
2549 : *
2550 : * @return OGRERR_NONE or OGRERR_UNSUPPORTED_SRS.
2551 : */
2552 :
2553 104 : OGRErr OGRSpatialReference::AutoIdentifyEPSG()
2554 :
2555 : {
2556 : /* -------------------------------------------------------------------- */
2557 : /* Do we have a GEOGCS node, but no authority? If so, try */
2558 : /* guessing it. */
2559 : /* -------------------------------------------------------------------- */
2560 104 : if( (IsProjected() || IsGeographic())
2561 : && GetAuthorityCode( "GEOGCS" ) == NULL )
2562 : {
2563 54 : int nGCS = GetEPSGGeogCS();
2564 54 : if( nGCS != -1 )
2565 20 : SetAuthority( "GEOGCS", "EPSG", nGCS );
2566 : }
2567 :
2568 : /* -------------------------------------------------------------------- */
2569 : /* Is this a UTM coordinate system with a common GEOGCS? */
2570 : /* -------------------------------------------------------------------- */
2571 : int nZone, bNorth;
2572 104 : if( (nZone = GetUTMZone( &bNorth )) != 0
2573 : && GetAuthorityCode( "PROJCS") == NULL )
2574 : {
2575 : const char *pszAuthName, *pszAuthCode;
2576 :
2577 72 : pszAuthName = GetAuthorityName( "PROJCS|GEOGCS" );
2578 72 : pszAuthCode = GetAuthorityCode( "PROJCS|GEOGCS" );
2579 :
2580 72 : if( pszAuthName == NULL || pszAuthCode == NULL )
2581 : {
2582 : /* don't exactly recognise datum */
2583 : }
2584 50 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4326 )
2585 : { // WGS84
2586 0 : if( bNorth )
2587 0 : SetAuthority( "PROJCS", "EPSG", 32600 + nZone );
2588 : else
2589 0 : SetAuthority( "PROJCS", "EPSG", 32700 + nZone );
2590 : }
2591 100 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4267
2592 : && nZone >= 3 && nZone <= 22 && bNorth )
2593 50 : SetAuthority( "PROJCS", "EPSG", 26700 + nZone ); // NAD27
2594 0 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4269
2595 : && nZone >= 3 && nZone <= 23 && bNorth )
2596 0 : SetAuthority( "PROJCS", "EPSG", 26900 + nZone ); // NAD83
2597 0 : else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4322 )
2598 : { // WGS72
2599 0 : if( bNorth )
2600 0 : SetAuthority( "PROJCS", "EPSG", 32200 + nZone );
2601 : else
2602 0 : SetAuthority( "PROJCS", "EPSG", 32300 + nZone );
2603 : }
2604 : }
2605 :
2606 : /* -------------------------------------------------------------------- */
2607 : /* Return. */
2608 : /* -------------------------------------------------------------------- */
2609 104 : if( IsProjected() && GetAuthorityCode("PROJCS") != NULL )
2610 50 : return OGRERR_NONE;
2611 54 : else if( IsGeographic() && GetAuthorityCode("GEOGCS") != NULL )
2612 4 : return OGRERR_NONE;
2613 : else
2614 50 : return OGRERR_UNSUPPORTED_SRS;
2615 : }
2616 :
2617 : /************************************************************************/
2618 : /* OSRAutoIdentifyEPSG() */
2619 : /************************************************************************/
2620 :
2621 : /**
2622 : * \brief Set EPSG authority info if possible.
2623 : *
2624 : * This function is the same as OGRSpatialReference::AutoIdentifyEPSG().
2625 : */
2626 :
2627 0 : OGRErr OSRAutoIdentifyEPSG( OGRSpatialReferenceH hSRS )
2628 :
2629 : {
2630 0 : VALIDATE_POINTER1( hSRS, "OSRAutoIdentifyEPSG", CE_Failure );
2631 :
2632 0 : return ((OGRSpatialReference *) hSRS)->AutoIdentifyEPSG();
2633 : }
2634 :
2635 : /************************************************************************/
2636 : /* EPSGTreatsAsLatLong() */
2637 : /************************************************************************/
2638 :
2639 : /**
2640 : * \brief This method returns TRUE if EPSG feels this geographic coordinate
2641 : * system should be treated as having lat/long coordinate ordering.
2642 : *
2643 : * Currently this returns TRUE for all geographic coordinate systems
2644 : * with an EPSG code set, and AXIS values set defining it as lat, long.
2645 : * Note that coordinate systems with an EPSG code and no axis settings
2646 : * will be assumed to not be lat/long.
2647 : *
2648 : * FALSE will be returned for all coordinate systems that are not geographic,
2649 : * or that do not have an EPSG code set.
2650 : *
2651 : * This method is the same as the C function OSREPSGTreatsAsLatLong().
2652 : *
2653 : * @return TRUE or FALSE.
2654 : */
2655 :
2656 424 : int OGRSpatialReference::EPSGTreatsAsLatLong()
2657 :
2658 : {
2659 424 : if( !IsGeographic() )
2660 42 : return FALSE;
2661 :
2662 382 : const char *pszAuth = GetAuthorityName( "GEOGCS" );
2663 :
2664 382 : if( pszAuth == NULL || !EQUAL(pszAuth,"EPSG") )
2665 0 : return FALSE;
2666 :
2667 382 : OGR_SRSNode *poFirstAxis = GetAttrNode( "GEOGCS|AXIS" );
2668 :
2669 382 : if( poFirstAxis == NULL )
2670 34 : return FALSE;
2671 :
2672 348 : if( poFirstAxis->GetChildCount() >= 2
2673 : && EQUAL(poFirstAxis->GetChild(1)->GetValue(),"NORTH") )
2674 348 : return TRUE;
2675 :
2676 0 : return FALSE;
2677 : }
2678 :
2679 : /************************************************************************/
2680 : /* OSREPSGTreatsAsLatLong() */
2681 : /************************************************************************/
2682 :
2683 : /**
2684 : * \brief This function returns TRUE if EPSG feels this geographic coordinate
2685 : * system should be treated as having lat/long coordinate ordering.
2686 : *
2687 : * This function is the same as OGRSpatialReference::OSREPSGTreatsAsLatLong().
2688 : */
2689 :
2690 8 : int OSREPSGTreatsAsLatLong( OGRSpatialReferenceH hSRS )
2691 :
2692 : {
2693 8 : VALIDATE_POINTER1( hSRS, "OSREPSGTreatsAsLatLong", CE_Failure );
2694 :
2695 8 : return ((OGRSpatialReference *) hSRS)->EPSGTreatsAsLatLong();
2696 : }
2697 :
|
