LCOV - code coverage report
Current view: top level - ogr - ogr_srs_cf1.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 353 520 67.9 %
Date: 2025-07-05 13:22:42 Functions: 11 11 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  OpenGIS Simple Features Reference Implementation
       4             :  * Purpose:  OGRSpatialReference translation to/from netCDF CF-1 georeferencing.
       5             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2004, Frank Warmerdam
       9             :  * Copyright (c) 2007-2024, Even Rouault <even.rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "ogr_core.h"
      15             : #include "ogrsf_frmts.h"
      16             : #include "ogr_srs_cf1.h"
      17             : 
      18             : /************************************************************************/
      19             : /*                         OSRImportFromCF1()                           */
      20             : /************************************************************************/
      21             : 
      22             : /**
      23             :  * \brief Import a CRS from netCDF CF-1 definitions.
      24             :  *
      25             :  * This function is the same as OGRSpatialReference::importFromCF1().
      26             :  * @since 3.9
      27             :  */
      28             : 
      29           4 : OGRErr OSRImportFromCF1(OGRSpatialReferenceH hSRS, CSLConstList papszKeyValues,
      30             :                         const char *pszUnits)
      31             : 
      32             : {
      33           4 :     VALIDATE_POINTER1(hSRS, "OSRImportFromCF1", OGRERR_FAILURE);
      34             : 
      35           4 :     return OGRSpatialReference::FromHandle(hSRS)->importFromCF1(papszKeyValues,
      36           4 :                                                                 pszUnits);
      37             : }
      38             : 
      39             : /************************************************************************/
      40             : /*                          FetchDoubleParam()                          */
      41             : /************************************************************************/
      42             : 
      43         322 : static double FetchDoubleParam(CSLConstList papszKeyValues,
      44             :                                const char *pszParam, double dfDefault)
      45             : 
      46             : {
      47         322 :     const char *pszValue = CSLFetchNameValue(papszKeyValues, pszParam);
      48         322 :     if (pszValue)
      49             :     {
      50         175 :         return CPLAtofM(pszValue);
      51             :     }
      52             : 
      53         147 :     return dfDefault;
      54             : }
      55             : 
      56             : /************************************************************************/
      57             : /*                           NCDFTokenizeArray()                        */
      58             : /************************************************************************/
      59             : 
      60             : // Parse a string, and return as a string list.
      61             : // If it an array of the form {a,b}, then tokenize it.
      62             : // Otherwise, return a copy.
      63          13 : static char **NCDFTokenizeArray(const char *pszValue)
      64             : {
      65          13 :     if (pszValue == nullptr || EQUAL(pszValue, ""))
      66           0 :         return nullptr;
      67             : 
      68          13 :     char **papszValues = nullptr;
      69          13 :     const int nLen = static_cast<int>(strlen(pszValue));
      70             : 
      71          13 :     if (pszValue[0] == '{' && nLen > 2 && pszValue[nLen - 1] == '}')
      72             :     {
      73           5 :         char *pszTemp = static_cast<char *>(CPLMalloc((nLen - 2) + 1));
      74           5 :         strncpy(pszTemp, pszValue + 1, nLen - 2);
      75           5 :         pszTemp[nLen - 2] = '\0';
      76           5 :         papszValues = CSLTokenizeString2(pszTemp, ",", CSLT_ALLOWEMPTYTOKENS);
      77           5 :         CPLFree(pszTemp);
      78             :     }
      79             :     else
      80             :     {
      81           8 :         papszValues = reinterpret_cast<char **>(CPLCalloc(2, sizeof(char *)));
      82           8 :         papszValues[0] = CPLStrdup(pszValue);
      83           8 :         papszValues[1] = nullptr;
      84             :     }
      85             : 
      86          13 :     return papszValues;
      87             : }
      88             : 
      89             : /************************************************************************/
      90             : /*                           FetchStandardParallels()                   */
      91             : /************************************************************************/
      92             : 
      93             : static std::vector<std::string>
      94          18 : FetchStandardParallels(CSLConstList papszKeyValues)
      95             : {
      96             :     // cf-1.0 tags
      97             :     const char *pszValue =
      98          18 :         CSLFetchNameValue(papszKeyValues, CF_PP_STD_PARALLEL);
      99             : 
     100          18 :     std::vector<std::string> ret;
     101          18 :     if (pszValue != nullptr)
     102             :     {
     103          30 :         CPLStringList aosValues;
     104          25 :         if (pszValue[0] != '{' &&
     105          18 :             (strchr(pszValue, ',') != nullptr ||
     106          23 :              CPLString(pszValue).Trim().find(' ') != std::string::npos))
     107             :         {
     108             :             // Some files like
     109             :             // ftp://data.knmi.nl/download/KNW-NetCDF-3D/1.0/noversion/2013/11/14/KNW-1.0_H37-ERA_NL_20131114.nc
     110             :             // do not use standard formatting for arrays, but just space
     111             :             // separated syntax
     112           2 :             aosValues = CSLTokenizeString2(pszValue, ", ", 0);
     113             :         }
     114             :         else
     115             :         {
     116          13 :             aosValues = NCDFTokenizeArray(pszValue);
     117             :         }
     118          37 :         for (int i = 0; i < aosValues.size(); i++)
     119             :         {
     120          22 :             ret.push_back(aosValues[i]);
     121             :         }
     122             :     }
     123             :     // Try gdal tags.
     124             :     else
     125             :     {
     126           3 :         pszValue = CSLFetchNameValue(papszKeyValues, CF_PP_STD_PARALLEL_1);
     127             : 
     128           3 :         if (pszValue != nullptr)
     129           0 :             ret.push_back(pszValue);
     130             : 
     131           3 :         pszValue = CSLFetchNameValue(papszKeyValues, CF_PP_STD_PARALLEL_2);
     132             : 
     133           3 :         if (pszValue != nullptr)
     134           0 :             ret.push_back(pszValue);
     135             :     }
     136             : 
     137          18 :     return ret;
     138             : }
     139             : 
     140             : /************************************************************************/
     141             : /*                          importFromCF1()                             */
     142             : /************************************************************************/
     143             : 
     144             : /**
     145             :  * \brief Import a CRS from netCDF CF-1 definitions.
     146             :  *
     147             :  * http://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings
     148             :  *
     149             :  * This function is the equivalent of the C function OSRImportFromCF1().
     150             :  *
     151             :  * @param papszKeyValues Key/value pairs from the grid mapping variable.
     152             :  * Multi-valued parameters (typically "standard_parallel") should be comma
     153             :  * separated.
     154             :  *
     155             :  * @param pszUnits Value of the "units" attribute of the X/Y arrays. May be nullptr
     156             :  *
     157             :  * @return OGRERR_NONE on success or an error code in case of failure.
     158             :  * @since 3.9
     159             :  */
     160             : 
     161          38 : OGRErr OGRSpatialReference::importFromCF1(CSLConstList papszKeyValues,
     162             :                                           const char *pszUnits)
     163             : {
     164             :     // Import from "spatial_ref" or "crs_wkt" attributes in priority
     165          38 :     const char *pszWKT = CSLFetchNameValue(papszKeyValues, NCDF_SPATIAL_REF);
     166          38 :     if (!pszWKT)
     167             :     {
     168          37 :         pszWKT = CSLFetchNameValue(papszKeyValues, NCDF_CRS_WKT);
     169             :     }
     170          38 :     if (pszWKT)
     171           2 :         return importFromWkt(pszWKT);
     172             : 
     173             :     const char *pszGridMappingName =
     174          36 :         CSLFetchNameValue(papszKeyValues, CF_GRD_MAPPING_NAME);
     175             : 
     176             :     // Some files such as
     177             :     // http://www.ecad.eu/download/ensembles/data/Grid_0.44deg_rot/tg_0.44deg_rot_v16.0.nc.gz
     178             :     // lack an explicit projection_var:grid_mapping_name attribute
     179          40 :     if (pszGridMappingName == nullptr &&
     180           4 :         CSLFetchNameValue(papszKeyValues, CF_PP_GRID_NORTH_POLE_LONGITUDE) !=
     181             :             nullptr)
     182             :     {
     183           3 :         pszGridMappingName = CF_PT_ROTATED_LATITUDE_LONGITUDE;
     184             :     }
     185             : 
     186          36 :     if (pszGridMappingName == nullptr)
     187           1 :         return OGRERR_FAILURE;
     188             : 
     189          35 :     bool bGotGeogCS = false;
     190             : 
     191             :     // Check for datum/spheroid information.
     192             :     double dfEarthRadius =
     193          35 :         FetchDoubleParam(papszKeyValues, CF_PP_EARTH_RADIUS, -1.0);
     194             : 
     195             :     const double dfLonPrimeMeridian =
     196          35 :         FetchDoubleParam(papszKeyValues, CF_PP_LONG_PRIME_MERIDIAN, 0.0);
     197             : 
     198             :     const char *pszPMName =
     199          35 :         CSLFetchNameValue(papszKeyValues, CF_PRIME_MERIDIAN_NAME);
     200             : 
     201             :     // Should try to find PM name from its value if not Greenwich.
     202          35 :     if (pszPMName == nullptr && !CPLIsEqual(dfLonPrimeMeridian, 0.0))
     203           0 :         pszPMName = "unknown";
     204             : 
     205             :     double dfInverseFlattening =
     206          35 :         FetchDoubleParam(papszKeyValues, CF_PP_INVERSE_FLATTENING, -1.0);
     207             : 
     208             :     double dfSemiMajorAxis =
     209          35 :         FetchDoubleParam(papszKeyValues, CF_PP_SEMI_MAJOR_AXIS, -1.0);
     210             : 
     211             :     const double dfSemiMinorAxis =
     212          35 :         FetchDoubleParam(papszKeyValues, CF_PP_SEMI_MINOR_AXIS, -1.0);
     213             : 
     214             :     // See if semi-major exists if radius doesn't.
     215          35 :     if (dfEarthRadius < 0.0)
     216          34 :         dfEarthRadius = dfSemiMajorAxis;
     217             : 
     218             :     // If still no radius, check old tag.
     219          35 :     if (dfEarthRadius < 0.0)
     220             :         dfEarthRadius =
     221          14 :             FetchDoubleParam(papszKeyValues, CF_PP_EARTH_RADIUS_OLD, -1.0);
     222             : 
     223             :     const char *pszEllipsoidName =
     224          35 :         CSLFetchNameValue(papszKeyValues, CF_REFERENCE_ELLIPSOID_NAME);
     225             : 
     226             :     const char *pszDatumName =
     227          35 :         CSLFetchNameValue(papszKeyValues, CF_HORIZONTAL_DATUM_NAME);
     228             : 
     229             :     const char *pszGeogName =
     230          35 :         CSLFetchNameValue(papszKeyValues, CF_GEOGRAPHIC_CRS_NAME);
     231          35 :     if (pszGeogName == nullptr)
     232          34 :         pszGeogName = "unknown";
     233             : 
     234          35 :     bool bRotatedPole = false;
     235             : 
     236             :     // Has radius value.
     237          35 :     if (dfEarthRadius > 0.0)
     238             :     {
     239             :         // Check for inv_flat tag.
     240          25 :         if (dfInverseFlattening < 0.0)
     241             :         {
     242             :             // No inv_flat tag, check for semi_minor.
     243           7 :             if (dfSemiMinorAxis < 0.0)
     244             :             {
     245             :                 // No way to get inv_flat, use sphere.
     246           4 :                 SetGeogCS(pszGeogName, pszDatumName,
     247             :                           pszEllipsoidName ? pszEllipsoidName : "Sphere",
     248             :                           dfEarthRadius, 0.0, pszPMName, dfLonPrimeMeridian);
     249           4 :                 bGotGeogCS = true;
     250             :             }
     251             :             else
     252             :             {
     253           3 :                 if (dfSemiMajorAxis < 0.0)
     254           0 :                     dfSemiMajorAxis = dfEarthRadius;
     255             :                 // set inv_flat using semi_minor/major
     256             :                 dfInverseFlattening =
     257           3 :                     OSRCalcInvFlattening(dfSemiMajorAxis, dfSemiMinorAxis);
     258             : 
     259           3 :                 SetGeogCS(pszGeogName, pszDatumName,
     260             :                           pszEllipsoidName ? pszEllipsoidName : "Spheroid",
     261             :                           dfEarthRadius, dfInverseFlattening, pszPMName,
     262             :                           dfLonPrimeMeridian);
     263           3 :                 bGotGeogCS = true;
     264             :             }
     265             :         }
     266             :         else
     267             :         {
     268          18 :             SetGeogCS(pszGeogName, pszDatumName,
     269             :                       pszEllipsoidName ? pszEllipsoidName : "Spheroid",
     270             :                       dfEarthRadius, dfInverseFlattening, pszPMName,
     271             :                       dfLonPrimeMeridian);
     272          18 :             bGotGeogCS = true;
     273             :         }
     274             : 
     275          25 :         if (bGotGeogCS)
     276          25 :             CPLDebug("GDAL_netCDF", "got spheroid from CF: (%f , %f)",
     277             :                      dfEarthRadius, dfInverseFlattening);
     278             :     }
     279             : 
     280             :     // Transverse Mercator.
     281          35 :     if (EQUAL(pszGridMappingName, CF_PT_TM))
     282             :     {
     283             :         const double dfScale =
     284           7 :             FetchDoubleParam(papszKeyValues, CF_PP_SCALE_FACTOR_MERIDIAN, 1.0);
     285             : 
     286             :         const double dfCenterLon =
     287           7 :             FetchDoubleParam(papszKeyValues, CF_PP_LONG_CENTRAL_MERIDIAN, 0.0);
     288             : 
     289             :         const double dfCenterLat =
     290           7 :             FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     291             : 
     292             :         const double dfFalseEasting =
     293           7 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     294             : 
     295             :         const double dfFalseNorthing =
     296           7 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     297             : 
     298           7 :         SetTM(dfCenterLat, dfCenterLon, dfScale, dfFalseEasting,
     299             :               dfFalseNorthing);
     300             : 
     301           7 :         if (!bGotGeogCS)
     302           0 :             SetWellKnownGeogCS("WGS84");
     303             :     }
     304             : 
     305             :     // Albers Equal Area.
     306          35 :     if (EQUAL(pszGridMappingName, CF_PT_AEA))
     307             :     {
     308             :         const double dfCenterLon =
     309           2 :             FetchDoubleParam(papszKeyValues, CF_PP_LONG_CENTRAL_MERIDIAN, 0.0);
     310             : 
     311             :         const double dfFalseEasting =
     312           2 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     313             : 
     314             :         const double dfFalseNorthing =
     315           2 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     316             : 
     317           4 :         const auto aosStdParallels = FetchStandardParallels(papszKeyValues);
     318             : 
     319           2 :         double dfStdP1 = 0;
     320           2 :         double dfStdP2 = 0;
     321           2 :         if (aosStdParallels.size() == 1)
     322             :         {
     323             :             // TODO CF-1 standard says it allows AEA to be encoded
     324             :             // with only 1 standard parallel.  How should this
     325             :             // actually map to a 2StdP OGC WKT version?
     326           0 :             CPLError(CE_Warning, CPLE_NotSupported,
     327             :                      "NetCDF driver import of AEA-1SP is not tested, "
     328             :                      "using identical std. parallels.");
     329           0 :             dfStdP1 = CPLAtofM(aosStdParallels[0].c_str());
     330           0 :             dfStdP2 = dfStdP1;
     331             :         }
     332           2 :         else if (aosStdParallels.size() == 2)
     333             :         {
     334           2 :             dfStdP1 = CPLAtofM(aosStdParallels[0].c_str());
     335           2 :             dfStdP2 = CPLAtofM(aosStdParallels[1].c_str());
     336             :         }
     337             :         // Old default.
     338             :         else
     339             :         {
     340             :             dfStdP1 =
     341           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_STD_PARALLEL_1, 0.0);
     342             : 
     343             :             dfStdP2 =
     344           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_STD_PARALLEL_2, 0.0);
     345             :         }
     346             : 
     347             :         const double dfCenterLat =
     348           2 :             FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     349             : 
     350           2 :         SetACEA(dfStdP1, dfStdP2, dfCenterLat, dfCenterLon, dfFalseEasting,
     351             :                 dfFalseNorthing);
     352             : 
     353           2 :         if (!bGotGeogCS)
     354           0 :             SetWellKnownGeogCS("WGS84");
     355             :     }
     356             : 
     357             :     // Cylindrical Equal Area
     358          33 :     else if (EQUAL(pszGridMappingName, CF_PT_CEA) ||
     359          33 :              EQUAL(pszGridMappingName, CF_PT_LCEA))
     360             :     {
     361           0 :         const auto aosStdParallels = FetchStandardParallels(papszKeyValues);
     362             : 
     363           0 :         double dfStdP1 = 0;
     364           0 :         if (!aosStdParallels.empty())
     365             :         {
     366           0 :             dfStdP1 = CPLAtofM(aosStdParallels[0].c_str());
     367             :         }
     368             :         else
     369             :         {
     370             :             // TODO: Add support for 'scale_factor_at_projection_origin'
     371             :             // variant to standard parallel.  Probably then need to calc
     372             :             // a std parallel equivalent.
     373           0 :             CPLError(CE_Failure, CPLE_NotSupported,
     374             :                      "NetCDF driver does not support import of CF-1 LCEA "
     375             :                      "'scale_factor_at_projection_origin' variant yet.");
     376             :         }
     377             : 
     378             :         const double dfCentralMeridian =
     379           0 :             FetchDoubleParam(papszKeyValues, CF_PP_LONG_CENTRAL_MERIDIAN, 0.0);
     380             : 
     381             :         const double dfFalseEasting =
     382           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     383             : 
     384             :         const double dfFalseNorthing =
     385           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     386             : 
     387           0 :         SetCEA(dfStdP1, dfCentralMeridian, dfFalseEasting, dfFalseNorthing);
     388             : 
     389           0 :         if (!bGotGeogCS)
     390           0 :             SetWellKnownGeogCS("WGS84");
     391             :     }
     392             : 
     393             :     // lambert_azimuthal_equal_area.
     394          33 :     else if (EQUAL(pszGridMappingName, CF_PT_LAEA))
     395             :     {
     396             :         const double dfCenterLon =
     397           1 :             FetchDoubleParam(papszKeyValues, CF_PP_LON_PROJ_ORIGIN, 0.0);
     398             : 
     399             :         const double dfCenterLat =
     400           1 :             FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     401             : 
     402             :         const double dfFalseEasting =
     403           1 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     404             : 
     405             :         const double dfFalseNorthing =
     406           1 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     407             : 
     408           1 :         SetLAEA(dfCenterLat, dfCenterLon, dfFalseEasting, dfFalseNorthing);
     409             : 
     410           1 :         if (!bGotGeogCS)
     411           0 :             SetWellKnownGeogCS("WGS84");
     412             : 
     413           2 :         if (GetAttrValue("DATUM") != nullptr &&
     414           1 :             EQUAL(GetAttrValue("DATUM"), "WGS_1984"))
     415             :         {
     416           0 :             SetProjCS("LAEA (WGS84)");
     417             :         }
     418             :     }
     419             : 
     420             :     // Azimuthal Equidistant.
     421          32 :     else if (EQUAL(pszGridMappingName, CF_PT_AE))
     422             :     {
     423             :         const double dfCenterLon =
     424           0 :             FetchDoubleParam(papszKeyValues, CF_PP_LON_PROJ_ORIGIN, 0.0);
     425             : 
     426             :         const double dfCenterLat =
     427           0 :             FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     428             : 
     429             :         const double dfFalseEasting =
     430           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     431             : 
     432             :         const double dfFalseNorthing =
     433           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     434             : 
     435           0 :         SetAE(dfCenterLat, dfCenterLon, dfFalseEasting, dfFalseNorthing);
     436             : 
     437           0 :         if (!bGotGeogCS)
     438           0 :             SetWellKnownGeogCS("WGS84");
     439             :     }
     440             : 
     441             :     // Lambert conformal conic.
     442          32 :     else if (EQUAL(pszGridMappingName, CF_PT_LCC))
     443             :     {
     444             :         const double dfCenterLon =
     445           8 :             FetchDoubleParam(papszKeyValues, CF_PP_LONG_CENTRAL_MERIDIAN, 0.0);
     446             : 
     447             :         const double dfCenterLat =
     448           8 :             FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     449             : 
     450             :         const double dfFalseEasting =
     451           8 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     452             : 
     453             :         const double dfFalseNorthing =
     454           8 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     455             : 
     456          16 :         const auto aosStdParallels = FetchStandardParallels(papszKeyValues);
     457             : 
     458             :         // 2SP variant.
     459           8 :         if (aosStdParallels.size() == 2)
     460             :         {
     461           5 :             const double dfStdP1 = CPLAtofM(aosStdParallels[0].c_str());
     462           5 :             const double dfStdP2 = CPLAtofM(aosStdParallels[1].c_str());
     463           5 :             SetLCC(dfStdP1, dfStdP2, dfCenterLat, dfCenterLon, dfFalseEasting,
     464             :                    dfFalseNorthing);
     465             :         }
     466             :         // 1SP variant (with standard_parallel or center lon).
     467             :         // See comments in netcdfdataset.h for this projection.
     468             :         else
     469             :         {
     470           3 :             double dfScale = FetchDoubleParam(papszKeyValues,
     471             :                                               CF_PP_SCALE_FACTOR_ORIGIN, -1.0);
     472             : 
     473             :             // CF definition, without scale factor.
     474           3 :             if (CPLIsEqual(dfScale, -1.0))
     475             :             {
     476             :                 double dfStdP1;
     477             :                 // With standard_parallel.
     478           3 :                 if (aosStdParallels.size() == 1)
     479           3 :                     dfStdP1 = CPLAtofM(aosStdParallels[0].c_str());
     480             :                 // With center lon instead.
     481             :                 else
     482           0 :                     dfStdP1 = dfCenterLat;
     483             :                 // dfStdP2 = dfStdP1;
     484             : 
     485             :                 // Test if we should actually compute scale factor.
     486           3 :                 if (!CPLIsEqual(dfStdP1, dfCenterLat))
     487             :                 {
     488           0 :                     CPLError(CE_Warning, CPLE_NotSupported,
     489             :                              "NetCDF driver import of LCC-1SP with "
     490             :                              "standard_parallel1 != "
     491             :                              "latitude_of_projection_origin "
     492             :                              "(which forces a computation of scale_factor) "
     493             :                              "is experimental (bug #3324)");
     494             :                     // Use Snyder eq. 15-4 to compute dfScale from
     495             :                     // dfStdP1 and dfCenterLat.  Only tested for
     496             :                     // dfStdP1=dfCenterLat and (25,26), needs more data
     497             :                     // for testing.  Other option: use the 2SP variant -
     498             :                     // how to compute new standard parallels?
     499           0 :                     dfScale =
     500           0 :                         (cos(dfStdP1) *
     501           0 :                          pow(tan(M_PI / 4 + dfStdP1 / 2), sin(dfStdP1))) /
     502           0 :                         (cos(dfCenterLat) * pow(tan(M_PI / 4 + dfCenterLat / 2),
     503             :                                                 sin(dfCenterLat)));
     504             :                 }
     505             :                 // Default is 1.0.
     506             :                 else
     507             :                 {
     508           3 :                     dfScale = 1.0;
     509             :                 }
     510             : 
     511           3 :                 SetLCC1SP(dfCenterLat, dfCenterLon, dfScale, dfFalseEasting,
     512             :                           dfFalseNorthing);
     513             :                 // Store dfStdP1 so we can output it to CF later.
     514           3 :                 SetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, dfStdP1);
     515             :             }
     516             :             // OGC/PROJ.4 definition with scale factor.
     517             :             else
     518             :             {
     519           0 :                 SetLCC1SP(dfCenterLat, dfCenterLon, dfScale, dfFalseEasting,
     520             :                           dfFalseNorthing);
     521             :             }
     522             :         }
     523             : 
     524           8 :         if (!bGotGeogCS)
     525           3 :             SetWellKnownGeogCS("WGS84");
     526             :     }
     527             : 
     528             :     // Is this Latitude/Longitude Grid explicitly?
     529             : 
     530          24 :     else if (EQUAL(pszGridMappingName, CF_PT_LATITUDE_LONGITUDE))
     531             :     {
     532           3 :         if (!bGotGeogCS)
     533           0 :             SetWellKnownGeogCS("WGS84");
     534             :     }
     535             : 
     536             :     // Mercator.
     537          21 :     else if (EQUAL(pszGridMappingName, CF_PT_MERCATOR))
     538             :     {
     539             : 
     540             :         // If there is a standard_parallel, know it is Mercator 2SP.
     541           0 :         const auto aosStdParallels = FetchStandardParallels(papszKeyValues);
     542             : 
     543           0 :         if (!aosStdParallels.empty())
     544             :         {
     545             :             // CF-1 Mercator 2SP always has lat centered at equator.
     546           0 :             const double dfStdP1 = CPLAtofM(aosStdParallels[0].c_str());
     547             : 
     548           0 :             const double dfCenterLat = 0.0;
     549             : 
     550             :             const double dfCenterLon =
     551           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_LON_PROJ_ORIGIN, 0.0);
     552             : 
     553             :             const double dfFalseEasting =
     554           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     555             : 
     556             :             const double dfFalseNorthing =
     557           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     558             : 
     559           0 :             SetMercator2SP(dfStdP1, dfCenterLat, dfCenterLon, dfFalseEasting,
     560             :                            dfFalseNorthing);
     561             :         }
     562             :         else
     563             :         {
     564             :             const double dfCenterLon =
     565           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_LON_PROJ_ORIGIN, 0.0);
     566             : 
     567             :             const double dfCenterLat =
     568           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     569             : 
     570           0 :             const double dfScale = FetchDoubleParam(
     571             :                 papszKeyValues, CF_PP_SCALE_FACTOR_ORIGIN, 1.0);
     572             : 
     573             :             const double dfFalseEasting =
     574           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     575             : 
     576             :             const double dfFalseNorthing =
     577           0 :                 FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     578             : 
     579           0 :             SetMercator(dfCenterLat, dfCenterLon, dfScale, dfFalseEasting,
     580             :                         dfFalseNorthing);
     581             :         }
     582             : 
     583           0 :         if (!bGotGeogCS)
     584           0 :             SetWellKnownGeogCS("WGS84");
     585             :     }
     586             : 
     587             :     // Orthographic.
     588          21 :     else if (EQUAL(pszGridMappingName, CF_PT_ORTHOGRAPHIC))
     589             :     {
     590             :         const double dfCenterLon =
     591           0 :             FetchDoubleParam(papszKeyValues, CF_PP_LON_PROJ_ORIGIN, 0.0);
     592             : 
     593             :         const double dfCenterLat =
     594           0 :             FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     595             : 
     596             :         const double dfFalseEasting =
     597           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     598             : 
     599             :         const double dfFalseNorthing =
     600           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     601             : 
     602           0 :         SetOrthographic(dfCenterLat, dfCenterLon, dfFalseEasting,
     603             :                         dfFalseNorthing);
     604             : 
     605           0 :         if (!bGotGeogCS)
     606           0 :             SetWellKnownGeogCS("WGS84");
     607             :     }
     608             : 
     609             :     // Polar Stereographic.
     610          21 :     else if (EQUAL(pszGridMappingName, CF_PT_POLAR_STEREO))
     611             :     {
     612          16 :         const auto aosStdParallels = FetchStandardParallels(papszKeyValues);
     613             : 
     614             :         const double dfCenterLon =
     615           8 :             FetchDoubleParam(papszKeyValues, CF_PP_VERT_LONG_FROM_POLE, 0.0);
     616             : 
     617             :         const double dfFalseEasting =
     618           8 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     619             : 
     620             :         const double dfFalseNorthing =
     621           8 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     622             : 
     623             :         // CF allows the use of standard_parallel (lat_ts) OR
     624             :         // scale_factor (k0), make sure we have standard_parallel, using
     625             :         // Snyder eq. 22-7 with k=1 and lat=standard_parallel.
     626           8 :         if (!aosStdParallels.empty())
     627             :         {
     628           5 :             const double dfStdP1 = CPLAtofM(aosStdParallels[0].c_str());
     629             : 
     630             :             // Polar Stereographic Variant B with latitude of standard
     631             :             // parallel
     632           5 :             SetPS(dfStdP1, dfCenterLon, 1.0, dfFalseEasting, dfFalseNorthing);
     633             :         }
     634             :         else
     635             :         {
     636             :             // Fetch latitude_of_projection_origin (+90/-90).
     637             :             double dfLatProjOrigin =
     638           3 :                 FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     639           3 :             if (!CPLIsEqual(dfLatProjOrigin, 90.0) &&
     640           0 :                 !CPLIsEqual(dfLatProjOrigin, -90.0))
     641             :             {
     642           0 :                 CPLError(CE_Failure, CPLE_NotSupported,
     643             :                          "Polar Stereographic must have a %s "
     644             :                          "parameter equal to +90 or -90.",
     645             :                          CF_PP_LAT_PROJ_ORIGIN);
     646           0 :                 dfLatProjOrigin = 90.0;
     647             :             }
     648             : 
     649           3 :             const double dfScale = FetchDoubleParam(
     650             :                 papszKeyValues, CF_PP_SCALE_FACTOR_ORIGIN, 1.0);
     651             : 
     652             :             // Polar Stereographic Variant A with scale factor at
     653             :             // natural origin and latitude of origin = +/- 90
     654           3 :             SetPS(dfLatProjOrigin, dfCenterLon, dfScale, dfFalseEasting,
     655             :                   dfFalseNorthing);
     656             :         }
     657             : 
     658           8 :         if (!bGotGeogCS)
     659           6 :             SetWellKnownGeogCS("WGS84");
     660             :     }
     661             : 
     662             :     // Stereographic.
     663          13 :     else if (EQUAL(pszGridMappingName, CF_PT_STEREO))
     664             :     {
     665             :         const double dfCenterLon =
     666           0 :             FetchDoubleParam(papszKeyValues, CF_PP_LON_PROJ_ORIGIN, 0.0);
     667             : 
     668             :         const double dfCenterLat =
     669           0 :             FetchDoubleParam(papszKeyValues, CF_PP_LAT_PROJ_ORIGIN, 0.0);
     670             : 
     671             :         const double dfScale =
     672           0 :             FetchDoubleParam(papszKeyValues, CF_PP_SCALE_FACTOR_ORIGIN, 1.0);
     673             : 
     674             :         const double dfFalseEasting =
     675           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     676             : 
     677             :         const double dfFalseNorthing =
     678           0 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     679             : 
     680           0 :         SetStereographic(dfCenterLat, dfCenterLon, dfScale, dfFalseEasting,
     681             :                          dfFalseNorthing);
     682             : 
     683           0 :         if (!bGotGeogCS)
     684           0 :             SetWellKnownGeogCS("WGS84");
     685             :     }
     686             : 
     687             :     // Geostationary.
     688          13 :     else if (EQUAL(pszGridMappingName, CF_PT_GEOS))
     689             :     {
     690             :         const double dfCenterLon =
     691           3 :             FetchDoubleParam(papszKeyValues, CF_PP_LON_PROJ_ORIGIN, 0.0);
     692             : 
     693           3 :         const double dfSatelliteHeight = FetchDoubleParam(
     694             :             papszKeyValues, CF_PP_PERSPECTIVE_POINT_HEIGHT, 35785831.0);
     695             : 
     696             :         const char *pszSweepAxisAngle =
     697           3 :             CSLFetchNameValue(papszKeyValues, CF_PP_SWEEP_ANGLE_AXIS);
     698             : 
     699             :         const double dfFalseEasting =
     700           3 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_EASTING, 0.0);
     701             : 
     702             :         const double dfFalseNorthing =
     703           3 :             FetchDoubleParam(papszKeyValues, CF_PP_FALSE_NORTHING, 0.0);
     704             : 
     705           3 :         SetGEOS(dfCenterLon, dfSatelliteHeight, dfFalseEasting,
     706             :                 dfFalseNorthing);
     707             : 
     708           3 :         if (!bGotGeogCS)
     709           1 :             SetWellKnownGeogCS("WGS84");
     710             : 
     711           3 :         if (pszSweepAxisAngle != nullptr && EQUAL(pszSweepAxisAngle, "x"))
     712             :         {
     713           2 :             char *pszProj4 = nullptr;
     714           2 :             exportToProj4(&pszProj4);
     715           4 :             CPLString osProj4 = pszProj4;
     716           2 :             osProj4 += " +sweep=x";
     717           2 :             SetExtension(GetRoot()->GetValue(), "PROJ4", osProj4);
     718           2 :             CPLFree(pszProj4);
     719             :         }
     720             :     }
     721             : 
     722          10 :     else if (EQUAL(pszGridMappingName, CF_PT_ROTATED_LATITUDE_LONGITUDE))
     723             :     {
     724           3 :         const double dfGridNorthPoleLong = FetchDoubleParam(
     725             :             papszKeyValues, CF_PP_GRID_NORTH_POLE_LONGITUDE, 0.0);
     726           3 :         const double dfGridNorthPoleLat = FetchDoubleParam(
     727             :             papszKeyValues, CF_PP_GRID_NORTH_POLE_LATITUDE, 0.0);
     728           3 :         const double dfNorthPoleGridLong = FetchDoubleParam(
     729             :             papszKeyValues, CF_PP_NORTH_POLE_GRID_LONGITUDE, 0.0);
     730             : 
     731           3 :         if (!bGotGeogCS)
     732           0 :             SetWellKnownGeogCS("WGS84");
     733             : 
     734           3 :         bRotatedPole = true;
     735           3 :         SetDerivedGeogCRSWithPoleRotationNetCDFCFConvention(
     736             :             "Rotated_pole", dfGridNorthPoleLat, dfGridNorthPoleLong,
     737             :             dfNorthPoleGridLong);
     738             :     }
     739             : 
     740          35 :     if (IsProjected() && !bRotatedPole)
     741             :     {
     742             :         const char *pszProjectedCRSName =
     743          29 :             CSLFetchNameValue(papszKeyValues, CF_PROJECTED_CRS_NAME);
     744          29 :         if (pszProjectedCRSName)
     745           0 :             SetProjCS(pszProjectedCRSName);
     746             :     }
     747             : 
     748             :     // Add units to PROJCS.
     749          35 :     if (IsGeographic() && !bRotatedPole)
     750             :     {
     751           3 :         SetAngularUnits(SRS_UA_DEGREE, CPLAtof(SRS_UA_DEGREE_CONV));
     752           3 :         SetAuthority("GEOGCS|UNIT", "EPSG", 9122);
     753             :     }
     754          32 :     else if (pszUnits != nullptr && !EQUAL(pszUnits, "") && !bRotatedPole)
     755             :     {
     756          21 :         if (EQUAL(pszUnits, "m") || EQUAL(pszUnits, "metre") ||
     757          10 :             EQUAL(pszUnits, "meter"))
     758             :         {
     759          11 :             SetLinearUnits("metre", 1.0);
     760          11 :             SetAuthority("PROJCS|UNIT", "EPSG", 9001);
     761             :         }
     762          10 :         else if (EQUAL(pszUnits, "km"))
     763             :         {
     764           6 :             SetLinearUnits("kilometre", 1000.0);
     765           6 :             SetAuthority("PROJCS|UNIT", "EPSG", 9036);
     766             :         }
     767           4 :         else if (EQUAL(pszUnits, "US_survey_foot") ||
     768           4 :                  EQUAL(pszUnits, "US_survey_feet"))
     769             :         {
     770           0 :             SetLinearUnits("US survey foot", CPLAtof(SRS_UL_US_FOOT_CONV));
     771           0 :             SetAuthority("PROJCS|UNIT", "EPSG", 9003);
     772             :         }
     773             :         else
     774             :         {
     775           4 :             CPLError(CE_Warning, CPLE_AppDefined,
     776             :                      "Unhandled X/Y axis unit %s. SRS will ignore "
     777             :                      "axis unit and be likely wrong.",
     778             :                      pszUnits);
     779             :         }
     780             :     }
     781             : 
     782          35 :     return OGRERR_NONE;
     783             : }
     784             : 
     785             : /* -------------------------------------------------------------------- */
     786             : /*         CF-1 to GDAL mappings                                        */
     787             : /* -------------------------------------------------------------------- */
     788             : 
     789             : /* Following are a series of mappings from CF-1 convention parameters
     790             :  * for each projection, to the equivalent in OGC WKT used internally by GDAL.
     791             :  * See: http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.5/apf.html
     792             :  */
     793             : 
     794             : /* A struct allowing us to map between GDAL(OGC WKT) and CF-1 attributes */
     795             : typedef struct
     796             : {
     797             :     const char *CF_ATT;
     798             :     const char *WKT_ATT;
     799             :     // TODO: mappings may need default values, like scale factor?
     800             :     // double defval;
     801             : } oNetcdfSRS_PP;
     802             : 
     803             : // default mappings, for the generic case
     804             : /* These 'generic' mappings are based on what was previously in the
     805             :    poNetCDFSRS struct. They will be used as a fallback in case none
     806             :    of the others match (i.e. you are exporting a projection that has
     807             :    no CF-1 equivalent).
     808             :    They are not used for known CF-1 projections since there is not a
     809             :    unique 2-way projection-independent
     810             :    mapping between OGC WKT params and CF-1 ones: it varies per-projection.
     811             : */
     812             : 
     813             : static const oNetcdfSRS_PP poGenericMappings[] = {
     814             :     /* scale_factor is handled as a special case, write 2 values */
     815             :     {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
     816             :     {CF_PP_STD_PARALLEL_2, SRS_PP_STANDARD_PARALLEL_2},
     817             :     {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
     818             :     {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_LONGITUDE_OF_CENTER},
     819             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_LONGITUDE_OF_ORIGIN},
     820             :     // Multiple mappings to LAT_PROJ_ORIGIN
     821             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
     822             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
     823             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     824             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     825             :     {nullptr, nullptr},
     826             : };
     827             : 
     828             : // Albers equal area
     829             : //
     830             : // grid_mapping_name = albers_conical_equal_area
     831             : // WKT: Albers_Conic_Equal_Area
     832             : // EPSG:9822
     833             : //
     834             : // Map parameters:
     835             : //
     836             : //    * standard_parallel - There may be 1 or 2 values.
     837             : //    * longitude_of_central_meridian
     838             : //    * latitude_of_projection_origin
     839             : //    * false_easting
     840             : //    * false_northing
     841             : //
     842             : static const oNetcdfSRS_PP poAEAMappings[] = {
     843             :     {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
     844             :     {CF_PP_STD_PARALLEL_2, SRS_PP_STANDARD_PARALLEL_2},
     845             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
     846             :     {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_LONGITUDE_OF_CENTER},
     847             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     848             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     849             :     {nullptr, nullptr}};
     850             : 
     851             : // Azimuthal equidistant
     852             : //
     853             : // grid_mapping_name = azimuthal_equidistant
     854             : // WKT: Azimuthal_Equidistant
     855             : //
     856             : // Map parameters:
     857             : //
     858             : //    * longitude_of_projection_origin
     859             : //    * latitude_of_projection_origin
     860             : //    * false_easting
     861             : //    * false_northing
     862             : //
     863             : static const oNetcdfSRS_PP poAEMappings[] = {
     864             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
     865             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_LONGITUDE_OF_CENTER},
     866             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     867             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     868             :     {nullptr, nullptr}};
     869             : 
     870             : // Lambert azimuthal equal area
     871             : //
     872             : // grid_mapping_name = lambert_azimuthal_equal_area
     873             : // WKT: Lambert_Azimuthal_Equal_Area
     874             : //
     875             : // Map parameters:
     876             : //
     877             : //    * longitude_of_projection_origin
     878             : //    * latitude_of_projection_origin
     879             : //    * false_easting
     880             : //    * false_northing
     881             : //
     882             : static const oNetcdfSRS_PP poLAEAMappings[] = {
     883             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
     884             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_LONGITUDE_OF_CENTER},
     885             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     886             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     887             :     {nullptr, nullptr}};
     888             : 
     889             : // Lambert conformal
     890             : //
     891             : // grid_mapping_name = lambert_conformal_conic
     892             : // WKT: Lambert_Conformal_Conic_1SP / Lambert_Conformal_Conic_2SP
     893             : //
     894             : // Map parameters:
     895             : //
     896             : //    * standard_parallel - There may be 1 or 2 values.
     897             : //    * longitude_of_central_meridian
     898             : //    * latitude_of_projection_origin
     899             : //    * false_easting
     900             : //    * false_northing
     901             : //
     902             : // See
     903             : // http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html
     904             : 
     905             : // Lambert conformal conic - 1SP
     906             : /* See bug # 3324
     907             :    It seems that the missing scale factor can be computed from
     908             :    standard_parallel1 and latitude_of_projection_origin. If both are equal (the
     909             :    common case) then scale factor=1, else use Snyder eq. 15-4. We save in the
     910             :    WKT standard_parallel1 for export to CF, but do not export scale factor. If a
     911             :    WKT has a scale factor != 1 and no standard_parallel1 then export is not CF,
     912             :    but we output scale factor for compat. is there a formula for that?
     913             : */
     914             : static const oNetcdfSRS_PP poLCC1SPMappings[] = {
     915             :     {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
     916             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
     917             :     {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
     918             :     {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR}, /* special case */
     919             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     920             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     921             :     {nullptr, nullptr}};
     922             : 
     923             : // Lambert conformal conic - 2SP
     924             : static const oNetcdfSRS_PP poLCC2SPMappings[] = {
     925             :     {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
     926             :     {CF_PP_STD_PARALLEL_2, SRS_PP_STANDARD_PARALLEL_2},
     927             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
     928             :     {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
     929             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     930             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     931             :     {nullptr, nullptr}};
     932             : 
     933             : // Lambert cylindrical equal area
     934             : //
     935             : // grid_mapping_name = lambert_cylindrical_equal_area
     936             : // WKT: Cylindrical_Equal_Area
     937             : // EPSG:9834 (Spherical) and EPSG:9835
     938             : //
     939             : // Map parameters:
     940             : //
     941             : //    * longitude_of_central_meridian
     942             : //    * either standard_parallel or scale_factor_at_projection_origin
     943             : //    * false_easting
     944             : //    * false_northing
     945             : //
     946             : // NB: CF-1 specifies a 'scale_factor_at_projection' alternative
     947             : //  to std_parallel ... but no reference to this in EPSG/remotesensing.org
     948             : //  ignore for now.
     949             : //
     950             : static const oNetcdfSRS_PP poLCEAMappings[] = {
     951             :     {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
     952             :     {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
     953             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     954             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     955             :     {nullptr, nullptr}};
     956             : 
     957             : // Latitude-Longitude
     958             : //
     959             : // grid_mapping_name = latitude_longitude
     960             : //
     961             : // Map parameters:
     962             : //
     963             : //    * None
     964             : //
     965             : // NB: handled as a special case - !isProjected()
     966             : 
     967             : // Mercator
     968             : //
     969             : // grid_mapping_name = mercator
     970             : // WKT: Mercator_1SP / Mercator_2SP
     971             : //
     972             : // Map parameters:
     973             : //
     974             : //    * longitude_of_projection_origin
     975             : //    * either standard_parallel or scale_factor_at_projection_origin
     976             : //    * false_easting
     977             : //    * false_northing
     978             : 
     979             : // Mercator 1 Standard Parallel (EPSG:9804)
     980             : static const oNetcdfSRS_PP poM1SPMappings[] = {
     981             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
     982             :     // LAT_PROJ_ORIGIN is always equator (0) in CF-1
     983             :     {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR},
     984             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     985             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     986             :     {nullptr, nullptr}};
     987             : 
     988             : // Mercator 2 Standard Parallel
     989             : static const oNetcdfSRS_PP poM2SPMappings[] = {
     990             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
     991             :     {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
     992             :     // From best understanding of this projection, only
     993             :     // actually specify one SP - it is the same N/S of equator.
     994             :     // {CF_PP_STD_PARALLEL_2, SRS_PP_LATITUDE_OF_ORIGIN},
     995             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
     996             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
     997             :     {nullptr, nullptr}};
     998             : 
     999             : // Orthographic
    1000             : // grid_mapping_name = orthographic
    1001             : // WKT: Orthographic
    1002             : //
    1003             : // Map parameters:
    1004             : //
    1005             : //    * longitude_of_projection_origin
    1006             : //    * latitude_of_projection_origin
    1007             : //    * false_easting
    1008             : //    * false_northing
    1009             : //
    1010             : static const oNetcdfSRS_PP poOrthoMappings[] = {
    1011             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
    1012             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
    1013             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
    1014             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
    1015             :     {nullptr, nullptr}};
    1016             : 
    1017             : // Polar stereographic
    1018             : //
    1019             : // grid_mapping_name = polar_stereographic
    1020             : // WKT: Polar_Stereographic
    1021             : //
    1022             : // Map parameters:
    1023             : //
    1024             : //    * straight_vertical_longitude_from_pole
    1025             : //    * latitude_of_projection_origin - Either +90. or -90.
    1026             : //    * Either standard_parallel or scale_factor_at_projection_origin
    1027             : //    * false_easting
    1028             : //    * false_northing
    1029             : 
    1030             : static const oNetcdfSRS_PP poPSmappings[] = {
    1031             :     /* {CF_PP_STD_PARALLEL_1, SRS_PP_LATITUDE_OF_ORIGIN}, */
    1032             :     /* {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR},   */
    1033             :     {CF_PP_VERT_LONG_FROM_POLE, SRS_PP_CENTRAL_MERIDIAN},
    1034             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
    1035             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
    1036             :     {nullptr, nullptr}};
    1037             : 
    1038             : // Rotated Pole
    1039             : //
    1040             : // grid_mapping_name = rotated_latitude_longitude
    1041             : // WKT: N/A
    1042             : //
    1043             : // Map parameters:
    1044             : //
    1045             : //    * grid_north_pole_latitude
    1046             : //    * grid_north_pole_longitude
    1047             : //    * north_pole_grid_longitude - This parameter is optional (default is 0.).
    1048             : 
    1049             : // No WKT equivalent
    1050             : 
    1051             : // Stereographic
    1052             : //
    1053             : // grid_mapping_name = stereographic
    1054             : // WKT: Stereographic (and/or Oblique_Stereographic??)
    1055             : //
    1056             : // Map parameters:
    1057             : //
    1058             : //    * longitude_of_projection_origin
    1059             : //    * latitude_of_projection_origin
    1060             : //    * scale_factor_at_projection_origin
    1061             : //    * false_easting
    1062             : //    * false_northing
    1063             : //
    1064             : // NB: see bug#4267 Stereographic vs. Oblique_Stereographic
    1065             : //
    1066             : static const oNetcdfSRS_PP poStMappings[] = {
    1067             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
    1068             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
    1069             :     {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR},
    1070             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
    1071             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
    1072             :     {nullptr, nullptr}};
    1073             : 
    1074             : // Transverse Mercator
    1075             : //
    1076             : // grid_mapping_name = transverse_mercator
    1077             : // WKT: Transverse_Mercator
    1078             : //
    1079             : // Map parameters:
    1080             : //
    1081             : //    * scale_factor_at_central_meridian
    1082             : //    * longitude_of_central_meridian
    1083             : //    * latitude_of_projection_origin
    1084             : //    * false_easting
    1085             : //    * false_northing
    1086             : //
    1087             : static const oNetcdfSRS_PP poTMMappings[] = {
    1088             :     {CF_PP_SCALE_FACTOR_MERIDIAN, SRS_PP_SCALE_FACTOR},
    1089             :     {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
    1090             :     {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
    1091             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
    1092             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
    1093             :     {nullptr, nullptr}};
    1094             : 
    1095             : // Vertical perspective
    1096             : //
    1097             : // grid_mapping_name = vertical_perspective
    1098             : // WKT: ???
    1099             : //
    1100             : // Map parameters:
    1101             : //
    1102             : //    * latitude_of_projection_origin
    1103             : //    * longitude_of_projection_origin
    1104             : //    * perspective_point_height
    1105             : //    * false_easting
    1106             : //    * false_northing
    1107             : //
    1108             : // TODO: see how to map this to OGR
    1109             : 
    1110             : static const oNetcdfSRS_PP poGEOSMappings[] = {
    1111             :     {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
    1112             :     {CF_PP_PERSPECTIVE_POINT_HEIGHT, SRS_PP_SATELLITE_HEIGHT},
    1113             :     {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
    1114             :     {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
    1115             :     /* { CF_PP_SWEEP_ANGLE_AXIS, .... } handled as a proj.4 extension */
    1116             :     {nullptr, nullptr}};
    1117             : 
    1118             : /* Mappings for various projections, including netcdf and GDAL projection names
    1119             :    and corresponding oNetcdfSRS_PP mapping struct.
    1120             :    A NULL mappings value means that the projection is not included in the CF
    1121             :    standard and the generic mapping (poGenericMappings) will be used. */
    1122             : typedef struct
    1123             : {
    1124             :     const char *CF_SRS;
    1125             :     const char *WKT_SRS;
    1126             :     const oNetcdfSRS_PP *mappings;
    1127             : } oNetcdfSRS_PT;
    1128             : 
    1129             : static const oNetcdfSRS_PT poNetcdfSRS_PT[] = {
    1130             :     {CF_PT_AEA, SRS_PT_ALBERS_CONIC_EQUAL_AREA, poAEAMappings},
    1131             :     {CF_PT_AE, SRS_PT_AZIMUTHAL_EQUIDISTANT, poAEMappings},
    1132             :     {"cassini_soldner", SRS_PT_CASSINI_SOLDNER, nullptr},
    1133             :     {CF_PT_LCEA, SRS_PT_CYLINDRICAL_EQUAL_AREA, poLCEAMappings},
    1134             :     {"eckert_iv", SRS_PT_ECKERT_IV, nullptr},
    1135             :     {"eckert_vi", SRS_PT_ECKERT_VI, nullptr},
    1136             :     {"equidistant_conic", SRS_PT_EQUIDISTANT_CONIC, nullptr},
    1137             :     {"equirectangular", SRS_PT_EQUIRECTANGULAR, nullptr},
    1138             :     {"gall_stereographic", SRS_PT_GALL_STEREOGRAPHIC, nullptr},
    1139             :     {CF_PT_GEOS, SRS_PT_GEOSTATIONARY_SATELLITE, poGEOSMappings},
    1140             :     {"goode_homolosine", SRS_PT_GOODE_HOMOLOSINE, nullptr},
    1141             :     {"gnomonic", SRS_PT_GNOMONIC, nullptr},
    1142             :     {"hotine_oblique_mercator", SRS_PT_HOTINE_OBLIQUE_MERCATOR, nullptr},
    1143             :     {"hotine_oblique_mercator_2P",
    1144             :      SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN, nullptr},
    1145             :     {"laborde_oblique_mercator", SRS_PT_LABORDE_OBLIQUE_MERCATOR, nullptr},
    1146             :     {CF_PT_LCC, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP, poLCC1SPMappings},
    1147             :     {CF_PT_LCC, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP, poLCC2SPMappings},
    1148             :     {CF_PT_LAEA, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA, poLAEAMappings},
    1149             :     {CF_PT_MERCATOR, SRS_PT_MERCATOR_1SP, poM1SPMappings},
    1150             :     {CF_PT_MERCATOR, SRS_PT_MERCATOR_2SP, poM2SPMappings},
    1151             :     {"miller_cylindrical", SRS_PT_MILLER_CYLINDRICAL, nullptr},
    1152             :     {"mollweide", SRS_PT_MOLLWEIDE, nullptr},
    1153             :     {"new_zealand_map_grid", SRS_PT_NEW_ZEALAND_MAP_GRID, nullptr},
    1154             :     /* for now map to STEREO, see bug #4267 */
    1155             :     {"oblique_stereographic", SRS_PT_OBLIQUE_STEREOGRAPHIC, nullptr},
    1156             :     /* {STEREO, SRS_PT_OBLIQUE_STEREOGRAPHIC, poStMappings },  */
    1157             :     {CF_PT_ORTHOGRAPHIC, SRS_PT_ORTHOGRAPHIC, poOrthoMappings},
    1158             :     {CF_PT_POLAR_STEREO, SRS_PT_POLAR_STEREOGRAPHIC, poPSmappings},
    1159             :     {"polyconic", SRS_PT_POLYCONIC, nullptr},
    1160             :     {"robinson", SRS_PT_ROBINSON, nullptr},
    1161             :     {"sinusoidal", SRS_PT_SINUSOIDAL, nullptr},
    1162             :     {CF_PT_STEREO, SRS_PT_STEREOGRAPHIC, poStMappings},
    1163             :     {"swiss_oblique_cylindrical", SRS_PT_SWISS_OBLIQUE_CYLINDRICAL, nullptr},
    1164             :     {CF_PT_TM, SRS_PT_TRANSVERSE_MERCATOR, poTMMappings},
    1165             :     {"TM_south_oriented", SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED, nullptr},
    1166             :     {nullptr, nullptr, nullptr},
    1167             : };
    1168             : 
    1169             : /* Write any needed projection attributes *
    1170             :  * poPROJCS: ptr to proj crd system
    1171             :  * pszProjection: name of projection system in GDAL WKT
    1172             :  *
    1173             :  * The function first looks for the oNetcdfSRS_PP mapping object
    1174             :  * that corresponds to the input projection name. If none is found
    1175             :  * the generic mapping is used.  In the case of specific mappings,
    1176             :  * the driver looks for each attribute listed in the mapping object
    1177             :  * and then looks up the value within the OGR_SRSNode. In the case
    1178             :  * of the generic mapping, the lookup is reversed (projection params,
    1179             :  * then mapping).  For more generic code, GDAL->NETCDF
    1180             :  * mappings and the associated value are saved in std::map objects.
    1181             :  */
    1182             : 
    1183             : // NOTE: modifications by ET to combine the specific and generic mappings.
    1184             : 
    1185             : static std::vector<std::pair<std::string, double>>
    1186          46 : NCDFGetProjAttribs(const OGR_SRSNode *poPROJCS, const char *pszProjection)
    1187             : {
    1188          46 :     const oNetcdfSRS_PP *poMap = nullptr;
    1189          46 :     int nMapIndex = -1;
    1190             : 
    1191             :     // Find the appropriate mapping.
    1192        1415 :     for (int iMap = 0; poNetcdfSRS_PT[iMap].WKT_SRS != nullptr; iMap++)
    1193             :     {
    1194        1415 :         if (EQUAL(pszProjection, poNetcdfSRS_PT[iMap].WKT_SRS))
    1195             :         {
    1196          46 :             nMapIndex = iMap;
    1197          46 :             poMap = poNetcdfSRS_PT[iMap].mappings;
    1198          46 :             break;
    1199             :         }
    1200             :     }
    1201             : 
    1202             :     // ET TODO if projection name is not found, should we do something special?
    1203          46 :     if (nMapIndex == -1)
    1204             :     {
    1205           0 :         CPLError(CE_Warning, CPLE_AppDefined,
    1206             :                  "projection name %s not found in the lookup tables!",
    1207             :                  pszProjection);
    1208             :     }
    1209             :     // If no mapping was found or assigned, set the generic one.
    1210          46 :     if (!poMap)
    1211             :     {
    1212           0 :         CPLError(CE_Warning, CPLE_AppDefined,
    1213             :                  "projection name %s in not part of the CF standard, "
    1214             :                  "will not be supported by CF!",
    1215             :                  pszProjection);
    1216           0 :         poMap = poGenericMappings;
    1217             :     }
    1218             : 
    1219             :     // Initialize local map objects.
    1220             : 
    1221             :     // Attribute <GDAL,NCDF> and Value <NCDF,value> mappings
    1222          92 :     std::map<std::string, std::string> oAttMap;
    1223         275 :     for (int iMap = 0; poMap[iMap].WKT_ATT != nullptr; iMap++)
    1224             :     {
    1225         229 :         oAttMap[poMap[iMap].WKT_ATT] = poMap[iMap].CF_ATT;
    1226             :     }
    1227             : 
    1228          46 :     const char *pszParamVal = nullptr;
    1229          92 :     std::map<std::string, double> oValMap;
    1230         595 :     for (int iChild = 0; iChild < poPROJCS->GetChildCount(); iChild++)
    1231             :     {
    1232         549 :         const OGR_SRSNode *poNode = poPROJCS->GetChild(iChild);
    1233         781 :         if (!EQUAL(poNode->GetValue(), "PARAMETER") ||
    1234         232 :             poNode->GetChildCount() != 2)
    1235         317 :             continue;
    1236         232 :         const char *pszParamStr = poNode->GetChild(0)->GetValue();
    1237         232 :         pszParamVal = poNode->GetChild(1)->GetValue();
    1238             : 
    1239         232 :         oValMap[pszParamStr] = CPLAtof(pszParamVal);
    1240             :     }
    1241             : 
    1242             :     // Results to write.
    1243          46 :     std::vector<std::pair<std::string, double>> oOutList;
    1244             : 
    1245             :     // Lookup mappings and fill output vector.
    1246          46 :     if (poMap != poGenericMappings)
    1247             :     {
    1248             :         // special case for PS (Polar Stereographic) grid.
    1249          46 :         if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
    1250             :         {
    1251           2 :             const double dfLat = oValMap[SRS_PP_LATITUDE_OF_ORIGIN];
    1252             : 
    1253           2 :             auto oScaleFactorIter = oValMap.find(SRS_PP_SCALE_FACTOR);
    1254           2 :             if (oScaleFactorIter != oValMap.end())
    1255             :             {
    1256             :                 // Polar Stereographic (variant A)
    1257           1 :                 const double dfScaleFactor = oScaleFactorIter->second;
    1258             :                 // dfLat should be +/- 90
    1259           1 :                 oOutList.push_back(
    1260           2 :                     std::make_pair(std::string(CF_PP_LAT_PROJ_ORIGIN), dfLat));
    1261           1 :                 oOutList.push_back(std::make_pair(
    1262           2 :                     std::string(CF_PP_SCALE_FACTOR_ORIGIN), dfScaleFactor));
    1263             :             }
    1264             :             else
    1265             :             {
    1266             :                 // Polar Stereographic (variant B)
    1267           1 :                 const double dfLatPole = (dfLat > 0) ? 90.0 : -90.0;
    1268           1 :                 oOutList.push_back(std::make_pair(
    1269           2 :                     std::string(CF_PP_LAT_PROJ_ORIGIN), dfLatPole));
    1270           1 :                 oOutList.push_back(
    1271           2 :                     std::make_pair(std::string(CF_PP_STD_PARALLEL), dfLat));
    1272             :             }
    1273             :         }
    1274             : 
    1275             :         // Specific mapping, loop over mapping values.
    1276         275 :         for (const auto &oAttIter : oAttMap)
    1277             :         {
    1278         229 :             const std::string &osGDALAtt = oAttIter.first;
    1279         229 :             const std::string &osNCDFAtt = oAttIter.second;
    1280         229 :             const auto oValIter = oValMap.find(osGDALAtt);
    1281             : 
    1282         229 :             if (oValIter != oValMap.end())
    1283             :             {
    1284         229 :                 double dfValue = oValIter->second;
    1285         229 :                 bool bWriteVal = true;
    1286             : 
    1287             :                 // special case for LCC-1SP
    1288             :                 //   See comments in netcdfdataset.h for this projection.
    1289         270 :                 if (EQUAL(SRS_PP_SCALE_FACTOR, osGDALAtt.c_str()) &&
    1290          41 :                     EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
    1291             :                 {
    1292             :                     // Default is to not write as it is not CF-1.
    1293           0 :                     bWriteVal = false;
    1294             :                     // Test if there is no standard_parallel1.
    1295           0 :                     if (!cpl::contains(oValMap, CF_PP_STD_PARALLEL_1))
    1296             :                     {
    1297             :                         // If scale factor != 1.0, write value for GDAL, but
    1298             :                         // this is not supported by CF-1.
    1299           0 :                         if (!CPLIsEqual(dfValue, 1.0))
    1300             :                         {
    1301           0 :                             CPLError(
    1302             :                                 CE_Failure, CPLE_NotSupported,
    1303             :                                 "NetCDF driver export of LCC-1SP with scale "
    1304             :                                 "factor != 1.0 and no standard_parallel1 is "
    1305             :                                 "not CF-1 (bug #3324).  Use the 2SP variant "
    1306             :                                 "which is supported by CF.");
    1307           0 :                             bWriteVal = true;
    1308             :                         }
    1309             :                         // Else copy standard_parallel1 from
    1310             :                         // latitude_of_origin, because scale_factor=1.0.
    1311             :                         else
    1312             :                         {
    1313             :                             const auto oValIter2 = oValMap.find(
    1314           0 :                                 std::string(SRS_PP_LATITUDE_OF_ORIGIN));
    1315           0 :                             if (oValIter2 != oValMap.end())
    1316             :                             {
    1317           0 :                                 oOutList.push_back(std::make_pair(
    1318           0 :                                     std::string(CF_PP_STD_PARALLEL_1),
    1319           0 :                                     oValIter2->second));
    1320             :                             }
    1321             :                             else
    1322             :                             {
    1323           0 :                                 CPLError(CE_Failure, CPLE_NotSupported,
    1324             :                                          "NetCDF driver export of LCC-1SP with "
    1325             :                                          "no standard_parallel1 "
    1326             :                                          "and no latitude_of_origin is not "
    1327             :                                          "supported (bug #3324).");
    1328             :                             }
    1329             :                         }
    1330             :                     }
    1331             :                 }
    1332         229 :                 if (bWriteVal)
    1333         229 :                     oOutList.push_back(std::make_pair(osNCDFAtt, dfValue));
    1334             :             }
    1335             : #ifdef NCDF_DEBUG
    1336             :             else
    1337             :             {
    1338             :                 CPLDebug("GDAL_netCDF", "NOT FOUND!");
    1339             :             }
    1340             : #endif
    1341             :         }
    1342             :     }
    1343             :     else
    1344             :     {
    1345             :         // Generic mapping, loop over projected values.
    1346           0 :         for (const auto &oValIter : oValMap)
    1347             :         {
    1348           0 :             const auto &osGDALAtt = oValIter.first;
    1349           0 :             const double dfValue = oValIter.second;
    1350             : 
    1351           0 :             const auto oAttIter = oAttMap.find(osGDALAtt);
    1352             : 
    1353           0 :             if (oAttIter != oAttMap.end())
    1354             :             {
    1355           0 :                 oOutList.push_back(std::make_pair(oAttIter->second, dfValue));
    1356             :             }
    1357             :             /* for SRS_PP_SCALE_FACTOR write 2 mappings */
    1358           0 :             else if (EQUAL(osGDALAtt.c_str(), SRS_PP_SCALE_FACTOR))
    1359             :             {
    1360           0 :                 oOutList.push_back(std::make_pair(
    1361           0 :                     std::string(CF_PP_SCALE_FACTOR_MERIDIAN), dfValue));
    1362           0 :                 oOutList.push_back(std::make_pair(
    1363           0 :                     std::string(CF_PP_SCALE_FACTOR_ORIGIN), dfValue));
    1364             :             }
    1365             :             /* if not found insert the GDAL name */
    1366             :             else
    1367             :             {
    1368           0 :                 oOutList.push_back(std::make_pair(osGDALAtt, dfValue));
    1369             :             }
    1370             :         }
    1371             :     }
    1372             : 
    1373          92 :     return oOutList;
    1374             : }
    1375             : 
    1376             : /************************************************************************/
    1377             : /*                           exportToCF1()                              */
    1378             : /************************************************************************/
    1379             : 
    1380             : /**
    1381             :  * \brief Export a CRS to netCDF CF-1 definitions.
    1382             :  *
    1383             :  * http://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings
    1384             :  *
    1385             :  * This function is the equivalent of the C function OSRExportToCF1().
    1386             :  *
    1387             :  * @param[out] ppszGridMappingName Pointer to the suggested name for the grid
    1388             :  * mapping variable. ppszGridMappingName may be nullptr.
    1389             :  * *ppszGridMappingName should be freed with CPLFree().
    1390             :  *
    1391             :  * @param[out] ppapszKeyValues Pointer to a null-terminated list of key/value pairs,
    1392             :  * to write into the grid mapping variable. ppapszKeyValues may be
    1393             :  * nullptr. *ppapszKeyValues should be freed with CSLDestroy()
    1394             :  * Values may be of type string, double or a list of 2 double values (comma
    1395             :  * separated).
    1396             :  *
    1397             :  * @param[out] ppszUnits Pointer to the value of the "units" attribute of the
    1398             :  * X/Y arrays. ppszGridMappingName may be nullptr. *ppszUnits should be freed with
    1399             :  * CPLFree().
    1400             :  *
    1401             :  * @param[in] papszOptions Options. Currently none supported
    1402             :  *
    1403             :  * @return OGRERR_NONE on success or an error code in case of failure.
    1404             :  * @since 3.9
    1405             :  */
    1406             : 
    1407             : OGRErr
    1408         223 : OGRSpatialReference::exportToCF1(char **ppszGridMappingName,
    1409             :                                  char ***ppapszKeyValues, char **ppszUnits,
    1410             :                                  CPL_UNUSED CSLConstList papszOptions) const
    1411             : {
    1412         223 :     if (ppszGridMappingName)
    1413         128 :         *ppszGridMappingName = nullptr;
    1414             : 
    1415         223 :     if (ppapszKeyValues)
    1416         132 :         *ppapszKeyValues = nullptr;
    1417             : 
    1418         223 :     if (ppszUnits)
    1419          41 :         *ppszUnits = nullptr;
    1420             : 
    1421         223 :     if (ppszGridMappingName || ppapszKeyValues)
    1422             :     {
    1423         132 :         bool bWriteWkt = true;
    1424             : 
    1425             :         struct Value
    1426             :         {
    1427             :             std::string key{};
    1428             :             std::string valueStr{};
    1429             :             std::vector<double> doubles{};
    1430             :         };
    1431             : 
    1432         132 :         std::vector<Value> oParams;
    1433             : 
    1434             :         const auto addParamString =
    1435         792 :             [&oParams](const char *key, const char *value)
    1436             :         {
    1437         792 :             Value v;
    1438         396 :             v.key = key;
    1439         396 :             v.valueStr = value;
    1440         396 :             oParams.emplace_back(std::move(v));
    1441         396 :         };
    1442             : 
    1443        1246 :         const auto addParamDouble = [&oParams](const char *key, double value)
    1444             :         {
    1445        1246 :             Value v;
    1446         623 :             v.key = key;
    1447         623 :             v.doubles.push_back(value);
    1448         623 :             oParams.emplace_back(std::move(v));
    1449         623 :         };
    1450             : 
    1451             :         const auto addParam2Double =
    1452           6 :             [&oParams](const char *key, double value1, double value2)
    1453             :         {
    1454           6 :             Value v;
    1455           3 :             v.key = key;
    1456           3 :             v.doubles.push_back(value1);
    1457           3 :             v.doubles.push_back(value2);
    1458           3 :             oParams.emplace_back(std::move(v));
    1459           3 :         };
    1460             : 
    1461         132 :         std::string osCFProjection;
    1462         132 :         if (IsProjected())
    1463             :         {
    1464             :             // Write CF-1.5 compliant Projected attributes.
    1465             : 
    1466          46 :             const OGR_SRSNode *poPROJCS = GetAttrNode("PROJCS");
    1467          46 :             if (poPROJCS == nullptr)
    1468           0 :                 return OGRERR_FAILURE;
    1469          46 :             const char *pszProjName = GetAttrValue("PROJECTION");
    1470          46 :             if (pszProjName == nullptr)
    1471           0 :                 return OGRERR_FAILURE;
    1472             : 
    1473             :             // Basic Projection info (grid_mapping and datum).
    1474        1415 :             for (int i = 0; poNetcdfSRS_PT[i].WKT_SRS != nullptr; i++)
    1475             :             {
    1476        1415 :                 if (EQUAL(poNetcdfSRS_PT[i].WKT_SRS, pszProjName))
    1477             :                 {
    1478          46 :                     osCFProjection = poNetcdfSRS_PT[i].CF_SRS;
    1479          46 :                     break;
    1480             :                 }
    1481             :             }
    1482          46 :             if (osCFProjection.empty())
    1483           0 :                 return OGRERR_FAILURE;
    1484             : 
    1485          46 :             addParamString(CF_GRD_MAPPING_NAME, osCFProjection.c_str());
    1486             : 
    1487             :             // Various projection attributes.
    1488             :             // PDS: keep in sync with SetProjection function
    1489          92 :             auto oOutList = NCDFGetProjAttribs(poPROJCS, pszProjName);
    1490             : 
    1491             :             /* Write all the values that were found */
    1492          46 :             double dfStdP[2] = {0, 0};
    1493          46 :             bool bFoundStdP1 = false;
    1494          46 :             bool bFoundStdP2 = false;
    1495         279 :             for (const auto &it : oOutList)
    1496             :             {
    1497         233 :                 const char *pszParamVal = it.first.c_str();
    1498         233 :                 double dfValue = it.second;
    1499             :                 /* Handle the STD_PARALLEL attrib */
    1500         233 :                 if (EQUAL(pszParamVal, CF_PP_STD_PARALLEL_1))
    1501             :                 {
    1502           3 :                     bFoundStdP1 = true;
    1503           3 :                     dfStdP[0] = dfValue;
    1504             :                 }
    1505         230 :                 else if (EQUAL(pszParamVal, CF_PP_STD_PARALLEL_2))
    1506             :                 {
    1507           3 :                     bFoundStdP2 = true;
    1508           3 :                     dfStdP[1] = dfValue;
    1509             :                 }
    1510             :                 else
    1511             :                 {
    1512         227 :                     addParamDouble(pszParamVal, dfValue);
    1513             :                 }
    1514             :             }
    1515             :             /* Now write the STD_PARALLEL attrib */
    1516          46 :             if (bFoundStdP1)
    1517             :             {
    1518             :                 /* one value  */
    1519           3 :                 if (!bFoundStdP2)
    1520             :                 {
    1521           0 :                     addParamDouble(CF_PP_STD_PARALLEL, dfStdP[0]);
    1522             :                 }
    1523             :                 else
    1524             :                 {
    1525             :                     // Two values.
    1526           3 :                     addParam2Double(CF_PP_STD_PARALLEL, dfStdP[0], dfStdP[1]);
    1527             :                 }
    1528             :             }
    1529             : 
    1530          46 :             if (EQUAL(pszProjName, SRS_PT_GEOSTATIONARY_SATELLITE))
    1531             :             {
    1532             :                 const char *pszPredefProj4 =
    1533           0 :                     GetExtension(GetRoot()->GetValue(), "PROJ4", nullptr);
    1534           0 :                 const char *pszSweepAxisAngle =
    1535           0 :                     (pszPredefProj4 != nullptr &&
    1536           0 :                      strstr(pszPredefProj4, "+sweep=x"))
    1537             :                         ? "x"
    1538             :                         : "y";
    1539           0 :                 addParamString(CF_PP_SWEEP_ANGLE_AXIS, pszSweepAxisAngle);
    1540             :             }
    1541             :         }
    1542          86 :         else if (IsDerivedGeographic())
    1543             :         {
    1544           0 :             const OGR_SRSNode *poConversion = GetAttrNode("DERIVINGCONVERSION");
    1545           0 :             if (poConversion == nullptr)
    1546           0 :                 return OGRERR_FAILURE;
    1547           0 :             const char *pszMethod = GetAttrValue("METHOD");
    1548           0 :             if (pszMethod == nullptr)
    1549           0 :                 return OGRERR_FAILURE;
    1550             : 
    1551           0 :             std::map<std::string, double> oValMap;
    1552           0 :             for (int iChild = 0; iChild < poConversion->GetChildCount();
    1553             :                  iChild++)
    1554             :             {
    1555           0 :                 const OGR_SRSNode *poNode = poConversion->GetChild(iChild);
    1556           0 :                 if (!EQUAL(poNode->GetValue(), "PARAMETER") ||
    1557           0 :                     poNode->GetChildCount() <= 2)
    1558           0 :                     continue;
    1559           0 :                 const char *pszParamStr = poNode->GetChild(0)->GetValue();
    1560           0 :                 const char *pszParamVal = poNode->GetChild(1)->GetValue();
    1561           0 :                 oValMap[pszParamStr] = CPLAtof(pszParamVal);
    1562             :             }
    1563             : 
    1564           0 :             constexpr const char *ROTATED_POLE_VAR_NAME = "rotated_pole";
    1565           0 :             if (EQUAL(pszMethod, "PROJ ob_tran o_proj=longlat"))
    1566             :             {
    1567             :                 // Not enough interoperable to be written as WKT
    1568           0 :                 bWriteWkt = false;
    1569             : 
    1570           0 :                 const double dfLon0 = oValMap["lon_0"];
    1571           0 :                 const double dfLonp = oValMap["o_lon_p"];
    1572           0 :                 const double dfLatp = oValMap["o_lat_p"];
    1573             : 
    1574           0 :                 osCFProjection = ROTATED_POLE_VAR_NAME;
    1575           0 :                 addParamString(CF_GRD_MAPPING_NAME,
    1576             :                                CF_PT_ROTATED_LATITUDE_LONGITUDE);
    1577           0 :                 addParamDouble(CF_PP_GRID_NORTH_POLE_LONGITUDE, dfLon0 - 180);
    1578           0 :                 addParamDouble(CF_PP_GRID_NORTH_POLE_LATITUDE, dfLatp);
    1579           0 :                 addParamDouble(CF_PP_NORTH_POLE_GRID_LONGITUDE, dfLonp);
    1580             :             }
    1581           0 :             else if (EQUAL(pszMethod, "Pole rotation (netCDF CF convention)"))
    1582             :             {
    1583             :                 // Not enough interoperable to be written as WKT
    1584           0 :                 bWriteWkt = false;
    1585             : 
    1586             :                 const double dfGridNorthPoleLat =
    1587           0 :                     oValMap["Grid north pole latitude (netCDF CF convention)"];
    1588             :                 const double dfGridNorthPoleLong =
    1589           0 :                     oValMap["Grid north pole longitude (netCDF CF convention)"];
    1590             :                 const double dfNorthPoleGridLong =
    1591           0 :                     oValMap["North pole grid longitude (netCDF CF convention)"];
    1592             : 
    1593           0 :                 osCFProjection = ROTATED_POLE_VAR_NAME;
    1594           0 :                 addParamString(CF_GRD_MAPPING_NAME,
    1595             :                                CF_PT_ROTATED_LATITUDE_LONGITUDE);
    1596           0 :                 addParamDouble(CF_PP_GRID_NORTH_POLE_LONGITUDE,
    1597             :                                dfGridNorthPoleLong);
    1598           0 :                 addParamDouble(CF_PP_GRID_NORTH_POLE_LATITUDE,
    1599             :                                dfGridNorthPoleLat);
    1600           0 :                 addParamDouble(CF_PP_NORTH_POLE_GRID_LONGITUDE,
    1601             :                                dfNorthPoleGridLong);
    1602             :             }
    1603           0 :             else if (EQUAL(pszMethod, "Pole rotation (GRIB convention)"))
    1604             :             {
    1605             :                 // Not enough interoperable to be written as WKT
    1606           0 :                 bWriteWkt = false;
    1607             : 
    1608             :                 const double dfLatSouthernPole =
    1609           0 :                     oValMap["Latitude of the southern pole (GRIB convention)"];
    1610             :                 const double dfLonSouthernPole =
    1611           0 :                     oValMap["Longitude of the southern pole (GRIB convention)"];
    1612             :                 const double dfAxisRotation =
    1613           0 :                     oValMap["Axis rotation (GRIB convention)"];
    1614             : 
    1615           0 :                 const double dfLon0 = dfLonSouthernPole;
    1616           0 :                 const double dfLonp = dfAxisRotation == 0 ? 0 : -dfAxisRotation;
    1617           0 :                 const double dfLatp =
    1618           0 :                     dfLatSouthernPole == 0 ? 0 : -dfLatSouthernPole;
    1619             : 
    1620           0 :                 osCFProjection = ROTATED_POLE_VAR_NAME;
    1621           0 :                 addParamString(CF_GRD_MAPPING_NAME,
    1622             :                                CF_PT_ROTATED_LATITUDE_LONGITUDE);
    1623           0 :                 addParamDouble(CF_PP_GRID_NORTH_POLE_LONGITUDE, dfLon0 - 180);
    1624           0 :                 addParamDouble(CF_PP_GRID_NORTH_POLE_LATITUDE, dfLatp);
    1625           0 :                 addParamDouble(CF_PP_NORTH_POLE_GRID_LONGITUDE, dfLonp);
    1626             :             }
    1627             :             else
    1628             :             {
    1629           0 :                 CPLError(CE_Failure, CPLE_NotSupported,
    1630             :                          "Unsupported method for DerivedGeographicCRS: %s",
    1631             :                          pszMethod);
    1632           0 :                 return OGRERR_FAILURE;
    1633             :             }
    1634             :         }
    1635             :         else
    1636             :         {
    1637             :             // Write CF-1.5 compliant Geographics attributes.
    1638             :             // Note: WKT information will not be preserved (e.g. WGS84).
    1639          86 :             osCFProjection = "crs";
    1640          86 :             addParamString(CF_GRD_MAPPING_NAME, CF_PT_LATITUDE_LONGITUDE);
    1641             :         }
    1642             : 
    1643         132 :         constexpr const char *CF_LNG_NAME = "long_name";
    1644         132 :         addParamString(CF_LNG_NAME, "CRS definition");
    1645             : 
    1646             :         // Write CF-1.5 compliant common attributes.
    1647             : 
    1648             :         // DATUM information.
    1649         132 :         addParamDouble(CF_PP_LONG_PRIME_MERIDIAN, GetPrimeMeridian());
    1650         132 :         addParamDouble(CF_PP_SEMI_MAJOR_AXIS, GetSemiMajor());
    1651         132 :         addParamDouble(CF_PP_INVERSE_FLATTENING, GetInvFlattening());
    1652             : 
    1653         132 :         if (bWriteWkt)
    1654             :         {
    1655         132 :             char *pszSpatialRef = nullptr;
    1656         132 :             exportToWkt(&pszSpatialRef);
    1657         132 :             if (pszSpatialRef && pszSpatialRef[0])
    1658             :             {
    1659         132 :                 addParamString(NCDF_CRS_WKT, pszSpatialRef);
    1660             :             }
    1661         132 :             CPLFree(pszSpatialRef);
    1662             :         }
    1663             : 
    1664         132 :         if (ppszGridMappingName)
    1665         128 :             *ppszGridMappingName = CPLStrdup(osCFProjection.c_str());
    1666             : 
    1667         132 :         if (ppapszKeyValues)
    1668             :         {
    1669         132 :             CPLStringList aosKeyValues;
    1670        1154 :             for (const auto &param : oParams)
    1671             :             {
    1672        1022 :                 if (!param.valueStr.empty())
    1673             :                 {
    1674             :                     aosKeyValues.AddNameValue(param.key.c_str(),
    1675         396 :                                               param.valueStr.c_str());
    1676             :                 }
    1677             :                 else
    1678             :                 {
    1679        1252 :                     std::string osVal;
    1680        1255 :                     for (const double dfVal : param.doubles)
    1681             :                     {
    1682         629 :                         if (!osVal.empty())
    1683           3 :                             osVal += ',';
    1684         629 :                         osVal += CPLSPrintf("%.17g", dfVal);
    1685             :                     }
    1686         626 :                     aosKeyValues.AddNameValue(param.key.c_str(), osVal.c_str());
    1687             :                 }
    1688             :             }
    1689         132 :             *ppapszKeyValues = aosKeyValues.StealList();
    1690             :         }
    1691             :     }
    1692             : 
    1693         223 :     if (ppszUnits)
    1694             :     {
    1695          41 :         const char *pszUnits = nullptr;
    1696          41 :         const char *pszUnitsToWrite = "";
    1697             : 
    1698          41 :         const double dfUnits = GetLinearUnits(&pszUnits);
    1699          41 :         if (fabs(dfUnits - 1.0) < 1e-15 || pszUnits == nullptr ||
    1700           1 :             EQUAL(pszUnits, "m") || EQUAL(pszUnits, "metre"))
    1701             :         {
    1702          40 :             pszUnitsToWrite = "m";
    1703             :         }
    1704           1 :         else if (fabs(dfUnits - 1000.0) < 1e-15)
    1705             :         {
    1706           0 :             pszUnitsToWrite = "km";
    1707             :         }
    1708           1 :         else if (fabs(dfUnits - CPLAtof(SRS_UL_US_FOOT_CONV)) < 1e-15 ||
    1709           1 :                  EQUAL(pszUnits, SRS_UL_US_FOOT) ||
    1710           0 :                  EQUAL(pszUnits, "US survey foot"))
    1711             :         {
    1712           1 :             pszUnitsToWrite = "US_survey_foot";
    1713             :         }
    1714          41 :         if (pszUnitsToWrite)
    1715          41 :             *ppszUnits = CPLStrdup(pszUnitsToWrite);
    1716             :     }
    1717             : 
    1718         223 :     return OGRERR_NONE;
    1719             : }
    1720             : 
    1721             : /************************************************************************/
    1722             : /*                          OSRExportToCF1()                            */
    1723             : /************************************************************************/
    1724             : /**
    1725             :  * \brief Export a CRS to netCDF CF-1 definitions.
    1726             :  *
    1727             :  * This function is the same as OGRSpatialReference::exportToCF1().
    1728             :  */
    1729           5 : OGRErr OSRExportToCF1(OGRSpatialReferenceH hSRS, char **ppszGridMappingName,
    1730             :                       char ***ppapszKeyValues, char **ppszUnits,
    1731             :                       CSLConstList papszOptions)
    1732             : 
    1733             : {
    1734           5 :     VALIDATE_POINTER1(hSRS, "OSRExportToCF1", OGRERR_FAILURE);
    1735             : 
    1736           5 :     return OGRSpatialReference::FromHandle(hSRS)->exportToCF1(
    1737           5 :         ppszGridMappingName, ppapszKeyValues, ppszUnits, papszOptions);
    1738             : }

Generated by: LCOV version 1.14