1 /******************************************************************************
2  *
3  * Project:  OpenGIS Simple Features Reference Implementation
4  * Purpose:  OGRSpatialReference translation to/from USGS georeferencing
5  *           information (used in GCTP package).
6  * Author:   Andrey Kiselev, dron@ak4719.spb.edu
7  *
8  ******************************************************************************
9  * Copyright (c) 2004, Andrey Kiselev <dron@ak4719.spb.edu>
10  * Copyright (c) 2008-2009, Even Rouault <even dot rouault at spatialys.com>
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 "cpl_port.h"
32 #include "ogr_srs_api.h"
33 
34 #include <cmath>
35 #include <cstddef>
36 
37 #include "cpl_conv.h"
38 #include "cpl_csv.h"
39 #include "cpl_error.h"
40 #include "cpl_string.h"
41 #include "ogr_core.h"
42 #include "ogr_p.h"
43 #include "ogr_spatialref.h"
44 
45 CPL_CVSID("$Id: ogr_srs_usgs.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
46 
47 /************************************************************************/
48 /*  GCTP projection codes.                                              */
49 /************************************************************************/
50 
51 constexpr long GEO    = 0L;   // Geographic
52 constexpr long UTM    = 1L;   // Universal Transverse Mercator (UTM)
53 constexpr long SPCS   = 2L;   // State Plane Coordinates
54 constexpr long ALBERS = 3L;   // Albers Conical Equal Area
55 constexpr long LAMCC  = 4L;   // Lambert Conformal Conic
56 constexpr long MERCAT = 5L;   // Mercator
57 constexpr long PS     = 6L;   // Polar Stereographic
58 constexpr long POLYC  = 7L;   // Polyconic
59 constexpr long EQUIDC = 8L;   // Equidistant Conic
60 constexpr long TM     = 9L;   // Transverse Mercator
61 constexpr long STEREO = 10L;  // Stereographic
62 constexpr long LAMAZ  = 11L;  // Lambert Azimuthal Equal Area
63 constexpr long AZMEQD = 12L;  // Azimuthal Equidistant
64 constexpr long GNOMON = 13L;  // Gnomonic
65 constexpr long ORTHO  = 14L;  // Orthographic
66 // constexpr long GVNSP  = 15L;  // General Vertical Near-Side Perspective
67 constexpr long SNSOID = 16L;  // Sinusiodal
68 constexpr long EQRECT = 17L;  // Equirectangular
69 constexpr long MILLER = 18L;  // Miller Cylindrical
70 constexpr long VGRINT = 19L;  // Van der Grinten
71 constexpr long HOM    = 20L;  // (Hotine) Oblique Mercator
72 constexpr long ROBIN  = 21L;  // Robinson
73 // constexpr long SOM    = 22L;  // Space Oblique Mercator (SOM)
74 // constexpr long ALASKA = 23L;  // Alaska Conformal
75 // constexpr long GOODE  = 24L;  // Interrupted Goode Homolosine
76 constexpr long MOLL   = 25L;  // Mollweide
77 // constexpr long IMOLL  = 26L;  // Interrupted Mollweide
78 // constexpr long HAMMER = 27L;  // Hammer
79 constexpr long WAGIV  = 28L;  // Wagner IV
80 constexpr long WAGVII = 29L;  // Wagner VII
81 // constexpr long OBEQA  = 30L;  // Oblated Equal Area
82 // constexpr long ISINUS1 = 31L; // Integerized Sinusoidal Grid (the same as 99)
83 // constexpr long CEA    = 97L;  // Cylindrical Equal Area (Grid corners set
84                                  // in meters for EASE grid)
85 // constexpr long BCEA   = 98L;  // Cylindrical Equal Area (Grid corners set
86                                  // in DMS degs for EASE grid)
87 // constexpr long ISINUS = 99L;  // Integerized Sinusoidal Grid
88                                  // (added by Raj Gejjagaraguppe ARC for MODIS)
89 
90 /************************************************************************/
91 /*  GCTP ellipsoid codes.                                               */
92 /************************************************************************/
93 
94 constexpr long CLARKE1866         = 0L;
95 // constexpr long CLARKE1880         = 1L;
96 // constexpr long BESSEL             = 2L;
97 // constexpr long INTERNATIONAL1967  = 3L;
98 // constexpr long INTERNATIONAL1909  = 4L;
99 // constexpr long WGS72              = 5L;
100 // constexpr long EVEREST            = 6L;
101 // constexpr long WGS66              = 7L;
102 constexpr long GRS1980            = 8L;
103 // constexpr long AIRY               = 9L;
104 // constexpr long MODIFIED_EVEREST   = 10L;
105 // constexpr long MODIFIED_AIRY      = 11L;
106 constexpr long WGS84              = 12L;
107 // constexpr long SOUTHEAST_ASIA     = 13L;
108 // constexpr long AUSTRALIAN_NATIONAL= 14L;
109 // constexpr long KRASSOVSKY         = 15L;
110 // constexpr long HOUGH              = 16L;
111 // constexpr long MERCURY1960        = 17L;
112 // constexpr long MODIFIED_MERCURY   = 18L;
113 // constexpr long SPHERE             = 19L;
114 
115 /************************************************************************/
116 /*  Correspondence between GCTP and EPSG ellipsoid codes.               */
117 /************************************************************************/
118 
119 constexpr int aoEllips[] =
120 {
121     7008,   // Clarke, 1866 (NAD1927)
122     7034,   // Clarke, 1880
123     7004,   // Bessel, 1841
124     0,      // FIXME: New International, 1967 --- skipped
125     7022,   // International, 1924 (Hayford, 1909) XXX?
126     7043,   // WGS, 1972
127     7042,   // Everest, 1830
128     7025,   // FIXME: WGS, 1966
129     7019,   // GRS, 1980 (NAD1983)
130     7001,   // Airy, 1830
131     7018,   // Modified Everest
132     7002,   // Modified Airy
133     7030,   // WGS, 1984 (GPS)
134     0,      // FIXME: Southeast Asia --- skipped
135     7003,   // Australian National, 1965
136     7024,   // Krassovsky, 1940
137     7053,   // Hough
138     0,      // FIXME: Mercury, 1960 --- skipped
139     0,      // FIXME: Modified Mercury, 1968 --- skipped
140     7047,   // Sphere, rad 6370997 m (normal sphere)
141     7006,   // Bessel, 1841 (Namibia)
142     7016,   // Everest (Sabah & Sarawak)
143     7044,   // Everest, 1956
144     7056,   // Everest, Malaysia 1969
145     7018,   // Everest, Malay & Singapr 1948
146     0,      // FIXME: Everest, Pakistan --- skipped
147     7022,   // Hayford (International 1924) XXX?
148     7020,   // Helmert 1906
149     7021,   // Indonesian, 1974
150     7036,   // South American, 1969
151     0       // FIXME: WGS 60 --- skipped
152 };
153 
154 #define NUMBER_OF_ELLIPSOIDS    static_cast<int>(CPL_ARRAYSIZE(aoEllips))
155 
156 /************************************************************************/
157 /*                         OSRImportFromUSGS()                          */
158 /************************************************************************/
159 
160 /**
161  * \brief Import coordinate system from USGS projection definition.
162  *
163  * This function is the same as OGRSpatialReference::importFromUSGS().
164  */
OSRImportFromUSGS(OGRSpatialReferenceH hSRS,long iProjsys,long iZone,double * padfPrjParams,long iDatum)165 OGRErr OSRImportFromUSGS( OGRSpatialReferenceH hSRS, long iProjsys,
166                           long iZone, double *padfPrjParams, long iDatum )
167 
168 {
169     VALIDATE_POINTER1( hSRS, "OSRImportFromUSGS", OGRERR_FAILURE );
170 
171     return OGRSpatialReference::FromHandle(hSRS)->importFromUSGS( iProjsys, iZone,
172                                                            padfPrjParams,
173                                                            iDatum );
174 }
175 
OGRSpatialReferenceUSGSUnpackNoOp(double dfVal)176 static double OGRSpatialReferenceUSGSUnpackNoOp( double dfVal )
177 {
178     return dfVal;
179 }
180 
OGRSpatialReferenceUSGSUnpackRadian(double dfVal)181 static double OGRSpatialReferenceUSGSUnpackRadian( double dfVal )
182 {
183     return dfVal * 180.0 / M_PI;
184 }
185 
186 /************************************************************************/
187 /*                          importFromUSGS()                            */
188 /************************************************************************/
189 
190 /**
191  * \brief Import coordinate system from USGS projection definition.
192  *
193  * This method will import projection definition in style, used by USGS GCTP
194  * software. GCTP operates on angles in packed DMS format (see
195  * CPLDecToPackedDMS() function for details), so all angle values (latitudes,
196  * longitudes, azimuths, etc.) specified in the padfPrjParams array should
197  * be in the packed DMS format, unless bAnglesInPackedDMSFormat is set to FALSE.
198  *
199  * This function is the equivalent of the C function OSRImportFromUSGS().
200  * Note that the bAnglesInPackedDMSFormat parameter is only present in the C++
201  * method. The C function assumes bAnglesInPackedFormat = TRUE.
202  *
203  * @param iProjSys Input projection system code, used in GCTP.
204  *
205  * @param iZone Input zone for UTM and State Plane projection systems. For
206  * Southern Hemisphere UTM use a negative zone code. iZone ignored for all
207  * other projections.
208  *
209  * @param padfPrjParams Array of 15 coordinate system parameters. These
210  * parameters differs for different projections.
211  *
212  * <h4>Projection Transformation Package Projection Parameters</h4>
213  * \code{.unparsed}
214  * ----------------------------------------------------------------------------
215  *                         |                    Array Element
216  *  Code & Projection Id   |---------------------------------------------------
217  *                         |   0  |   1  |  2   |  3   |   4   |    5    |6 | 7
218  * ----------------------------------------------------------------------------
219  *  0 Geographic           |      |      |      |      |       |         |  |
220  *  1 U T M                |Lon/Z |Lat/Z |      |      |       |         |  |
221  *  2 State Plane          |      |      |      |      |       |         |  |
222  *  3 Albers Equal Area    |SMajor|SMinor|STDPR1|STDPR2|CentMer|OriginLat|FE|FN
223  *  4 Lambert Conformal C  |SMajor|SMinor|STDPR1|STDPR2|CentMer|OriginLat|FE|FN
224  *  5 Mercator             |SMajor|SMinor|      |      |CentMer|TrueScale|FE|FN
225  *  6 Polar Stereographic  |SMajor|SMinor|      |      |LongPol|TrueScale|FE|FN
226  *  7 Polyconic            |SMajor|SMinor|      |      |CentMer|OriginLat|FE|FN
227  *  8 Equid. Conic A       |SMajor|SMinor|STDPAR|      |CentMer|OriginLat|FE|FN
228  *    Equid. Conic B       |SMajor|SMinor|STDPR1|STDPR2|CentMer|OriginLat|FE|FN
229  *  9 Transverse Mercator  |SMajor|SMinor|Factor|      |CentMer|OriginLat|FE|FN
230  * 10 Stereographic        |Sphere|      |      |      |CentLon|CenterLat|FE|FN
231  * 11 Lambert Azimuthal    |Sphere|      |      |      |CentLon|CenterLat|FE|FN
232  * 12 Azimuthal            |Sphere|      |      |      |CentLon|CenterLat|FE|FN
233  * 13 Gnomonic             |Sphere|      |      |      |CentLon|CenterLat|FE|FN
234  * 14 Orthographic         |Sphere|      |      |      |CentLon|CenterLat|FE|FN
235  * 15 Gen. Vert. Near Per  |Sphere|      |Height|      |CentLon|CenterLat|FE|FN
236  * 16 Sinusoidal           |Sphere|      |      |      |CentMer|         |FE|FN
237  * 17 Equirectangular      |Sphere|      |      |      |CentMer|TrueScale|FE|FN
238  * 18 Miller Cylindrical   |Sphere|      |      |      |CentMer|         |FE|FN
239  * 19 Van der Grinten      |Sphere|      |      |      |CentMer|OriginLat|FE|FN
240  * 20 Hotin Oblique Merc A |SMajor|SMinor|Factor|      |       |OriginLat|FE|FN
241  *    Hotin Oblique Merc B |SMajor|SMinor|Factor|AziAng|AzmthPt|OriginLat|FE|FN
242  * 21 Robinson             |Sphere|      |      |      |CentMer|         |FE|FN
243  * 22 Space Oblique Merc A |SMajor|SMinor|      |IncAng|AscLong|         |FE|FN
244  *    Space Oblique Merc B |SMajor|SMinor|Satnum|Path  |       |         |FE|FN
245  * 23 Alaska Conformal     |SMajor|SMinor|      |      |       |         |FE|FN
246  * 24 Interrupted Goode    |Sphere|      |      |      |       |         |  |
247  * 25 Mollweide            |Sphere|      |      |      |CentMer|         |FE|FN
248  * 26 Interrupt Mollweide  |Sphere|      |      |      |       |         |  |
249  * 27 Hammer               |Sphere|      |      |      |CentMer|         |FE|FN
250  * 28 Wagner IV            |Sphere|      |      |      |CentMer|         |FE|FN
251  * 29 Wagner VII           |Sphere|      |      |      |CentMer|         |FE|FN
252  * 30 Oblated Equal Area   |Sphere|      |Shapem|Shapen|CentLon|CenterLat|FE|FN
253  * ----------------------------------------------------------------------------
254  *
255  *       ----------------------------------------------------
256  *                               |      Array Element       |
257  *         Code & Projection Id  |---------------------------
258  *                               |  8  |  9 |  10 | 11 | 12 |
259  *       ----------------------------------------------------
260  *        0 Geographic           |     |    |     |    |    |
261  *        1 U T M                |     |    |     |    |    |
262  *        2 State Plane          |     |    |     |    |    |
263  *        3 Albers Equal Area    |     |    |     |    |    |
264  *        4 Lambert Conformal C  |     |    |     |    |    |
265  *        5 Mercator             |     |    |     |    |    |
266  *        6 Polar Stereographic  |     |    |     |    |    |
267  *        7 Polyconic            |     |    |     |    |    |
268  *        8 Equid. Conic A       |zero |    |     |    |    |
269  *          Equid. Conic B       |one  |    |     |    |    |
270  *        9 Transverse Mercator  |     |    |     |    |    |
271  *       10 Stereographic        |     |    |     |    |    |
272  *       11 Lambert Azimuthal    |     |    |     |    |    |
273  *       12 Azimuthal            |     |    |     |    |    |
274  *       13 Gnomonic             |     |    |     |    |    |
275  *       14 Orthographic         |     |    |     |    |    |
276  *       15 Gen. Vert. Near Per  |     |    |     |    |    |
277  *       16 Sinusoidal           |     |    |     |    |    |
278  *       17 Equirectangular      |     |    |     |    |    |
279  *       18 Miller Cylindrical   |     |    |     |    |    |
280  *       19 Van der Grinten      |     |    |     |    |    |
281  *       20 Hotin Oblique Merc A |Long1|Lat1|Long2|Lat2|zero|
282  *          Hotin Oblique Merc B |     |    |     |    |one |
283  *       21 Robinson             |     |    |     |    |    |
284  *       22 Space Oblique Merc A |PSRev|LRat|PFlag|    |zero|
285  *          Space Oblique Merc B |     |    |     |    |one |
286  *       23 Alaska Conformal     |     |    |     |    |    |
287  *       24 Interrupted Goode    |     |    |     |    |    |
288  *       25 Mollweide            |     |    |     |    |    |
289  *       26 Interrupt Mollweide  |     |    |     |    |    |
290  *       27 Hammer               |     |    |     |    |    |
291  *       28 Wagner IV            |     |    |     |    |    |
292  *       29 Wagner VII           |     |    |     |    |    |
293  *       30 Oblated Equal Area   |Angle|    |     |    |    |
294  *       ----------------------------------------------------
295  *
296  *   where
297  *
298  *    Lon/Z     Longitude of any point in the UTM zone or zero.  If zero,
299  *              a zone code must be specified.
300  *    Lat/Z     Latitude of any point in the UTM zone or zero.  If zero, a
301  *              zone code must be specified.
302  *    SMajor    Semi-major axis of ellipsoid.  If zero, Clarke 1866 in meters
303  *              is assumed.
304  *    SMinor    Eccentricity squared of the ellipsoid if less than zero,
305  *              if zero, a spherical form is assumed, or if greater than
306  *              zero, the semi-minor axis of ellipsoid.
307  *    Sphere    Radius of reference sphere.  If zero, 6370997 meters is used.
308  *    STDPAR    Latitude of the standard parallel
309  *    STDPR1    Latitude of the first standard parallel
310  *    STDPR2    Latitude of the second standard parallel
311  *    CentMer   Longitude of the central meridian
312  *    OriginLat Latitude of the projection origin
313  *    FE        False easting in the same units as the semi-major axis
314  *    FN        False northing in the same units as the semi-major axis
315  *    TrueScale Latitude of true scale
316  *    LongPol   Longitude down below pole of map
317  *    Factor    Scale factor at central meridian (Transverse Mercator) or
318  *              center of projection (Hotine Oblique Mercator)
319  *    CentLon   Longitude of center of projection
320  *    CenterLat Latitude of center of projection
321  *    Height    Height of perspective point
322  *    Long1     Longitude of first point on center line (Hotine Oblique
323  *              Mercator, format A)
324  *    Long2     Longitude of second point on center line (Hotine Oblique
325  *              Mercator, format A)
326  *    Lat1      Latitude of first point on center line (Hotine Oblique
327  *              Mercator, format A)
328  *    Lat2      Latitude of second point on center line (Hotine Oblique
329  *              Mercator, format A)
330  *    AziAng    Azimuth angle east of north of center line (Hotine Oblique
331  *              Mercator, format B)
332  *    AzmthPt   Longitude of point on central meridian where azimuth occurs
333  *              (Hotine Oblique Mercator, format B)
334  *    IncAng    Inclination of orbit at ascending node, counter-clockwise
335  *              from equator (SOM, format A)
336  *    AscLong   Longitude of ascending orbit at equator (SOM, format A)
337  *    PSRev     Period of satellite revolution in minutes (SOM, format A)
338  *    LRat      Landsat ratio to compensate for confusion at northern end
339  *              of orbit (SOM, format A -- use 0.5201613)
340  *    PFlag     End of path flag for Landsat:  0 = start of path,
341  *              1 = end of path (SOM, format A)
342  *    Satnum    Landsat Satellite Number (SOM, format B)
343  *    Path      Landsat Path Number (Use WRS-1 for Landsat 1, 2 and 3 and
344  *              WRS-2 for Landsat 4, 5 and 6.)  (SOM, format B)
345  *    Shapem    Oblated Equal Area oval shape parameter m
346  *    Shapen    Oblated Equal Area oval shape parameter n
347  *    Angle     Oblated Equal Area oval rotation angle
348  *
349  * Array elements 13 and 14 are set to zero. All array elements with blank
350  * fields are set to zero too.
351  * \endcode
352  *
353  * @param iDatum Input spheroid.<p>
354  *
355  * If the datum code is negative, the first two values in the parameter array
356  * (param) are used to define the values as follows:
357  *
358  * <ul>
359  *
360  * <li> If padfPrjParams[0] is a non-zero value and padfPrjParams[1] is
361  * greater than one, the semimajor axis is set to padfPrjParams[0] and
362  * the semiminor axis is set to padfPrjParams[1].
363  *
364  * <li> If padfPrjParams[0] is nonzero and padfPrjParams[1] is greater than
365  * zero but less than or equal to one, the semimajor axis is set to
366  * padfPrjParams[0] and the semiminor axis is computed from the eccentricity
367  * squared value padfPrjParams[1]:<p>
368  *
369  * semiminor = sqrt(1.0 - ES) * semimajor<p>
370  *
371  * where<p>
372  *
373  * ES = eccentricity squared
374  *
375  * <li> If padfPrjParams[0] is nonzero and padfPrjParams[1] is equal to zero,
376  * the semimajor axis and semiminor axis are set to padfPrjParams[0].
377  *
378  * <li> If padfPrjParams[0] equals zero and padfPrjParams[1] is greater than
379  * zero, the default Clarke 1866 is used to assign values to the semimajor
380  * axis and semiminor axis.
381  *
382  * <li> If padfPrjParams[0] and padfPrjParams[1] equals zero, the semimajor
383  * axis is set to 6370997.0 and the semiminor axis is set to zero.
384  *
385  * </ul>
386  *
387  * If a datum code is zero or greater, the semimajor and semiminor axis are
388  * defined by the datum code as found in the following table:
389  *
390  * <h4>Supported Datums</h4>
391  * \code{.unparsed}
392  *       0: Clarke 1866 (default)
393  *       1: Clarke 1880
394  *       2: Bessel
395  *       3: International 1967
396  *       4: International 1909
397  *       5: WGS 72
398  *       6: Everest
399  *       7: WGS 66
400  *       8: GRS 1980/WGS 84
401  *       9: Airy
402  *      10: Modified Everest
403  *      11: Modified Airy
404  *      12: WGS 84
405  *      13: Southeast Asia
406  *      14: Australian National
407  *      15: Krassovsky
408  *      16: Hough
409  *      17: Mercury 1960
410  *      18: Modified Mercury 1968
411  *      19: Sphere of Radius 6370997 meters
412  * \endcode
413  *
414  * @param nUSGSAngleFormat one of USGS_ANGLE_DECIMALDEGREES,
415  *    USGS_ANGLE_PACKEDDMS, or USGS_ANGLE_RADIANS (default is
416  *    USGS_ANGLE_PACKEDDMS).
417  *
418  * @return OGRERR_NONE on success or an error code in case of failure.
419  */
420 
importFromUSGS(long iProjSys,long iZone,double * padfPrjParams,long iDatum,int nUSGSAngleFormat)421 OGRErr OGRSpatialReference::importFromUSGS( long iProjSys, long iZone,
422                                             double *padfPrjParams,
423                                             long iDatum,
424                                             int nUSGSAngleFormat )
425 
426 {
427     if( !padfPrjParams )
428         return OGRERR_CORRUPT_DATA;
429 
430     double (*pfnUnpackAnglesFn)(double) = nullptr;
431 
432     if( nUSGSAngleFormat == USGS_ANGLE_DECIMALDEGREES )
433         pfnUnpackAnglesFn = OGRSpatialReferenceUSGSUnpackNoOp;
434     else if( nUSGSAngleFormat == USGS_ANGLE_RADIANS )
435         pfnUnpackAnglesFn = OGRSpatialReferenceUSGSUnpackRadian;
436     else
437         pfnUnpackAnglesFn = CPLPackedDMSToDec;
438 
439 /* -------------------------------------------------------------------- */
440 /*      Operate on the basis of the projection code.                    */
441 /* -------------------------------------------------------------------- */
442     switch( iProjSys )
443     {
444         case GEO:
445             break;
446 
447         case UTM:
448             {
449                 int bNorth = TRUE;
450 
451                 if( !iZone )
452                 {
453                     if( padfPrjParams[2] != 0.0 )
454                     {
455                         iZone = static_cast<long>(padfPrjParams[2]);
456                     }
457                     else if( padfPrjParams[0] != 0.0 &&
458                              padfPrjParams[1] != 0.0 )
459                     {
460                         const double dfUnpackedAngle =
461                             pfnUnpackAnglesFn(padfPrjParams[0]);
462                         iZone = static_cast<long>(
463                             ((dfUnpackedAngle + 180.0) / 6.0) + 1.0);
464                         if( dfUnpackedAngle < 0 )
465                             bNorth = FALSE;
466                     }
467                 }
468 
469                 if( iZone < -60 || iZone > 60 )
470                     return OGRERR_CORRUPT_DATA;
471 
472                 if( iZone < 0 )
473                 {
474                     iZone = -iZone;
475                     bNorth = FALSE;
476                 }
477                 SetUTM( static_cast<int>(iZone), bNorth );
478             }
479             break;
480 
481         case SPCS:
482             {
483                 int bNAD83 = TRUE;
484 
485                 if( iDatum == 0 )
486                     bNAD83 = FALSE;
487                 else if( iDatum != 8 )
488                     CPLError( CE_Warning, CPLE_AppDefined,
489                               "Wrong datum for State Plane projection %d. "
490                               "Should be 0 or 8.", static_cast<int>(iDatum) );
491 
492                 SetStatePlane( static_cast<int>(iZone), bNAD83 );
493             }
494             break;
495 
496         case ALBERS:
497             SetACEA( pfnUnpackAnglesFn(padfPrjParams[2]),
498                      pfnUnpackAnglesFn(padfPrjParams[3]),
499                      pfnUnpackAnglesFn(padfPrjParams[5]),
500                      pfnUnpackAnglesFn(padfPrjParams[4]),
501                      padfPrjParams[6], padfPrjParams[7] );
502             break;
503 
504         case LAMCC:
505             SetLCC( pfnUnpackAnglesFn(padfPrjParams[2]),
506                     pfnUnpackAnglesFn(padfPrjParams[3]),
507                     pfnUnpackAnglesFn(padfPrjParams[5]),
508                     pfnUnpackAnglesFn(padfPrjParams[4]),
509                     padfPrjParams[6], padfPrjParams[7] );
510             break;
511 
512         case MERCAT:
513             SetMercator( pfnUnpackAnglesFn(padfPrjParams[5]),
514                          pfnUnpackAnglesFn(padfPrjParams[4]),
515                          1.0,
516                          padfPrjParams[6], padfPrjParams[7] );
517             break;
518 
519         case PS:
520             SetPS( pfnUnpackAnglesFn(padfPrjParams[5]),
521                    pfnUnpackAnglesFn(padfPrjParams[4]),
522                    1.0,
523                    padfPrjParams[6], padfPrjParams[7] );
524 
525             break;
526 
527         case POLYC:
528             SetPolyconic( pfnUnpackAnglesFn(padfPrjParams[5]),
529                           pfnUnpackAnglesFn(padfPrjParams[4]),
530                           padfPrjParams[6], padfPrjParams[7] );
531             break;
532 
533         case EQUIDC:
534             if( padfPrjParams[8] != 0.0 )
535             {
536                 SetEC( pfnUnpackAnglesFn(padfPrjParams[2]),
537                        pfnUnpackAnglesFn(padfPrjParams[3]),
538                        pfnUnpackAnglesFn(padfPrjParams[5]),
539                        pfnUnpackAnglesFn(padfPrjParams[4]),
540                        padfPrjParams[6], padfPrjParams[7] );
541             }
542             else
543             {
544                 SetEC( pfnUnpackAnglesFn(padfPrjParams[2]),
545                        pfnUnpackAnglesFn(padfPrjParams[2]),
546                        pfnUnpackAnglesFn(padfPrjParams[5]),
547                        pfnUnpackAnglesFn(padfPrjParams[4]),
548                        padfPrjParams[6], padfPrjParams[7] );
549             }
550             break;
551 
552         case TM:
553             SetTM( pfnUnpackAnglesFn(padfPrjParams[5]),
554                    pfnUnpackAnglesFn(padfPrjParams[4]),
555                    padfPrjParams[2],
556                    padfPrjParams[6], padfPrjParams[7] );
557             break;
558 
559         case STEREO:
560             SetStereographic( pfnUnpackAnglesFn(padfPrjParams[5]),
561                               pfnUnpackAnglesFn(padfPrjParams[4]),
562                               1.0,
563                               padfPrjParams[6], padfPrjParams[7] );
564             break;
565 
566         case LAMAZ:
567             SetLAEA( pfnUnpackAnglesFn(padfPrjParams[5]),
568                      pfnUnpackAnglesFn(padfPrjParams[4]),
569                      padfPrjParams[6], padfPrjParams[7] );
570             break;
571 
572         case AZMEQD:
573             SetAE( pfnUnpackAnglesFn(padfPrjParams[5]),
574                    pfnUnpackAnglesFn(padfPrjParams[4]),
575                    padfPrjParams[6], padfPrjParams[7] );
576             break;
577 
578         case GNOMON:
579             SetGnomonic( pfnUnpackAnglesFn(padfPrjParams[5]),
580                          pfnUnpackAnglesFn(padfPrjParams[4]),
581                          padfPrjParams[6], padfPrjParams[7] );
582             break;
583 
584         case ORTHO:
585             SetOrthographic( pfnUnpackAnglesFn(padfPrjParams[5]),
586                              pfnUnpackAnglesFn(padfPrjParams[4]),
587                              padfPrjParams[6], padfPrjParams[7] );
588             break;
589 
590         // FIXME: GVNSP --- General Vertical Near-Side Perspective skipped.
591 
592         case SNSOID:
593             SetSinusoidal( pfnUnpackAnglesFn(padfPrjParams[4]),
594                            padfPrjParams[6], padfPrjParams[7] );
595             break;
596 
597         case EQRECT:
598             SetEquirectangular2( 0.0,
599                                  pfnUnpackAnglesFn(padfPrjParams[4]),
600                                  pfnUnpackAnglesFn(padfPrjParams[5]),
601                                  padfPrjParams[6], padfPrjParams[7] );
602             break;
603 
604         case MILLER:
605             SetMC( pfnUnpackAnglesFn(padfPrjParams[5]),
606                    pfnUnpackAnglesFn(padfPrjParams[4]),
607                    padfPrjParams[6], padfPrjParams[7] );
608             break;
609 
610         case VGRINT:
611             SetVDG( pfnUnpackAnglesFn(padfPrjParams[4]),
612                     padfPrjParams[6], padfPrjParams[7] );
613             break;
614 
615         case HOM:
616             if( padfPrjParams[12] != 0.0 )
617             {
618                 SetHOM( pfnUnpackAnglesFn(padfPrjParams[5]),
619                         pfnUnpackAnglesFn(padfPrjParams[4]),
620                         pfnUnpackAnglesFn(padfPrjParams[3]),
621                         0.0, padfPrjParams[2],
622                         padfPrjParams[6], padfPrjParams[7] );
623             }
624             else
625             {
626                 SetHOM2PNO( pfnUnpackAnglesFn(padfPrjParams[5]),
627                             pfnUnpackAnglesFn(padfPrjParams[9]),
628                             pfnUnpackAnglesFn(padfPrjParams[8]),
629                             pfnUnpackAnglesFn(padfPrjParams[11]),
630                             pfnUnpackAnglesFn(padfPrjParams[10]),
631                             padfPrjParams[2],
632                             padfPrjParams[6], padfPrjParams[7] );
633             }
634             break;
635 
636         case ROBIN:
637             SetRobinson( pfnUnpackAnglesFn(padfPrjParams[4]),
638                          padfPrjParams[6], padfPrjParams[7] );
639             break;
640 
641         // FIXME: SOM --- Space Oblique Mercator skipped.
642         // FIXME: ALASKA --- Alaska Conformal skipped.
643         // FIXME: GOODE --- Interrupted Goode skipped.
644 
645         case MOLL:
646             SetMollweide( pfnUnpackAnglesFn(padfPrjParams[4]),
647                           padfPrjParams[6], padfPrjParams[7] );
648             break;
649 
650         // FIXME: IMOLL --- Interrupted Mollweide skipped.
651         // FIXME: HAMMER --- Hammer skipped.
652 
653         case WAGIV:
654             SetWagner( 4, 0.0, padfPrjParams[6], padfPrjParams[7] );
655             break;
656 
657         case WAGVII:
658             SetWagner( 7, 0.0, padfPrjParams[6], padfPrjParams[7] );
659             break;
660 
661         // FIXME: OBEQA --- Oblated Equal Area skipped.
662         // FIXME: ISINUS1 --- Integerized Sinusoidal Grid (the same as 99).
663         // FIXME: CEA --- Cylindrical Equal Area skipped (Grid corners set in
664         // meters for EASE grid).
665         // FIXME: BCEA --- Cylindrical Equal Area skipped (Grid corners set in
666         // DMS degs for EASE grid).
667         // FIXME: ISINUS --- Integrized Sinusoidal skipped.
668 
669         default:
670             CPLDebug( "OSR_USGS", "Unsupported projection: %ld", iProjSys );
671             SetLocalCS( CPLString().Printf("GCTP projection number %ld",
672                                            iProjSys) );
673             break;
674     }
675 
676 /* -------------------------------------------------------------------- */
677 /*      Try to translate the datum/spheroid.                            */
678 /* -------------------------------------------------------------------- */
679 
680     if( !IsLocal() )
681     {
682         char *pszName = nullptr;
683         double dfSemiMajor = 0.0;
684         double dfInvFlattening = 0.0;
685 
686         if( iDatum < 0 ) // Use specified ellipsoid parameters.
687         {
688             if( padfPrjParams[0] > 0.0 )
689             {
690                 if( padfPrjParams[1] > 1.0 )
691                 {
692                     dfInvFlattening = OSRCalcInvFlattening(padfPrjParams[0],
693                                                            padfPrjParams[1]);
694                 }
695                 else if( padfPrjParams[1] > 0.0 )
696                 {
697                     dfInvFlattening =
698                         1.0 / ( 1.0 - sqrt(1.0 - padfPrjParams[1]) );
699                 }
700                 else
701                 {
702                     dfInvFlattening = 0.0;
703                 }
704 
705                 SetGeogCS( "Unknown datum based upon the custom spheroid",
706                            "Not specified (based on custom spheroid)",
707                            "Custom spheroid", padfPrjParams[0], dfInvFlattening,
708                            nullptr, 0, nullptr, 0 );
709             }
710             else if( padfPrjParams[1] > 0.0 )  // Clarke 1866.
711             {
712                 if( OSRGetEllipsoidInfo( 7008, &pszName, &dfSemiMajor,
713                                           &dfInvFlattening ) == OGRERR_NONE )
714                 {
715                     SetGeogCS( CPLString().Printf(
716                                     "Unknown datum based upon the %s ellipsoid",
717                                     pszName ),
718                                CPLString().Printf(
719                                     "Not specified (based on %s spheroid)",
720                                     pszName ),
721                                pszName, dfSemiMajor, dfInvFlattening,
722                                nullptr, 0.0, nullptr, 0.0 );
723                     SetAuthority( "SPHEROID", "EPSG", 7008 );
724                 }
725             }
726             else                              // Sphere, rad 6370997 m
727             {
728                 if( OSRGetEllipsoidInfo( 7047, &pszName, &dfSemiMajor,
729                                          &dfInvFlattening ) == OGRERR_NONE )
730                 {
731                     SetGeogCS( CPLString().Printf(
732                                     "Unknown datum based upon the %s ellipsoid",
733                                     pszName ),
734                                CPLString().Printf(
735                                     "Not specified (based on %s spheroid)",
736                                     pszName ),
737                                pszName, dfSemiMajor, dfInvFlattening,
738                                nullptr, 0.0, nullptr, 0.0 );
739                     SetAuthority( "SPHEROID", "EPSG", 7047 );
740                 }
741             }
742         }
743         else if( iDatum < NUMBER_OF_ELLIPSOIDS && aoEllips[iDatum] )
744         {
745             if( OSRGetEllipsoidInfo( aoEllips[iDatum], &pszName,
746                                      &dfSemiMajor,
747                                      &dfInvFlattening ) == OGRERR_NONE )
748             {
749                 SetGeogCS( CPLString().Printf(
750                                "Unknown datum based upon the %s ellipsoid",
751                                pszName),
752                            CPLString().Printf(
753                                "Not specified (based on %s spheroid)",
754                                pszName),
755                            pszName, dfSemiMajor, dfInvFlattening,
756                            nullptr, 0.0, nullptr, 0.0 );
757                 SetAuthority( "SPHEROID", "EPSG", aoEllips[iDatum] );
758             }
759             else
760             {
761                 CPLError( CE_Warning, CPLE_AppDefined,
762                           "Failed to lookup datum code %d. "
763                           "Falling back to use WGS84.",
764                           static_cast<int>(iDatum) );
765                 SetWellKnownGeogCS("WGS84");
766             }
767         }
768         else
769         {
770             CPLError( CE_Warning, CPLE_AppDefined,
771                       "Wrong datum code %d. Supported datums 0--%d only.  "
772                       "Setting WGS84 as a fallback.",
773                       static_cast<int>(iDatum), NUMBER_OF_ELLIPSOIDS );
774             SetWellKnownGeogCS( "WGS84" );
775         }
776 
777         CPLFree( pszName );
778     }
779 
780 /* -------------------------------------------------------------------- */
781 /*      Grid units translation                                          */
782 /* -------------------------------------------------------------------- */
783     if( IsLocal() || IsProjected() )
784         SetLinearUnits( SRS_UL_METER, 1.0 );
785 
786     return OGRERR_NONE;
787 }
788 
789 /************************************************************************/
790 /*                          OSRExportToUSGS()                           */
791 /************************************************************************/
792 /**
793  * \brief Export coordinate system in USGS GCTP projection definition.
794  *
795  * This function is the same as OGRSpatialReference::exportToUSGS().
796  */
797 
OSRExportToUSGS(OGRSpatialReferenceH hSRS,long * piProjSys,long * piZone,double ** ppadfPrjParams,long * piDatum)798 OGRErr OSRExportToUSGS( OGRSpatialReferenceH hSRS,
799                         long *piProjSys, long *piZone,
800                         double **ppadfPrjParams, long *piDatum )
801 
802 {
803     VALIDATE_POINTER1( hSRS, "OSRExportToUSGS", OGRERR_FAILURE );
804 
805     *ppadfPrjParams = nullptr;
806 
807     return OGRSpatialReference::FromHandle(hSRS)->exportToUSGS( piProjSys, piZone,
808                                                          ppadfPrjParams,
809                                                          piDatum );
810 }
811 
812 /************************************************************************/
813 /*                           exportToUSGS()                             */
814 /************************************************************************/
815 
816 /**
817  * \brief Export coordinate system in USGS GCTP projection definition.
818  *
819  * This method is the equivalent of the C function OSRExportToUSGS().
820  *
821  * @param piProjSys Pointer to variable, where the projection system code will
822  * be returned.
823  *
824  * @param piZone Pointer to variable, where the zone for UTM and State Plane
825  * projection systems will be returned.
826  *
827  * @param ppadfPrjParams Pointer to which dynamically allocated array of
828  * 15 projection parameters will be assigned. See importFromUSGS() for
829  * the list of parameters. Caller responsible to free this array.
830  *
831  * @param piDatum Pointer to variable, where the datum code will
832  * be returned.
833  *
834  * @return OGRERR_NONE on success or an error code on failure.
835  */
836 
exportToUSGS(long * piProjSys,long * piZone,double ** ppadfPrjParams,long * piDatum) const837 OGRErr OGRSpatialReference::exportToUSGS( long *piProjSys, long *piZone,
838                                           double **ppadfPrjParams,
839                                           long *piDatum ) const
840 
841 {
842     const char *pszProjection = GetAttrValue("PROJECTION");
843 
844 /* -------------------------------------------------------------------- */
845 /*      Fill all projection parameters with zero.                       */
846 /* -------------------------------------------------------------------- */
847     *ppadfPrjParams = static_cast<double *>(CPLMalloc(15 * sizeof(double)));
848     for( int i = 0; i < 15; i++ )
849         (*ppadfPrjParams)[i] = 0.0;
850 
851     *piZone = 0L;
852 
853 /* ==================================================================== */
854 /*      Handle the projection definition.                               */
855 /* ==================================================================== */
856     if( IsLocal() )
857         *piProjSys = GEO;
858 
859     else if( pszProjection == nullptr )
860     {
861 #ifdef DEBUG
862         CPLDebug( "OSR_USGS",
863                   "Empty projection definition, considered as Geographic" );
864 #endif
865         *piProjSys = GEO;
866     }
867 
868     else if( EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
869     {
870         *piProjSys = ALBERS;
871         (*ppadfPrjParams)[2] = CPLDecToPackedDMS(
872             GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 ) );
873         (*ppadfPrjParams)[3] = CPLDecToPackedDMS(
874             GetNormProjParm( SRS_PP_STANDARD_PARALLEL_2, 0.0 ) );
875         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
876             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
877         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
878             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
879         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
880         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
881     }
882 
883     else if( EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) )
884     {
885         *piProjSys = LAMCC;
886         (*ppadfPrjParams)[2] = CPLDecToPackedDMS(
887             GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 ) );
888         (*ppadfPrjParams)[3] = CPLDecToPackedDMS(
889             GetNormProjParm( SRS_PP_STANDARD_PARALLEL_2, 0.0 ) );
890         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
891             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
892         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
893             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
894         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
895         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
896     }
897 
898     else if( EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) )
899     {
900         *piProjSys = MERCAT;
901         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
902             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
903         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
904             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
905         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
906         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
907     }
908 
909     else if( EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC) )
910     {
911         *piProjSys = PS;
912         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
913             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
914         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
915             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
916         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
917         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
918     }
919 
920     else if( EQUAL(pszProjection, SRS_PT_POLYCONIC) )
921     {
922         *piProjSys = POLYC;
923         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
924             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
925         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
926             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
927         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
928         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
929     }
930 
931     else if( EQUAL(pszProjection, SRS_PT_EQUIDISTANT_CONIC) )
932     {
933         *piProjSys = EQUIDC;
934         (*ppadfPrjParams)[2] = CPLDecToPackedDMS(
935             GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 ) );
936         (*ppadfPrjParams)[3] = CPLDecToPackedDMS(
937             GetNormProjParm( SRS_PP_STANDARD_PARALLEL_2, 0.0 ) );
938         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
939             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
940         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
941             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
942         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
943         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
944         (*ppadfPrjParams)[8] = 1.0;
945     }
946 
947     else if( EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR) )
948     {
949         int bNorth;
950 
951         *piZone = GetUTMZone( &bNorth );
952 
953         if( *piZone != 0 )
954         {
955             *piProjSys = UTM;
956             if( !bNorth )
957                 *piZone = - *piZone;
958         }
959         else
960         {
961             *piProjSys = TM;
962             (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
963             (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
964                 GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
965             (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
966                 GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
967             (*ppadfPrjParams)[6] =
968                 GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
969             (*ppadfPrjParams)[7] =
970                 GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
971         }
972     }
973 
974     else if( EQUAL(pszProjection, SRS_PT_STEREOGRAPHIC) )
975     {
976         *piProjSys = STEREO;
977         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
978             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
979         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
980             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
981         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
982         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
983     }
984 
985     else if( EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
986     {
987         *piProjSys = LAMAZ;
988         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
989             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
990         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
991             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
992         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
993         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
994     }
995 
996     else if( EQUAL(pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT) )
997     {
998         *piProjSys = AZMEQD;
999         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1000             GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 ) );
1001         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
1002             GetNormProjParm( SRS_PP_LATITUDE_OF_CENTER, 0.0 ) );
1003         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1004         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1005     }
1006 
1007     else if( EQUAL(pszProjection, SRS_PT_GNOMONIC) )
1008     {
1009         *piProjSys = GNOMON;
1010         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1011             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
1012         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
1013             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
1014         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1015         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1016     }
1017 
1018     else if( EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC) )
1019     {
1020         *piProjSys = ORTHO;
1021         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1022             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
1023         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
1024             GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) );
1025         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1026         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1027     }
1028 
1029     else if( EQUAL(pszProjection, SRS_PT_SINUSOIDAL) )
1030     {
1031         *piProjSys = SNSOID;
1032         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1033             GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 ) );
1034         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1035         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1036     }
1037 
1038     else if( EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR) )
1039     {
1040         *piProjSys = EQRECT;
1041         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1042             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
1043         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
1044             GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 ) );
1045         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1046         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1047     }
1048 
1049     else if( EQUAL(pszProjection, SRS_PT_MILLER_CYLINDRICAL) )
1050     {
1051         *piProjSys = MILLER;
1052         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1053             GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 ) );
1054         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
1055             GetNormProjParm( SRS_PP_LATITUDE_OF_CENTER, 0.0 ) );
1056         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1057         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1058     }
1059 
1060     else if( EQUAL(pszProjection, SRS_PT_VANDERGRINTEN) )
1061     {
1062         *piProjSys = VGRINT;
1063         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1064             GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 ) );
1065         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1066         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1067     }
1068 
1069     else if( EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR) )
1070     {
1071         *piProjSys = HOM;
1072         (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
1073         (*ppadfPrjParams)[3] = CPLDecToPackedDMS(
1074             GetNormProjParm( SRS_PP_AZIMUTH, 0.0 ) );
1075         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1076             GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 ) );
1077         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
1078             GetNormProjParm( SRS_PP_LATITUDE_OF_CENTER, 0.0 ) );
1079         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1080         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1081         (*ppadfPrjParams)[12] = 1.0;
1082     }
1083 
1084     else if( EQUAL(pszProjection,
1085                    SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN) )
1086     {
1087         *piProjSys = HOM;
1088         (*ppadfPrjParams)[2] = GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 );
1089         (*ppadfPrjParams)[5] = CPLDecToPackedDMS(
1090             GetNormProjParm( SRS_PP_LATITUDE_OF_CENTER, 0.0 ) );
1091         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1092         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1093         (*ppadfPrjParams)[8] = CPLDecToPackedDMS(
1094             GetNormProjParm( SRS_PP_LONGITUDE_OF_POINT_1, 0.0 ) );
1095         (*ppadfPrjParams)[9] = CPLDecToPackedDMS(
1096             GetNormProjParm( SRS_PP_LATITUDE_OF_POINT_1, 0.0 ) );
1097         (*ppadfPrjParams)[10] = CPLDecToPackedDMS(
1098             GetNormProjParm( SRS_PP_LONGITUDE_OF_POINT_2, 0.0 ) );
1099         (*ppadfPrjParams)[11] = CPLDecToPackedDMS(
1100             GetNormProjParm( SRS_PP_LATITUDE_OF_POINT_2, 0.0 ) );
1101         (*ppadfPrjParams)[12] = 0.0;
1102     }
1103 
1104     else if( EQUAL(pszProjection, SRS_PT_ROBINSON) )
1105     {
1106         *piProjSys = ROBIN;
1107         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1108             GetNormProjParm( SRS_PP_LONGITUDE_OF_CENTER, 0.0 ) );
1109         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1110         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1111     }
1112 
1113     else if( EQUAL(pszProjection, SRS_PT_MOLLWEIDE) )
1114     {
1115         *piProjSys = MOLL;
1116         (*ppadfPrjParams)[4] = CPLDecToPackedDMS(
1117             GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
1118         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1119         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1120     }
1121 
1122     else if( EQUAL(pszProjection, SRS_PT_WAGNER_IV) )
1123     {
1124         *piProjSys = WAGIV;
1125         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1126         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1127     }
1128 
1129     else if( EQUAL(pszProjection, SRS_PT_WAGNER_VII) )
1130     {
1131         *piProjSys = WAGVII;
1132         (*ppadfPrjParams)[6] = GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 );
1133         (*ppadfPrjParams)[7] = GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 );
1134     }
1135     // Projection unsupported by GCTP.
1136     else
1137     {
1138         CPLDebug( "OSR_USGS",
1139                   "Projection \"%s\" unsupported by USGS GCTP. "
1140                   "Geographic system will be used.", pszProjection );
1141         *piProjSys = GEO;
1142     }
1143 
1144 /* -------------------------------------------------------------------- */
1145 /*      Translate the datum.                                            */
1146 /* -------------------------------------------------------------------- */
1147     const char *pszDatum = GetAttrValue( "DATUM" );
1148 
1149     if( pszDatum )
1150     {
1151         if( EQUAL( pszDatum, SRS_DN_NAD27 ) )
1152         {
1153             *piDatum = CLARKE1866;
1154         }
1155         else if( EQUAL( pszDatum, SRS_DN_NAD83 ) )
1156         {
1157             *piDatum = GRS1980;
1158         }
1159         else if( EQUAL( pszDatum, SRS_DN_WGS84 ) )
1160         {
1161             *piDatum = WGS84;
1162         }
1163         // If not found well known datum, translate ellipsoid.
1164         else
1165         {
1166             const double dfSemiMajor = GetSemiMajor();
1167             const double dfInvFlattening = GetInvFlattening();
1168 
1169 #ifdef DEBUG
1170             CPLDebug( "OSR_USGS",
1171                       "Datum \"%s\" unsupported by USGS GCTP. "
1172                       "Try to translate ellipsoid definition.", pszDatum );
1173 #endif
1174 
1175             int i = 0;  // Used after for.
1176             for( ; i < NUMBER_OF_ELLIPSOIDS; i++ )
1177             {
1178                 double dfSM = 0.0;
1179                 double dfIF = 0.0;
1180 
1181                 if( OSRGetEllipsoidInfo( aoEllips[i], nullptr,
1182                                          &dfSM, &dfIF ) == OGRERR_NONE
1183                     && CPLIsEqual( dfSemiMajor, dfSM )
1184                     && CPLIsEqual( dfInvFlattening, dfIF ) )
1185                 {
1186                     *piDatum = i;
1187                     break;
1188                 }
1189             }
1190 
1191             if( i == NUMBER_OF_ELLIPSOIDS )  // Didn't found matches; set
1192             {                                // custom ellipsoid parameters.
1193 #ifdef DEBUG
1194                 CPLDebug( "OSR_USGS",
1195                           "Ellipsoid \"%s\" unsupported by USGS GCTP. "
1196                           "Custom ellipsoid definition will be used.",
1197                           pszDatum );
1198 #endif
1199                 *piDatum = -1;
1200                 (*ppadfPrjParams)[0] = dfSemiMajor;
1201                 if( std::abs( dfInvFlattening ) < 0.000000000001 )
1202                 {
1203                     (*ppadfPrjParams)[1] = dfSemiMajor;
1204                 }
1205                 else
1206                 {
1207                     (*ppadfPrjParams)[1] =
1208                         dfSemiMajor * (1.0 - 1.0/dfInvFlattening);
1209                 }
1210             }
1211         }
1212     }
1213     else
1214     {
1215         *piDatum = -1;
1216     }
1217 
1218     return OGRERR_NONE;
1219 }
1220