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