LCOV - code coverage report
Current view: top level - ogr - ogr_srs_esri.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 190 329 57.8 %
Date: 2024-11-21 22:18:42 Functions: 5 7 71.4 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  OpenGIS Simple Features Reference Implementation
       4             :  * Purpose:  OGRSpatialReference translation to/from ESRI .prj definitions.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2000, Frank Warmerdam
       9             :  * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
      10             :  * Copyright (c) 2013, Kyle Shannon <kyle at pobox dot com>
      11             :  *
      12             :  * SPDX-License-Identifier: MIT
      13             :  ****************************************************************************/
      14             : 
      15             : #include "cpl_port.h"
      16             : #include "ogr_spatialref.h"
      17             : #include "ogr_srs_esri_names.h"
      18             : 
      19             : #include <cmath>
      20             : #include <climits>
      21             : #include <cstddef>
      22             : #include <cstdio>
      23             : #include <cstdlib>
      24             : #include <cstring>
      25             : #include <algorithm>
      26             : #include <limits>
      27             : #include <string>
      28             : 
      29             : #include "cpl_conv.h"
      30             : #include "cpl_csv.h"
      31             : #include "cpl_error.h"
      32             : #include "cpl_multiproc.h"
      33             : #include "cpl_string.h"
      34             : #include "cpl_vsi.h"
      35             : #include "ogr_core.h"
      36             : #include "ogr_p.h"
      37             : #include "ogr_srs_api.h"
      38             : 
      39             : /* -------------------------------------------------------------------- */
      40             : /*      Table relating USGS and ESRI state plane zones.                 */
      41             : /* -------------------------------------------------------------------- */
      42             : constexpr int anUsgsEsriZones[] = {
      43             :     101,  3101, 102,  3126, 201,  3151, 202,  3176, 203,  3201, 301,  3226,
      44             :     302,  3251, 401,  3276, 402,  3301, 403,  3326, 404,  3351, 405,  3376,
      45             :     406,  3401, 407,  3426, 501,  3451, 502,  3476, 503,  3501, 600,  3526,
      46             :     700,  3551, 901,  3601, 902,  3626, 903,  3576, 1001, 3651, 1002, 3676,
      47             :     1101, 3701, 1102, 3726, 1103, 3751, 1201, 3776, 1202, 3801, 1301, 3826,
      48             :     1302, 3851, 1401, 3876, 1402, 3901, 1501, 3926, 1502, 3951, 1601, 3976,
      49             :     1602, 4001, 1701, 4026, 1702, 4051, 1703, 6426, 1801, 4076, 1802, 4101,
      50             :     1900, 4126, 2001, 4151, 2002, 4176, 2101, 4201, 2102, 4226, 2103, 4251,
      51             :     2111, 6351, 2112, 6376, 2113, 6401, 2201, 4276, 2202, 4301, 2203, 4326,
      52             :     2301, 4351, 2302, 4376, 2401, 4401, 2402, 4426, 2403, 4451, 2500, 0,
      53             :     2501, 4476, 2502, 4501, 2503, 4526, 2600, 0,    2601, 4551, 2602, 4576,
      54             :     2701, 4601, 2702, 4626, 2703, 4651, 2800, 4676, 2900, 4701, 3001, 4726,
      55             :     3002, 4751, 3003, 4776, 3101, 4801, 3102, 4826, 3103, 4851, 3104, 4876,
      56             :     3200, 4901, 3301, 4926, 3302, 4951, 3401, 4976, 3402, 5001, 3501, 5026,
      57             :     3502, 5051, 3601, 5076, 3602, 5101, 3701, 5126, 3702, 5151, 3800, 5176,
      58             :     3900, 0,    3901, 5201, 3902, 5226, 4001, 5251, 4002, 5276, 4100, 5301,
      59             :     4201, 5326, 4202, 5351, 4203, 5376, 4204, 5401, 4205, 5426, 4301, 5451,
      60             :     4302, 5476, 4303, 5501, 4400, 5526, 4501, 5551, 4502, 5576, 4601, 5601,
      61             :     4602, 5626, 4701, 5651, 4702, 5676, 4801, 5701, 4802, 5726, 4803, 5751,
      62             :     4901, 5776, 4902, 5801, 4903, 5826, 4904, 5851, 5001, 6101, 5002, 6126,
      63             :     5003, 6151, 5004, 6176, 5005, 6201, 5006, 6226, 5007, 6251, 5008, 6276,
      64             :     5009, 6301, 5010, 6326, 5101, 5876, 5102, 5901, 5103, 5926, 5104, 5951,
      65             :     5105, 5976, 5201, 6001, 5200, 6026, 5200, 6076, 5201, 6051, 5202, 6051,
      66             :     5300, 0,    5400, 0};
      67             : 
      68             : /************************************************************************/
      69             : /*                           ESRIToUSGSZone()                           */
      70             : /*                                                                      */
      71             : /*      Convert ESRI style state plane zones to USGS style state        */
      72             : /*      plane zones.                                                    */
      73             : /************************************************************************/
      74             : 
      75           0 : static int ESRIToUSGSZone(int nESRIZone)
      76             : 
      77             : {
      78             :     // anUsgsEsriZones is a series of ints where 2 consecutive integers
      79             :     // are used to map from USGS to ESRI state plane zones.
      80             :     // TODO(schwehr): Would be better as a std::map.
      81           0 :     const int nPairs = sizeof(anUsgsEsriZones) / (2 * sizeof(int));
      82             : 
      83           0 :     for (int i = 0; i < nPairs; i++)
      84             :     {
      85           0 :         if (anUsgsEsriZones[i * 2 + 1] == nESRIZone)
      86           0 :             return anUsgsEsriZones[i * 2];
      87             :     }
      88             : 
      89           0 :     return 0;
      90             : }
      91             : 
      92             : /************************************************************************/
      93             : /*                         OSRImportFromESRI()                          */
      94             : /************************************************************************/
      95             : 
      96             : /**
      97             :  * \brief Import coordinate system from ESRI .prj format(s).
      98             :  *
      99             :  * This function is the same as the C++ method
     100             :  * OGRSpatialReference::importFromESRI().
     101             :  */
     102          11 : OGRErr OSRImportFromESRI(OGRSpatialReferenceH hSRS, char **papszPrj)
     103             : 
     104             : {
     105          11 :     VALIDATE_POINTER1(hSRS, "OSRImportFromESRI", OGRERR_FAILURE);
     106             : 
     107          11 :     return reinterpret_cast<OGRSpatialReference *>(hSRS)->importFromESRI(
     108          11 :         papszPrj);
     109             : }
     110             : 
     111             : /************************************************************************/
     112             : /*                              OSR_GDV()                               */
     113             : /*                                                                      */
     114             : /*      Fetch a particular parameter out of the parameter list, or      */
     115             : /*      the indicated default if it isn't available.  This is a         */
     116             : /*      helper function for importFromESRI().                           */
     117             : /************************************************************************/
     118             : 
     119         134 : static double OSR_GDV(char **papszNV, const char *pszField,
     120             :                       double dfDefaultValue)
     121             : 
     122             : {
     123         134 :     if (papszNV == nullptr || papszNV[0] == nullptr)
     124           0 :         return dfDefaultValue;
     125             : 
     126         134 :     if (STARTS_WITH_CI(pszField, "PARAM_"))
     127             :     {
     128         102 :         int iLine = 0;  // Used after for loop.
     129         790 :         for (; papszNV[iLine] != nullptr &&
     130         790 :                !STARTS_WITH_CI(papszNV[iLine], "Paramet");
     131             :              iLine++)
     132             :         {
     133             :         }
     134             : 
     135         441 :         for (int nOffset = atoi(pszField + 6);
     136         441 :              papszNV[iLine] != nullptr && nOffset > 0; iLine++)
     137             :         {
     138         339 :             if (strlen(papszNV[iLine]) > 0)
     139         339 :                 nOffset--;
     140             :         }
     141             : 
     142         102 :         while (papszNV[iLine] != nullptr && strlen(papszNV[iLine]) == 0)
     143           0 :             iLine++;
     144             : 
     145         102 :         if (papszNV[iLine] != nullptr)
     146             :         {
     147         102 :             char *const pszLine = papszNV[iLine];
     148             : 
     149             :             // Trim comments.
     150        7027 :             for (int i = 0; pszLine[i] != '\0'; i++)
     151             :             {
     152        6925 :                 if (pszLine[i] == '/' && pszLine[i + 1] == '*')
     153         101 :                     pszLine[i] = '\0';
     154             :             }
     155             : 
     156         102 :             double dfValue = 0.0;
     157         102 :             char **papszTokens = CSLTokenizeString(papszNV[iLine]);
     158         102 :             if (CSLCount(papszTokens) == 3)
     159             :             {
     160             :                 // http://agdcftp1.wr.usgs.gov/pub/projects/lcc/akcan_lcc/akcan.tar.gz
     161             :                 // contains weird values for the second. Ignore it and
     162             :                 // the result looks correct.
     163          62 :                 double dfSecond = CPLAtof(papszTokens[2]);
     164          62 :                 if (dfSecond < 0.0 || dfSecond >= 60.0)
     165           0 :                     dfSecond = 0.0;
     166             : 
     167          62 :                 dfValue = std::abs(CPLAtof(papszTokens[0])) +
     168          62 :                           CPLAtof(papszTokens[1]) / 60.0 + dfSecond / 3600.0;
     169             : 
     170          62 :                 if (CPLAtof(papszTokens[0]) < 0.0)
     171          15 :                     dfValue *= -1;
     172             :             }
     173          40 :             else if (CSLCount(papszTokens) > 0)
     174             :             {
     175          40 :                 dfValue = CPLAtof(papszTokens[0]);
     176             :             }
     177             :             else
     178             :             {
     179           0 :                 dfValue = dfDefaultValue;
     180             :             }
     181             : 
     182         102 :             CSLDestroy(papszTokens);
     183             : 
     184         102 :             return dfValue;
     185             :         }
     186             : 
     187           0 :         return dfDefaultValue;
     188             :     }
     189             : 
     190          32 :     int iLine = 0;  // Used after for loop.
     191         156 :     for (; papszNV[iLine] != nullptr &&
     192         152 :            !EQUALN(papszNV[iLine], pszField, strlen(pszField));
     193             :          iLine++)
     194             :     {
     195             :     }
     196             : 
     197          32 :     if (papszNV[iLine] == nullptr)
     198           4 :         return dfDefaultValue;
     199             : 
     200          28 :     return CPLAtof(papszNV[iLine] + strlen(pszField));
     201             : }
     202             : 
     203             : /************************************************************************/
     204             : /*                              OSR_GDS()                               */
     205             : /************************************************************************/
     206             : 
     207         117 : static CPLString OSR_GDS(char **papszNV, const char *pszField,
     208             :                          const char *pszDefaultValue)
     209             : 
     210             : {
     211         117 :     if (papszNV == nullptr || papszNV[0] == nullptr)
     212           0 :         return pszDefaultValue;
     213             : 
     214         117 :     int iLine = 0;  // Used after for loop.
     215         370 :     for (; papszNV[iLine] != nullptr &&
     216         359 :            !EQUALN(papszNV[iLine], pszField, strlen(pszField));
     217             :          iLine++)
     218             :     {
     219             :     }
     220             : 
     221         117 :     if (papszNV[iLine] == nullptr)
     222          11 :         return pszDefaultValue;
     223             : 
     224         106 :     char **papszTokens = CSLTokenizeString(papszNV[iLine]);
     225             : 
     226             :     CPLString osResult =
     227         212 :         CSLCount(papszTokens) > 1 ? papszTokens[1] : pszDefaultValue;
     228             : 
     229         106 :     CSLDestroy(papszTokens);
     230         106 :     return osResult;
     231             : }
     232             : 
     233             : /************************************************************************/
     234             : /*                          importFromESRI()                            */
     235             : /************************************************************************/
     236             : 
     237             : /**
     238             :  * \brief Import coordinate system from ESRI .prj format(s).
     239             :  *
     240             :  * This function will read the text loaded from an ESRI .prj file, and
     241             :  * translate it into an OGRSpatialReference definition.  This should support
     242             :  * many (but by no means all) old style (Arc/Info 7.x) .prj files, as well
     243             :  * as the newer pseudo-OGC WKT .prj files.  Note that new style .prj files
     244             :  * are in OGC WKT format, but require some manipulation to correct datum
     245             :  * names, and units on some projection parameters.  This is addressed within
     246             :  * importFromESRI() by an automatic call to morphFromESRI().
     247             :  *
     248             :  * Currently only GEOGRAPHIC, UTM, STATEPLANE, GREATBRITIAN_GRID, ALBERS,
     249             :  * EQUIDISTANT_CONIC, TRANSVERSE (mercator), POLAR, LAMBERT (Conic Conformal),
     250             :  * LAMBERT_AZIMUTHAL, MERCATOR and POLYCONIC
     251             :  * projections are supported from old style files.
     252             :  *
     253             :  * At this time there is no equivalent exportToESRI() method.  Writing old
     254             :  * style .prj files is not supported by OGRSpatialReference. However the
     255             :  * morphToESRI() and exportToWkt() methods can be used to generate output
     256             :  * suitable to write to new style (Arc 8) .prj files.
     257             :  *
     258             :  * This function is the equivalent of the C function OSRImportFromESRI().
     259             :  *
     260             :  * @param papszPrj NULL terminated list of strings containing the definition.
     261             :  *
     262             :  * @return OGRERR_NONE on success or an error code in case of failure.
     263             :  */
     264             : 
     265         943 : OGRErr OGRSpatialReference::importFromESRI(char **papszPrj)
     266             : 
     267             : {
     268         943 :     if (papszPrj == nullptr || papszPrj[0] == nullptr)
     269          29 :         return OGRERR_CORRUPT_DATA;
     270             : 
     271             :     /* -------------------------------------------------------------------- */
     272             :     /*      ArcGIS and related products now use a variant of Well Known     */
     273             :     /*      Text.  Try to recognize this and ingest it.  WKT is usually     */
     274             :     /*      all on one line, but we will accept multi-line formats and      */
     275             :     /*      concatenate.                                                    */
     276             :     /* -------------------------------------------------------------------- */
     277         914 :     if (STARTS_WITH_CI(papszPrj[0], "GEOGCS") ||
     278         563 :         STARTS_WITH_CI(papszPrj[0], "PROJCS") ||
     279          43 :         STARTS_WITH_CI(papszPrj[0], "LOCAL_CS")
     280             :         // Also accept COMPD_CS, even if it is unclear that it is valid
     281             :         // traditional ESRI WKT. But people might use such PRJ file
     282             :         // See https://github.com/OSGeo/gdal/issues/1881
     283          43 :         || STARTS_WITH_CI(papszPrj[0], "COMPD_CS"))
     284             :     {
     285        1742 :         std::string osWKT(papszPrj[0]);
     286         873 :         for (int i = 1; papszPrj[i] != nullptr; i++)
     287             :         {
     288           2 :             osWKT += papszPrj[i];
     289             :         }
     290         871 :         return importFromWkt(osWKT.c_str());
     291             :     }
     292             : 
     293             :     /* -------------------------------------------------------------------- */
     294             :     /*      Operate on the basis of the projection name.                    */
     295             :     /* -------------------------------------------------------------------- */
     296          86 :     CPLString osProj = OSR_GDS(papszPrj, "Projection", "");
     297          43 :     bool bDatumApplied = false;
     298          43 :     bool bHasRadiusOfSphereOfReference = false;
     299          43 :     double dfRadiusOfSphereOfReference = 0.0;
     300             : 
     301          43 :     if (EQUAL(osProj, ""))
     302             :     {
     303           5 :         CPLDebug("OGR_ESRI", "Can't find Projection");
     304           5 :         return OGRERR_CORRUPT_DATA;
     305             :     }
     306          38 :     else if (EQUAL(osProj, "GEOGRAPHIC"))
     307             :     {
     308             :         // Nothing to do.
     309             :     }
     310          35 :     else if (EQUAL(osProj, "utm"))
     311             :     {
     312          12 :         const double dfOsrGdv = OSR_GDV(papszPrj, "zone", 0.0);
     313          12 :         if (dfOsrGdv > 0 && dfOsrGdv < 61)
     314             :         {
     315          12 :             const double dfYShift = OSR_GDV(papszPrj, "Yshift", 0.0);
     316             : 
     317          12 :             SetUTM(static_cast<int>(dfOsrGdv), dfYShift == 0.0);
     318             :         }
     319             :         else
     320             :         {
     321           0 :             const double dfCentralMeridian = OSR_GDV(papszPrj, "PARAM_1", 0.0);
     322           0 :             const double dfRefLat = OSR_GDV(papszPrj, "PARAM_2", 0.0);
     323           0 :             if (dfCentralMeridian >= -180.0 && dfCentralMeridian <= 180.0)
     324             :             {
     325           0 :                 const int nZone = static_cast<int>(
     326           0 :                     (dfCentralMeridian + 183.0) / 6.0 + 0.0000001);
     327           0 :                 SetUTM(nZone, dfRefLat >= 0.0);
     328             :             }
     329             :         }
     330             :     }
     331          23 :     else if (EQUAL(osProj, "STATEPLANE"))
     332             :     {
     333           4 :         const double dfZone = OSR_GDV(papszPrj, "zone", 0.0);
     334             : 
     335           4 :         if (dfZone < std::numeric_limits<int>::min() ||
     336           4 :             dfZone > std::numeric_limits<int>::max() || std::isnan(dfZone))
     337             :         {
     338           0 :             CPLError(CE_Failure, CPLE_AppDefined, "zone out of range: %f",
     339             :                      dfZone);
     340           0 :             return OGRERR_CORRUPT_DATA;
     341             :         }
     342             : 
     343           4 :         int nZone = static_cast<int>(dfZone);
     344             : 
     345           4 :         if (nZone != 0)
     346           0 :             nZone = ESRIToUSGSZone(nZone);
     347             :         else
     348             :         {
     349           4 :             const double dfFipszone = OSR_GDV(papszPrj, "fipszone", 0.0);
     350             : 
     351           4 :             if (dfFipszone < std::numeric_limits<int>::min() ||
     352           8 :                 dfFipszone > std::numeric_limits<int>::max() ||
     353           4 :                 std::isnan(dfFipszone))
     354             :             {
     355           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     356             :                          "fipszone out of range: %f", dfFipszone);
     357           0 :                 return OGRERR_CORRUPT_DATA;
     358             :             }
     359             : 
     360           4 :             nZone = static_cast<int>(dfFipszone);
     361             :         }
     362             : 
     363           4 :         if (nZone != 0)
     364             :         {
     365           4 :             bDatumApplied = true;
     366           4 :             if (EQUAL(OSR_GDS(papszPrj, "Datum", "NAD83"), "NAD27"))
     367           0 :                 SetStatePlane(nZone, FALSE);
     368             :             else
     369           4 :                 SetStatePlane(nZone, TRUE);
     370             :         }
     371             :     }
     372          38 :     else if (EQUAL(osProj, "GREATBRITIAN_GRID") ||
     373          19 :              EQUAL(osProj, "GREATBRITAIN_GRID"))
     374             :     {
     375           0 :         const char *pszWkt =
     376             :             "PROJCS[\"OSGB 1936 / British National Grid\","
     377             :             "GEOGCS[\"OSGB 1936\",DATUM[\"OSGB_1936\","
     378             :             "SPHEROID[\"Airy 1830\",6377563.396,299.3249646]],"
     379             :             "PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],"
     380             :             "PROJECTION[\"Transverse_Mercator\"],"
     381             :             "PARAMETER[\"latitude_of_origin\",49],"
     382             :             "PARAMETER[\"central_meridian\",-2],"
     383             :             "PARAMETER[\"scale_factor\",0.999601272],"
     384             :             "PARAMETER[\"false_easting\",400000],"
     385             :             "PARAMETER[\"false_northing\",-100000],UNIT[\"metre\",1]]";
     386             : 
     387           0 :         bDatumApplied = true;
     388           0 :         importFromWkt(pszWkt);
     389             :     }
     390          19 :     else if (EQUAL(osProj, "ALBERS"))
     391             :     {
     392          13 :         SetACEA(OSR_GDV(papszPrj, "PARAM_1", 0.0),
     393             :                 OSR_GDV(papszPrj, "PARAM_2", 0.0),
     394             :                 OSR_GDV(papszPrj, "PARAM_4", 0.0),
     395             :                 OSR_GDV(papszPrj, "PARAM_3", 0.0),
     396             :                 OSR_GDV(papszPrj, "PARAM_5", 0.0),
     397             :                 OSR_GDV(papszPrj, "PARAM_6", 0.0));
     398             :     }
     399           6 :     else if (EQUAL(osProj, "LAMBERT"))
     400             :     {
     401           0 :         SetLCC(OSR_GDV(papszPrj, "PARAM_1", 0.0),
     402             :                OSR_GDV(papszPrj, "PARAM_2", 0.0),
     403             :                OSR_GDV(papszPrj, "PARAM_4", 0.0),
     404             :                OSR_GDV(papszPrj, "PARAM_3", 0.0),
     405             :                OSR_GDV(papszPrj, "PARAM_5", 0.0),
     406             :                OSR_GDV(papszPrj, "PARAM_6", 0.0));
     407             :     }
     408           6 :     else if (EQUAL(osProj, "LAMBERT_AZIMUTHAL"))
     409             :     {
     410          27 :         for (int iLine = 0; papszPrj[iLine] != nullptr; iLine++)
     411             :         {
     412          26 :             if (strstr(papszPrj[iLine], "radius of the sphere of reference"))
     413             :             {
     414           2 :                 bHasRadiusOfSphereOfReference = true;
     415           2 :                 break;
     416             :             }
     417             :         }
     418           3 :         if (bHasRadiusOfSphereOfReference)
     419             :         {
     420             :             // Cf "Workstation" variation of
     421             :             // https://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Lambert_Azimuthal_Equal_Area
     422             :             // that is supposed to be applied to spheres only.
     423             :             // It is also documented at page 71 of
     424             :             // https://kartoweb.itc.nl/geometrics/Map%20projections/Understanding%20Map%20Projections.pdf
     425             :             // ("Understanding Map Projections", Melita Kennedy, ArcInfo 8)
     426             :             // We don't particularly enforce the restriction to spheres, so if
     427             :             // the "Parameters" line had non-spherical axis dimensions, they
     428             :             // would be used.
     429             :             // EPSG has a EPSG:1027 projection method "Lambert Azimuthal Equal Area (Spherical)"
     430             :             // that uses the radius of the authalic sphere, for non-spherical
     431             :             // ellipsoids, but it is not obvious that it would be appropriate here.
     432             :             // Examples:
     433             :             // - https://community.esri.com/t5/data-management-questions/problem-when-projecting-a-raster-file/td-p/400814/page/2
     434             :             // - https://lists.osgeo.org/pipermail/gdal-dev/2024-May/058990.html
     435           2 :             dfRadiusOfSphereOfReference = OSR_GDV(papszPrj, "PARAM_1", 0.0);
     436           2 :             SetLAEA(OSR_GDV(papszPrj, "PARAM_3", 0.0),
     437             :                     OSR_GDV(papszPrj, "PARAM_2", 0.0),
     438             :                     OSR_GDV(papszPrj, "PARAM_4", 0.0),
     439             :                     OSR_GDV(papszPrj, "PARAM_5", 0.0));
     440             :         }
     441             :         else
     442             :         {
     443             :             // Example: https://trac.osgeo.org/gdal/ticket/4302
     444           1 :             SetLAEA(OSR_GDV(papszPrj, "PARAM_2", 0.0),
     445             :                     OSR_GDV(papszPrj, "PARAM_1", 0.0),
     446             :                     OSR_GDV(papszPrj, "PARAM_3", 0.0),
     447             :                     OSR_GDV(papszPrj, "PARAM_4", 0.0));
     448             :         }
     449             :     }
     450           3 :     else if (EQUAL(osProj, "EQUIDISTANT_CONIC"))
     451             :     {
     452           1 :         const double dfStdPCount = OSR_GDV(papszPrj, "PARAM_1", 0.0);
     453             :         // TODO(schwehr): What is a reasonable range for StdPCount?
     454           1 :         if (dfStdPCount < 0 || dfStdPCount > std::numeric_limits<int>::max() ||
     455           0 :             std::isnan(dfStdPCount))
     456             :         {
     457           1 :             CPLError(CE_Failure, CPLE_AppDefined, "StdPCount out of range: %lf",
     458             :                      dfStdPCount);
     459           1 :             return OGRERR_CORRUPT_DATA;
     460             :         }
     461           0 :         const int nStdPCount = static_cast<int>(dfStdPCount);
     462             : 
     463           0 :         if (nStdPCount == 1)
     464             :         {
     465           0 :             SetEC(OSR_GDV(papszPrj, "PARAM_2", 0.0),
     466             :                   OSR_GDV(papszPrj, "PARAM_2", 0.0),
     467             :                   OSR_GDV(papszPrj, "PARAM_4", 0.0),
     468             :                   OSR_GDV(papszPrj, "PARAM_3", 0.0),
     469             :                   OSR_GDV(papszPrj, "PARAM_5", 0.0),
     470             :                   OSR_GDV(papszPrj, "PARAM_6", 0.0));
     471             :         }
     472             :         else
     473             :         {
     474           0 :             SetEC(OSR_GDV(papszPrj, "PARAM_2", 0.0),
     475             :                   OSR_GDV(papszPrj, "PARAM_3", 0.0),
     476             :                   OSR_GDV(papszPrj, "PARAM_5", 0.0),
     477             :                   OSR_GDV(papszPrj, "PARAM_4", 0.0),
     478             :                   OSR_GDV(papszPrj, "PARAM_5", 0.0),
     479             :                   OSR_GDV(papszPrj, "PARAM_7", 0.0));
     480             :         }
     481             :     }
     482           2 :     else if (EQUAL(osProj, "TRANSVERSE"))
     483             :     {
     484           1 :         SetTM(OSR_GDV(papszPrj, "PARAM_3", 0.0),
     485             :               OSR_GDV(papszPrj, "PARAM_2", 0.0),
     486             :               OSR_GDV(papszPrj, "PARAM_1", 0.0),
     487             :               OSR_GDV(papszPrj, "PARAM_4", 0.0),
     488             :               OSR_GDV(papszPrj, "PARAM_5", 0.0));
     489             :     }
     490           1 :     else if (EQUAL(osProj, "POLAR"))
     491             :     {
     492           0 :         SetPS(OSR_GDV(papszPrj, "PARAM_2", 0.0),
     493             :               OSR_GDV(papszPrj, "PARAM_1", 0.0), 1.0,
     494             :               OSR_GDV(papszPrj, "PARAM_3", 0.0),
     495             :               OSR_GDV(papszPrj, "PARAM_4", 0.0));
     496             :     }
     497           1 :     else if (EQUAL(osProj, "MERCATOR"))
     498             :     {
     499           1 :         SetMercator2SP(OSR_GDV(papszPrj, "PARAM_2", 0.0), 0.0,
     500             :                        OSR_GDV(papszPrj, "PARAM_1", 0.0),
     501             :                        OSR_GDV(papszPrj, "PARAM_3", 0.0),
     502             :                        OSR_GDV(papszPrj, "PARAM_4", 0.0));
     503             :     }
     504           0 :     else if (EQUAL(osProj, SRS_PT_MERCATOR_AUXILIARY_SPHERE))
     505             :     {
     506             :         // This is EPSG:3875 Pseudo Mercator. We might as well import it from
     507             :         // the EPSG spec.
     508           0 :         importFromEPSG(3857);
     509           0 :         bDatumApplied = true;
     510             :     }
     511           0 :     else if (EQUAL(osProj, "POLYCONIC"))
     512             :     {
     513           0 :         SetPolyconic(OSR_GDV(papszPrj, "PARAM_2", 0.0),
     514             :                      OSR_GDV(papszPrj, "PARAM_1", 0.0),
     515             :                      OSR_GDV(papszPrj, "PARAM_3", 0.0),
     516             :                      OSR_GDV(papszPrj, "PARAM_4", 0.0));
     517             :     }
     518             :     else
     519             :     {
     520           0 :         CPLDebug("OGR_ESRI", "Unsupported projection: %s", osProj.c_str());
     521           0 :         SetLocalCS(osProj);
     522             :     }
     523             : 
     524             :     /* -------------------------------------------------------------------- */
     525             :     /*      Try to translate the datum/spheroid.                            */
     526             :     /* -------------------------------------------------------------------- */
     527          37 :     if (!IsLocal() && !bDatumApplied)
     528             :     {
     529          66 :         const CPLString osDatum = OSR_GDS(papszPrj, "Datum", "");
     530             : 
     531          61 :         if (EQUAL(osDatum, "NAD27") || EQUAL(osDatum, "NAD83") ||
     532          61 :             EQUAL(osDatum, "WGS84") || EQUAL(osDatum, "WGS72"))
     533             :         {
     534          21 :             SetWellKnownGeogCS(osDatum);
     535             :         }
     536          12 :         else if (EQUAL(osDatum, "EUR") || EQUAL(osDatum, "ED50"))
     537             :         {
     538           0 :             SetWellKnownGeogCS("EPSG:4230");
     539             :         }
     540          12 :         else if (EQUAL(osDatum, "GDA94"))
     541             :         {
     542           9 :             SetWellKnownGeogCS("EPSG:4283");
     543             :         }
     544             :         else
     545             :         {
     546           6 :             CPLString osSpheroid = OSR_GDS(papszPrj, "Spheroid", "");
     547             : 
     548           6 :             if (EQUAL(osSpheroid, "INT1909") ||
     549           3 :                 EQUAL(osSpheroid, "INTERNATIONAL1909"))
     550             :             {
     551           0 :                 OGRSpatialReference oGCS;
     552           0 :                 oGCS.importFromEPSG(4022);
     553           0 :                 CopyGeogCSFrom(&oGCS);
     554             :             }
     555           3 :             else if (EQUAL(osSpheroid, "AIRY"))
     556             :             {
     557           0 :                 OGRSpatialReference oGCS;
     558           0 :                 oGCS.importFromEPSG(4001);
     559           0 :                 CopyGeogCSFrom(&oGCS);
     560             :             }
     561           3 :             else if (EQUAL(osSpheroid, "CLARKE1866"))
     562             :             {
     563           0 :                 OGRSpatialReference oGCS;
     564           0 :                 oGCS.importFromEPSG(4008);
     565           0 :                 CopyGeogCSFrom(&oGCS);
     566             :             }
     567           3 :             else if (EQUAL(osSpheroid, "GRS80"))
     568             :             {
     569           0 :                 OGRSpatialReference oGCS;
     570           0 :                 oGCS.importFromEPSG(4019);
     571           0 :                 CopyGeogCSFrom(&oGCS);
     572             :             }
     573           3 :             else if (EQUAL(osSpheroid, "KRASOVSKY") ||
     574           6 :                      EQUAL(osSpheroid, "KRASSOVSKY") ||
     575           3 :                      EQUAL(osSpheroid, "KRASSOWSKY"))
     576             :             {
     577           0 :                 OGRSpatialReference oGCS;
     578           0 :                 oGCS.importFromEPSG(4024);
     579           0 :                 CopyGeogCSFrom(&oGCS);
     580             :             }
     581           3 :             else if (EQUAL(osSpheroid, "Bessel"))
     582             :             {
     583           0 :                 OGRSpatialReference oGCS;
     584           0 :                 oGCS.importFromEPSG(4004);
     585           0 :                 CopyGeogCSFrom(&oGCS);
     586             :             }
     587             :             else
     588             :             {
     589           3 :                 bool bFoundParameters = false;
     590          14 :                 for (int iLine = 0; papszPrj[iLine] != nullptr; iLine++)
     591             :                 {
     592          14 :                     if (STARTS_WITH_CI(papszPrj[iLine], "Parameters"))
     593             :                     {
     594           6 :                         char **papszTokens = CSLTokenizeString(
     595           3 :                             papszPrj[iLine] + strlen("Parameters"));
     596           3 :                         if (CSLCount(papszTokens) == 2)
     597             :                         {
     598           2 :                             OGRSpatialReference oGCS;
     599           2 :                             const double dfSemiMajor = CPLAtof(papszTokens[0]);
     600           2 :                             const double dfSemiMinor = CPLAtof(papszTokens[1]);
     601             :                             const double dfInvFlattening =
     602           2 :                                 OSRCalcInvFlattening(dfSemiMajor, dfSemiMinor);
     603           2 :                             oGCS.SetGeogCS("unknown", "unknown", "unknown",
     604             :                                            dfSemiMajor, dfInvFlattening);
     605           2 :                             CopyGeogCSFrom(&oGCS);
     606           2 :                             bFoundParameters = true;
     607             :                         }
     608           3 :                         CSLDestroy(papszTokens);
     609           3 :                         break;
     610             :                     }
     611             :                 }
     612           3 :                 if (!bFoundParameters)
     613             :                 {
     614           1 :                     if (bHasRadiusOfSphereOfReference)
     615             :                     {
     616           2 :                         OGRSpatialReference oGCS;
     617           1 :                         oGCS.SetGeogCS("unknown", "unknown", "unknown",
     618             :                                        dfRadiusOfSphereOfReference, 0);
     619           1 :                         CopyGeogCSFrom(&oGCS);
     620             :                     }
     621             :                     else
     622             :                     {
     623             :                         // If unknown, default to WGS84 so there is something there.
     624           0 :                         SetWellKnownGeogCS("WGS84");
     625             :                     }
     626             :                 }
     627             :             }
     628             :         }
     629             :     }
     630             : 
     631             :     /* -------------------------------------------------------------------- */
     632             :     /*      Linear units translation                                        */
     633             :     /* -------------------------------------------------------------------- */
     634          37 :     if (IsLocal() || IsProjected())
     635             :     {
     636          34 :         const double dfOldUnits = GetLinearUnits();
     637          68 :         const CPLString osValue = OSR_GDS(papszPrj, "Units", "");
     638          68 :         CPLString osOldAuth;
     639             :         {
     640          34 :             const char *pszOldAuth = GetAuthorityCode(nullptr);
     641          34 :             if (pszOldAuth)
     642           4 :                 osOldAuth = pszOldAuth;
     643             :         }
     644             : 
     645          34 :         if (EQUAL(osValue, ""))
     646           0 :             SetLinearUnitsAndUpdateParameters(SRS_UL_METER, 1.0);
     647          34 :         else if (EQUAL(osValue, "FEET"))
     648           2 :             SetLinearUnitsAndUpdateParameters(SRS_UL_US_FOOT,
     649             :                                               CPLAtof(SRS_UL_US_FOOT_CONV));
     650          32 :         else if (CPLAtof(osValue) != 0.0)
     651           1 :             SetLinearUnitsAndUpdateParameters("user-defined",
     652           1 :                                               1.0 / CPLAtof(osValue));
     653             :         else
     654          31 :             SetLinearUnitsAndUpdateParameters(osValue, 1.0);
     655             : 
     656             :         // Reinstall authority if linear units value has not changed (bug #1697)
     657          34 :         const double dfNewUnits = GetLinearUnits();
     658          38 :         if (IsProjected() && !osOldAuth.empty() && dfOldUnits != 0.0 &&
     659           4 :             std::abs(dfNewUnits / dfOldUnits - 1) < 1e-8)
     660             :         {
     661           1 :             SetAuthority("PROJCS", "EPSG", atoi(osOldAuth));
     662             :         }
     663             :     }
     664             : 
     665          37 :     return OGRERR_NONE;
     666             : }
     667             : 
     668             : /************************************************************************/
     669             : /*                       FindCodeFromDict()                             */
     670             : /*                                                                      */
     671             : /*      Find the code from a dict file.                                 */
     672             : /************************************************************************/
     673           0 : static int FindCodeFromDict(const char *pszDictFile, const char *CSName,
     674             :                             char *code)
     675             : {
     676             :     /* -------------------------------------------------------------------- */
     677             :     /*      Find and open file.                                             */
     678             :     /* -------------------------------------------------------------------- */
     679           0 :     const char *pszFilename = CPLFindFile("gdal", pszDictFile);
     680           0 :     if (pszFilename == nullptr)
     681           0 :         return OGRERR_UNSUPPORTED_SRS;
     682             : 
     683           0 :     VSILFILE *fp = VSIFOpenL(pszFilename, "rb");
     684           0 :     if (fp == nullptr)
     685           0 :         return OGRERR_UNSUPPORTED_SRS;
     686             : 
     687             :     /* -------------------------------------------------------------------- */
     688             :     /*      Process lines.                                                  */
     689             :     /* -------------------------------------------------------------------- */
     690           0 :     OGRErr eErr = OGRERR_UNSUPPORTED_SRS;
     691           0 :     const char *pszLine = nullptr;
     692             : 
     693           0 :     while ((pszLine = CPLReadLineL(fp)) != nullptr)
     694             :     {
     695           0 :         if (pszLine[0] == '#')
     696           0 :             continue;
     697             : 
     698           0 :         if (strstr(pszLine, CSName))
     699             :         {
     700           0 :             const char *pComma = strchr(pszLine, ',');
     701           0 :             if (pComma)
     702             :             {
     703           0 :                 strncpy(code, pszLine, pComma - pszLine);
     704           0 :                 code[pComma - pszLine] = '\0';
     705           0 :                 eErr = OGRERR_NONE;
     706             :             }
     707           0 :             break;
     708             :         }
     709             :     }
     710             : 
     711             :     /* -------------------------------------------------------------------- */
     712             :     /*      Cleanup                                                         */
     713             :     /* -------------------------------------------------------------------- */
     714           0 :     VSIFCloseL(fp);
     715             : 
     716           0 :     return eErr;
     717             : }
     718             : 
     719             : /************************************************************************/
     720             : /*                    ImportFromESRIStatePlaneWKT()                     */
     721             : /*                                                                      */
     722             : /*      Search a ESRI State Plane WKT and import it.                    */
     723             : /************************************************************************/
     724             : 
     725           3 : OGRErr OGRSpatialReference::ImportFromESRIStatePlaneWKT(int code,
     726             :                                                         const char *datumName,
     727             :                                                         const char *unitsName,
     728             :                                                         int pcsCode,
     729             :                                                         const char *csName)
     730             : {
     731             :     // If the CS name is known.
     732           3 :     if (code == 0 && !datumName && !unitsName && pcsCode == 32767 && csName)
     733             :     {
     734           0 :         char codeS[10] = {};
     735           0 :         if (FindCodeFromDict("esri_StatePlane_extra.wkt", csName, codeS) !=
     736             :             OGRERR_NONE)
     737           0 :             return OGRERR_FAILURE;
     738           0 :         return importFromDict("esri_StatePlane_extra.wkt", codeS);
     739             :     }
     740             : 
     741           3 :     int searchCode = -1;
     742           3 :     if (unitsName == nullptr)
     743           0 :         unitsName = "";
     744             : 
     745             :     // Find state plane prj str by pcs code only.
     746           3 :     if (code == 0 && !datumName && pcsCode != 32767)
     747             :     {
     748           0 :         int unitCode = 1;
     749           0 :         if (EQUAL(unitsName, "international_feet"))
     750           0 :             unitCode = 3;
     751           0 :         else if (strstr(unitsName, "feet") || strstr(unitsName, "foot"))
     752           0 :             unitCode = 2;
     753             : 
     754           0 :         for (int i = 0; statePlanePcsCodeToZoneCode[i] != 0; i += 2)
     755             :         {
     756           0 :             if (pcsCode == statePlanePcsCodeToZoneCode[i])
     757             :             {
     758           0 :                 searchCode = statePlanePcsCodeToZoneCode[i + 1];
     759           0 :                 const int unitIndex = searchCode % 10;
     760           0 :                 if ((unitCode == 1 && !(unitIndex == 0 || unitIndex == 1)) ||
     761           0 :                     (unitCode == 2 &&
     762           0 :                      !(unitIndex == 2 || unitIndex == 3 || unitIndex == 4)) ||
     763           0 :                     (unitCode == 3 && !(unitIndex == 5 || unitIndex == 6)))
     764             :                 {
     765           0 :                     searchCode -= unitIndex;
     766           0 :                     switch (unitIndex)
     767             :                     {
     768           0 :                         case 0:
     769             :                         case 3:
     770             :                         case 5:
     771           0 :                             if (unitCode == 2)
     772           0 :                                 searchCode += 3;
     773           0 :                             else if (unitCode == 3)
     774           0 :                                 searchCode += 5;
     775           0 :                             break;
     776           0 :                         case 1:
     777             :                         case 2:
     778             :                         case 6:
     779           0 :                             if (unitCode == 1)
     780           0 :                                 searchCode += 1;
     781           0 :                             if (unitCode == 2)
     782           0 :                                 searchCode += 2;
     783           0 :                             else if (unitCode == 3)
     784           0 :                                 searchCode += 6;
     785           0 :                             break;
     786           0 :                         case 4:
     787             :                             // FIXME? The following cond is not possible:
     788             :                             // if( unitCode == 2 )
     789             :                             //     searchCode += 4;
     790           0 :                             break;
     791             :                     }
     792             :                 }
     793           0 :                 break;
     794             :             }
     795           0 :         }
     796             :     }
     797             :     else  // Find state plane prj str by all inputs.
     798             :     {
     799           3 :         if (code < 0 || code > INT_MAX / 10)
     800           0 :             return OGRERR_FAILURE;
     801             : 
     802             :         // Need to have a special EPSG-ESRI zone code mapping first.
     803         360 :         for (int i = 0; statePlaneZoneMapping[i] != 0; i += 3)
     804             :         {
     805         357 :             if (code == statePlaneZoneMapping[i] &&
     806           0 :                 (statePlaneZoneMapping[i + 1] == -1 ||
     807           0 :                  pcsCode == statePlaneZoneMapping[i + 1]))
     808             :             {
     809           0 :                 code = statePlaneZoneMapping[i + 2];
     810           0 :                 break;
     811             :             }
     812             :         }
     813           3 :         searchCode = code * 10;
     814           3 :         if (!datumName)
     815             :         {
     816           0 :             CPLError(CE_Failure, CPLE_AppDefined, "datumName is NULL.");
     817           0 :             return OGRERR_FAILURE;
     818             :         }
     819           3 :         if (EQUAL(datumName, "HARN"))
     820             :         {
     821           0 :             if (EQUAL(unitsName, "international_feet"))
     822           0 :                 searchCode += 5;
     823           0 :             else if (strstr(unitsName, "feet") || strstr(unitsName, "foot"))
     824           0 :                 searchCode += 3;
     825             :         }
     826           3 :         else if (strstr(datumName, "NAD") && strstr(datumName, "83"))
     827             :         {
     828           3 :             if (EQUAL(unitsName, "meters"))
     829           0 :                 searchCode += 1;
     830           3 :             else if (EQUAL(unitsName, "international_feet"))
     831           0 :                 searchCode += 6;
     832           3 :             else if (strstr(unitsName, "feet") || strstr(unitsName, "foot"))
     833           3 :                 searchCode += 2;
     834             :         }
     835           0 :         else if (strstr(datumName, "NAD") && strstr(datumName, "27") &&
     836           0 :                  !EQUAL(unitsName, "meters"))
     837             :         {
     838           0 :             searchCode += 4;
     839             :         }
     840             :         else
     841           0 :             searchCode = -1;
     842             :     }
     843           3 :     if (searchCode > 0)
     844             :     {
     845           3 :         char codeS[20] = {};
     846           3 :         snprintf(codeS, sizeof(codeS), "%d", static_cast<int>(searchCode));
     847           3 :         return importFromDict("esri_StatePlane_extra.wkt", codeS);
     848             :     }
     849           0 :     return OGRERR_FAILURE;
     850             : }

Generated by: LCOV version 1.14