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