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 = static_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 48 : OGRErr OGRSpatialReference::importFromCF1(CSLConstList papszKeyValues,
162 : const char *pszUnits)
163 : {
164 : // Import from "spatial_ref" or "crs_wkt" attributes in priority
165 48 : const char *pszWKT = CSLFetchNameValue(papszKeyValues, NCDF_SPATIAL_REF);
166 48 : if (!pszWKT)
167 : {
168 38 : pszWKT = CSLFetchNameValue(papszKeyValues, NCDF_CRS_WKT);
169 : }
170 48 : if (pszWKT)
171 11 : return importFromWkt(pszWKT);
172 :
173 : const char *pszGridMappingName =
174 37 : 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 42 : if (pszGridMappingName == nullptr &&
180 5 : CSLFetchNameValue(papszKeyValues, CF_PP_GRID_NORTH_POLE_LONGITUDE) !=
181 : nullptr)
182 : {
183 3 : pszGridMappingName = CF_PT_ROTATED_LATITUDE_LONGITUDE;
184 : }
185 :
186 37 : if (pszGridMappingName == nullptr)
187 2 : 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") || EQUAL(pszUnits, "meters") ||
758 9 : EQUAL(pszUnits, "metres"))
759 : {
760 12 : SetLinearUnits("metre", 1.0);
761 12 : SetAuthority("PROJCS|UNIT", "EPSG", 9001);
762 : }
763 9 : else if (EQUAL(pszUnits, "km"))
764 : {
765 6 : SetLinearUnits("kilometre", 1000.0);
766 6 : SetAuthority("PROJCS|UNIT", "EPSG", 9036);
767 : }
768 3 : else if (EQUAL(pszUnits, "US_survey_foot") ||
769 3 : EQUAL(pszUnits, "US_survey_feet"))
770 : {
771 0 : SetLinearUnits("US survey foot", CPLAtof(SRS_UL_US_FOOT_CONV));
772 0 : SetAuthority("PROJCS|UNIT", "EPSG", 9003);
773 : }
774 : else
775 : {
776 3 : CPLError(CE_Warning, CPLE_AppDefined,
777 : "Unhandled X/Y axis unit %s. SRS will ignore "
778 : "axis unit and be likely wrong.",
779 : pszUnits);
780 : }
781 : }
782 :
783 35 : return OGRERR_NONE;
784 : }
785 :
786 : /* -------------------------------------------------------------------- */
787 : /* CF-1 to GDAL mappings */
788 : /* -------------------------------------------------------------------- */
789 :
790 : /* Following are a series of mappings from CF-1 convention parameters
791 : * for each projection, to the equivalent in OGC WKT used internally by GDAL.
792 : * See: http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.5/apf.html
793 : */
794 :
795 : /* A struct allowing us to map between GDAL(OGC WKT) and CF-1 attributes */
796 : typedef struct
797 : {
798 : const char *CF_ATT;
799 : const char *WKT_ATT;
800 : // TODO: mappings may need default values, like scale factor?
801 : // double defval;
802 : } oNetcdfSRS_PP;
803 :
804 : // default mappings, for the generic case
805 : /* These 'generic' mappings are based on what was previously in the
806 : poNetCDFSRS struct. They will be used as a fallback in case none
807 : of the others match (i.e. you are exporting a projection that has
808 : no CF-1 equivalent).
809 : They are not used for known CF-1 projections since there is not a
810 : unique 2-way projection-independent
811 : mapping between OGC WKT params and CF-1 ones: it varies per-projection.
812 : */
813 :
814 : static const oNetcdfSRS_PP poGenericMappings[] = {
815 : /* scale_factor is handled as a special case, write 2 values */
816 : {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
817 : {CF_PP_STD_PARALLEL_2, SRS_PP_STANDARD_PARALLEL_2},
818 : {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
819 : {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_LONGITUDE_OF_CENTER},
820 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_LONGITUDE_OF_ORIGIN},
821 : // Multiple mappings to LAT_PROJ_ORIGIN
822 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
823 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
824 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
825 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
826 : {nullptr, nullptr},
827 : };
828 :
829 : // Albers equal area
830 : //
831 : // grid_mapping_name = albers_conical_equal_area
832 : // WKT: Albers_Conic_Equal_Area
833 : // EPSG:9822
834 : //
835 : // Map parameters:
836 : //
837 : // * standard_parallel - There may be 1 or 2 values.
838 : // * longitude_of_central_meridian
839 : // * latitude_of_projection_origin
840 : // * false_easting
841 : // * false_northing
842 : //
843 : static const oNetcdfSRS_PP poAEAMappings[] = {
844 : {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
845 : {CF_PP_STD_PARALLEL_2, SRS_PP_STANDARD_PARALLEL_2},
846 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
847 : {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_LONGITUDE_OF_CENTER},
848 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
849 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
850 : {nullptr, nullptr}};
851 :
852 : // Azimuthal equidistant
853 : //
854 : // grid_mapping_name = azimuthal_equidistant
855 : // WKT: Azimuthal_Equidistant
856 : //
857 : // Map parameters:
858 : //
859 : // * longitude_of_projection_origin
860 : // * latitude_of_projection_origin
861 : // * false_easting
862 : // * false_northing
863 : //
864 : static const oNetcdfSRS_PP poAEMappings[] = {
865 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
866 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_LONGITUDE_OF_CENTER},
867 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
868 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
869 : {nullptr, nullptr}};
870 :
871 : // Lambert azimuthal equal area
872 : //
873 : // grid_mapping_name = lambert_azimuthal_equal_area
874 : // WKT: Lambert_Azimuthal_Equal_Area
875 : //
876 : // Map parameters:
877 : //
878 : // * longitude_of_projection_origin
879 : // * latitude_of_projection_origin
880 : // * false_easting
881 : // * false_northing
882 : //
883 : static const oNetcdfSRS_PP poLAEAMappings[] = {
884 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER},
885 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_LONGITUDE_OF_CENTER},
886 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
887 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
888 : {nullptr, nullptr}};
889 :
890 : // Lambert conformal
891 : //
892 : // grid_mapping_name = lambert_conformal_conic
893 : // WKT: Lambert_Conformal_Conic_1SP / Lambert_Conformal_Conic_2SP
894 : //
895 : // Map parameters:
896 : //
897 : // * standard_parallel - There may be 1 or 2 values.
898 : // * longitude_of_central_meridian
899 : // * latitude_of_projection_origin
900 : // * false_easting
901 : // * false_northing
902 : //
903 : // See
904 : // http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html
905 :
906 : // Lambert conformal conic - 1SP
907 : /* See bug # 3324
908 : It seems that the missing scale factor can be computed from
909 : standard_parallel1 and latitude_of_projection_origin. If both are equal (the
910 : common case) then scale factor=1, else use Snyder eq. 15-4. We save in the
911 : WKT standard_parallel1 for export to CF, but do not export scale factor. If a
912 : WKT has a scale factor != 1 and no standard_parallel1 then export is not CF,
913 : but we output scale factor for compat. is there a formula for that?
914 : */
915 : static const oNetcdfSRS_PP poLCC1SPMappings[] = {
916 : {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
917 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
918 : {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
919 : {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR}, /* special case */
920 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
921 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
922 : {nullptr, nullptr}};
923 :
924 : // Lambert conformal conic - 2SP
925 : static const oNetcdfSRS_PP poLCC2SPMappings[] = {
926 : {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
927 : {CF_PP_STD_PARALLEL_2, SRS_PP_STANDARD_PARALLEL_2},
928 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
929 : {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
930 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
931 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
932 : {nullptr, nullptr}};
933 :
934 : // Lambert cylindrical equal area
935 : //
936 : // grid_mapping_name = lambert_cylindrical_equal_area
937 : // WKT: Cylindrical_Equal_Area
938 : // EPSG:9834 (Spherical) and EPSG:9835
939 : //
940 : // Map parameters:
941 : //
942 : // * longitude_of_central_meridian
943 : // * either standard_parallel or scale_factor_at_projection_origin
944 : // * false_easting
945 : // * false_northing
946 : //
947 : // NB: CF-1 specifies a 'scale_factor_at_projection' alternative
948 : // to std_parallel ... but no reference to this in EPSG/remotesensing.org
949 : // ignore for now.
950 : //
951 : static const oNetcdfSRS_PP poLCEAMappings[] = {
952 : {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
953 : {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
954 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
955 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
956 : {nullptr, nullptr}};
957 :
958 : // Latitude-Longitude
959 : //
960 : // grid_mapping_name = latitude_longitude
961 : //
962 : // Map parameters:
963 : //
964 : // * None
965 : //
966 : // NB: handled as a special case - !isProjected()
967 :
968 : // Mercator
969 : //
970 : // grid_mapping_name = mercator
971 : // WKT: Mercator_1SP / Mercator_2SP
972 : //
973 : // Map parameters:
974 : //
975 : // * longitude_of_projection_origin
976 : // * either standard_parallel or scale_factor_at_projection_origin
977 : // * false_easting
978 : // * false_northing
979 :
980 : // Mercator 1 Standard Parallel (EPSG:9804)
981 : static const oNetcdfSRS_PP poM1SPMappings[] = {
982 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
983 : // LAT_PROJ_ORIGIN is always equator (0) in CF-1
984 : {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR},
985 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
986 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
987 : {nullptr, nullptr}};
988 :
989 : // Mercator 2 Standard Parallel
990 : static const oNetcdfSRS_PP poM2SPMappings[] = {
991 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
992 : {CF_PP_STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1},
993 : // From best understanding of this projection, only
994 : // actually specify one SP - it is the same N/S of equator.
995 : // {CF_PP_STD_PARALLEL_2, SRS_PP_LATITUDE_OF_ORIGIN},
996 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
997 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
998 : {nullptr, nullptr}};
999 :
1000 : // Orthographic
1001 : // grid_mapping_name = orthographic
1002 : // WKT: Orthographic
1003 : //
1004 : // Map parameters:
1005 : //
1006 : // * longitude_of_projection_origin
1007 : // * latitude_of_projection_origin
1008 : // * false_easting
1009 : // * false_northing
1010 : //
1011 : static const oNetcdfSRS_PP poOrthoMappings[] = {
1012 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
1013 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
1014 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
1015 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
1016 : {nullptr, nullptr}};
1017 :
1018 : // Polar stereographic
1019 : //
1020 : // grid_mapping_name = polar_stereographic
1021 : // WKT: Polar_Stereographic
1022 : //
1023 : // Map parameters:
1024 : //
1025 : // * straight_vertical_longitude_from_pole
1026 : // * latitude_of_projection_origin - Either +90. or -90.
1027 : // * Either standard_parallel or scale_factor_at_projection_origin
1028 : // * false_easting
1029 : // * false_northing
1030 :
1031 : static const oNetcdfSRS_PP poPSmappings[] = {
1032 : /* {CF_PP_STD_PARALLEL_1, SRS_PP_LATITUDE_OF_ORIGIN}, */
1033 : /* {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR}, */
1034 : {CF_PP_VERT_LONG_FROM_POLE, SRS_PP_CENTRAL_MERIDIAN},
1035 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
1036 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
1037 : {nullptr, nullptr}};
1038 :
1039 : // Rotated Pole
1040 : //
1041 : // grid_mapping_name = rotated_latitude_longitude
1042 : // WKT: N/A
1043 : //
1044 : // Map parameters:
1045 : //
1046 : // * grid_north_pole_latitude
1047 : // * grid_north_pole_longitude
1048 : // * north_pole_grid_longitude - This parameter is optional (default is 0.).
1049 :
1050 : // No WKT equivalent
1051 :
1052 : // Stereographic
1053 : //
1054 : // grid_mapping_name = stereographic
1055 : // WKT: Stereographic (and/or Oblique_Stereographic??)
1056 : //
1057 : // Map parameters:
1058 : //
1059 : // * longitude_of_projection_origin
1060 : // * latitude_of_projection_origin
1061 : // * scale_factor_at_projection_origin
1062 : // * false_easting
1063 : // * false_northing
1064 : //
1065 : // NB: see bug#4267 Stereographic vs. Oblique_Stereographic
1066 : //
1067 : static const oNetcdfSRS_PP poStMappings[] = {
1068 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
1069 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
1070 : {CF_PP_SCALE_FACTOR_ORIGIN, SRS_PP_SCALE_FACTOR},
1071 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
1072 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
1073 : {nullptr, nullptr}};
1074 :
1075 : // Transverse Mercator
1076 : //
1077 : // grid_mapping_name = transverse_mercator
1078 : // WKT: Transverse_Mercator
1079 : //
1080 : // Map parameters:
1081 : //
1082 : // * scale_factor_at_central_meridian
1083 : // * longitude_of_central_meridian
1084 : // * latitude_of_projection_origin
1085 : // * false_easting
1086 : // * false_northing
1087 : //
1088 : static const oNetcdfSRS_PP poTMMappings[] = {
1089 : {CF_PP_SCALE_FACTOR_MERIDIAN, SRS_PP_SCALE_FACTOR},
1090 : {CF_PP_LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN},
1091 : {CF_PP_LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN},
1092 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
1093 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
1094 : {nullptr, nullptr}};
1095 :
1096 : // Vertical perspective
1097 : //
1098 : // grid_mapping_name = vertical_perspective
1099 : // WKT: ???
1100 : //
1101 : // Map parameters:
1102 : //
1103 : // * latitude_of_projection_origin
1104 : // * longitude_of_projection_origin
1105 : // * perspective_point_height
1106 : // * false_easting
1107 : // * false_northing
1108 : //
1109 : // TODO: see how to map this to OGR
1110 :
1111 : static const oNetcdfSRS_PP poGEOSMappings[] = {
1112 : {CF_PP_LON_PROJ_ORIGIN, SRS_PP_CENTRAL_MERIDIAN},
1113 : {CF_PP_PERSPECTIVE_POINT_HEIGHT, SRS_PP_SATELLITE_HEIGHT},
1114 : {CF_PP_FALSE_EASTING, SRS_PP_FALSE_EASTING},
1115 : {CF_PP_FALSE_NORTHING, SRS_PP_FALSE_NORTHING},
1116 : /* { CF_PP_SWEEP_ANGLE_AXIS, .... } handled as a proj.4 extension */
1117 : {nullptr, nullptr}};
1118 :
1119 : /* Mappings for various projections, including netcdf and GDAL projection names
1120 : and corresponding oNetcdfSRS_PP mapping struct.
1121 : A NULL mappings value means that the projection is not included in the CF
1122 : standard and the generic mapping (poGenericMappings) will be used. */
1123 : typedef struct
1124 : {
1125 : const char *CF_SRS;
1126 : const char *WKT_SRS;
1127 : const oNetcdfSRS_PP *mappings;
1128 : } oNetcdfSRS_PT;
1129 :
1130 : static const oNetcdfSRS_PT poNetcdfSRS_PT[] = {
1131 : {CF_PT_AEA, SRS_PT_ALBERS_CONIC_EQUAL_AREA, poAEAMappings},
1132 : {CF_PT_AE, SRS_PT_AZIMUTHAL_EQUIDISTANT, poAEMappings},
1133 : {"cassini_soldner", SRS_PT_CASSINI_SOLDNER, nullptr},
1134 : {CF_PT_LCEA, SRS_PT_CYLINDRICAL_EQUAL_AREA, poLCEAMappings},
1135 : {"eckert_iv", SRS_PT_ECKERT_IV, nullptr},
1136 : {"eckert_vi", SRS_PT_ECKERT_VI, nullptr},
1137 : {"equidistant_conic", SRS_PT_EQUIDISTANT_CONIC, nullptr},
1138 : {"equirectangular", SRS_PT_EQUIRECTANGULAR, nullptr},
1139 : {"gall_stereographic", SRS_PT_GALL_STEREOGRAPHIC, nullptr},
1140 : {CF_PT_GEOS, SRS_PT_GEOSTATIONARY_SATELLITE, poGEOSMappings},
1141 : {"goode_homolosine", SRS_PT_GOODE_HOMOLOSINE, nullptr},
1142 : {"gnomonic", SRS_PT_GNOMONIC, nullptr},
1143 : {"hotine_oblique_mercator", SRS_PT_HOTINE_OBLIQUE_MERCATOR, nullptr},
1144 : {"hotine_oblique_mercator_2P",
1145 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN, nullptr},
1146 : {"laborde_oblique_mercator", SRS_PT_LABORDE_OBLIQUE_MERCATOR, nullptr},
1147 : {CF_PT_LCC, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP, poLCC1SPMappings},
1148 : {CF_PT_LCC, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP, poLCC2SPMappings},
1149 : {CF_PT_LAEA, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA, poLAEAMappings},
1150 : {CF_PT_MERCATOR, SRS_PT_MERCATOR_1SP, poM1SPMappings},
1151 : {CF_PT_MERCATOR, SRS_PT_MERCATOR_2SP, poM2SPMappings},
1152 : {"miller_cylindrical", SRS_PT_MILLER_CYLINDRICAL, nullptr},
1153 : {"mollweide", SRS_PT_MOLLWEIDE, nullptr},
1154 : {"new_zealand_map_grid", SRS_PT_NEW_ZEALAND_MAP_GRID, nullptr},
1155 : /* for now map to STEREO, see bug #4267 */
1156 : {"oblique_stereographic", SRS_PT_OBLIQUE_STEREOGRAPHIC, nullptr},
1157 : /* {STEREO, SRS_PT_OBLIQUE_STEREOGRAPHIC, poStMappings }, */
1158 : {CF_PT_ORTHOGRAPHIC, SRS_PT_ORTHOGRAPHIC, poOrthoMappings},
1159 : {CF_PT_POLAR_STEREO, SRS_PT_POLAR_STEREOGRAPHIC, poPSmappings},
1160 : {"polyconic", SRS_PT_POLYCONIC, nullptr},
1161 : {"robinson", SRS_PT_ROBINSON, nullptr},
1162 : {"sinusoidal", SRS_PT_SINUSOIDAL, nullptr},
1163 : {CF_PT_STEREO, SRS_PT_STEREOGRAPHIC, poStMappings},
1164 : {"swiss_oblique_cylindrical", SRS_PT_SWISS_OBLIQUE_CYLINDRICAL, nullptr},
1165 : {CF_PT_TM, SRS_PT_TRANSVERSE_MERCATOR, poTMMappings},
1166 : {"TM_south_oriented", SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED, nullptr},
1167 : {nullptr, nullptr, nullptr},
1168 : };
1169 :
1170 : /* Write any needed projection attributes *
1171 : * poPROJCS: ptr to proj crd system
1172 : * pszProjection: name of projection system in GDAL WKT
1173 : *
1174 : * The function first looks for the oNetcdfSRS_PP mapping object
1175 : * that corresponds to the input projection name. If none is found
1176 : * the generic mapping is used. In the case of specific mappings,
1177 : * the driver looks for each attribute listed in the mapping object
1178 : * and then looks up the value within the OGR_SRSNode. In the case
1179 : * of the generic mapping, the lookup is reversed (projection params,
1180 : * then mapping). For more generic code, GDAL->NETCDF
1181 : * mappings and the associated value are saved in std::map objects.
1182 : */
1183 :
1184 : // NOTE: modifications by ET to combine the specific and generic mappings.
1185 :
1186 : static std::vector<std::pair<std::string, double>>
1187 61 : NCDFGetProjAttribs(const OGR_SRSNode *poPROJCS, const char *pszProjection)
1188 : {
1189 61 : const oNetcdfSRS_PP *poMap = nullptr;
1190 61 : int nMapIndex = -1;
1191 :
1192 : // Find the appropriate mapping.
1193 1895 : for (int iMap = 0; poNetcdfSRS_PT[iMap].WKT_SRS != nullptr; iMap++)
1194 : {
1195 1895 : if (EQUAL(pszProjection, poNetcdfSRS_PT[iMap].WKT_SRS))
1196 : {
1197 61 : nMapIndex = iMap;
1198 61 : poMap = poNetcdfSRS_PT[iMap].mappings;
1199 61 : break;
1200 : }
1201 : }
1202 :
1203 : // ET TODO if projection name is not found, should we do something special?
1204 61 : if (nMapIndex == -1)
1205 : {
1206 0 : CPLError(CE_Warning, CPLE_AppDefined,
1207 : "projection name %s not found in the lookup tables!",
1208 : pszProjection);
1209 : }
1210 : // If no mapping was found or assigned, set the generic one.
1211 61 : if (!poMap)
1212 : {
1213 0 : CPLError(CE_Warning, CPLE_AppDefined,
1214 : "projection name %s in not part of the CF standard, "
1215 : "will not be supported by CF!",
1216 : pszProjection);
1217 0 : poMap = poGenericMappings;
1218 : }
1219 :
1220 : // Initialize local map objects.
1221 :
1222 : // Attribute <GDAL,NCDF> and Value <NCDF,value> mappings
1223 122 : std::map<std::string, std::string> oAttMap;
1224 365 : for (int iMap = 0; poMap[iMap].WKT_ATT != nullptr; iMap++)
1225 : {
1226 304 : oAttMap[poMap[iMap].WKT_ATT] = poMap[iMap].CF_ATT;
1227 : }
1228 :
1229 61 : const char *pszParamVal = nullptr;
1230 122 : std::map<std::string, double> oValMap;
1231 790 : for (int iChild = 0; iChild < poPROJCS->GetChildCount(); iChild++)
1232 : {
1233 729 : const OGR_SRSNode *poNode = poPROJCS->GetChild(iChild);
1234 1036 : if (!EQUAL(poNode->GetValue(), "PARAMETER") ||
1235 307 : poNode->GetChildCount() != 2)
1236 422 : continue;
1237 307 : const char *pszParamStr = poNode->GetChild(0)->GetValue();
1238 307 : pszParamVal = poNode->GetChild(1)->GetValue();
1239 :
1240 307 : oValMap[pszParamStr] = CPLAtof(pszParamVal);
1241 : }
1242 :
1243 : // Results to write.
1244 61 : std::vector<std::pair<std::string, double>> oOutList;
1245 :
1246 : // Lookup mappings and fill output vector.
1247 61 : if (poMap != poGenericMappings)
1248 : {
1249 : // special case for PS (Polar Stereographic) grid.
1250 61 : if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
1251 : {
1252 2 : const double dfLat = oValMap[SRS_PP_LATITUDE_OF_ORIGIN];
1253 :
1254 2 : auto oScaleFactorIter = oValMap.find(SRS_PP_SCALE_FACTOR);
1255 2 : if (oScaleFactorIter != oValMap.end())
1256 : {
1257 : // Polar Stereographic (variant A)
1258 1 : const double dfScaleFactor = oScaleFactorIter->second;
1259 : // dfLat should be +/- 90
1260 1 : oOutList.push_back(
1261 2 : std::make_pair(std::string(CF_PP_LAT_PROJ_ORIGIN), dfLat));
1262 1 : oOutList.push_back(std::make_pair(
1263 2 : std::string(CF_PP_SCALE_FACTOR_ORIGIN), dfScaleFactor));
1264 : }
1265 : else
1266 : {
1267 : // Polar Stereographic (variant B)
1268 1 : const double dfLatPole = (dfLat > 0) ? 90.0 : -90.0;
1269 1 : oOutList.push_back(std::make_pair(
1270 2 : std::string(CF_PP_LAT_PROJ_ORIGIN), dfLatPole));
1271 1 : oOutList.push_back(
1272 2 : std::make_pair(std::string(CF_PP_STD_PARALLEL), dfLat));
1273 : }
1274 : }
1275 :
1276 : // Specific mapping, loop over mapping values.
1277 365 : for (const auto &oAttIter : oAttMap)
1278 : {
1279 304 : const std::string &osGDALAtt = oAttIter.first;
1280 304 : const std::string &osNCDFAtt = oAttIter.second;
1281 304 : const auto oValIter = oValMap.find(osGDALAtt);
1282 :
1283 304 : if (oValIter != oValMap.end())
1284 : {
1285 304 : double dfValue = oValIter->second;
1286 304 : bool bWriteVal = true;
1287 :
1288 : // special case for LCC-1SP
1289 : // See comments in netcdfdataset.h for this projection.
1290 360 : if (EQUAL(SRS_PP_SCALE_FACTOR, osGDALAtt.c_str()) &&
1291 56 : EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
1292 : {
1293 : // Default is to not write as it is not CF-1.
1294 0 : bWriteVal = false;
1295 : // Test if there is no standard_parallel1.
1296 0 : if (!cpl::contains(oValMap, CF_PP_STD_PARALLEL_1))
1297 : {
1298 : // If scale factor != 1.0, write value for GDAL, but
1299 : // this is not supported by CF-1.
1300 0 : if (!CPLIsEqual(dfValue, 1.0))
1301 : {
1302 0 : CPLError(
1303 : CE_Failure, CPLE_NotSupported,
1304 : "NetCDF driver export of LCC-1SP with scale "
1305 : "factor != 1.0 and no standard_parallel1 is "
1306 : "not CF-1 (bug #3324). Use the 2SP variant "
1307 : "which is supported by CF.");
1308 0 : bWriteVal = true;
1309 : }
1310 : // Else copy standard_parallel1 from
1311 : // latitude_of_origin, because scale_factor=1.0.
1312 : else
1313 : {
1314 : const auto oValIter2 = oValMap.find(
1315 0 : std::string(SRS_PP_LATITUDE_OF_ORIGIN));
1316 0 : if (oValIter2 != oValMap.end())
1317 : {
1318 0 : oOutList.push_back(std::make_pair(
1319 0 : std::string(CF_PP_STD_PARALLEL_1),
1320 0 : oValIter2->second));
1321 : }
1322 : else
1323 : {
1324 0 : CPLError(CE_Failure, CPLE_NotSupported,
1325 : "NetCDF driver export of LCC-1SP with "
1326 : "no standard_parallel1 "
1327 : "and no latitude_of_origin is not "
1328 : "supported (bug #3324).");
1329 : }
1330 : }
1331 : }
1332 : }
1333 304 : if (bWriteVal)
1334 304 : oOutList.push_back(std::make_pair(osNCDFAtt, dfValue));
1335 : }
1336 : #ifdef NCDF_DEBUG
1337 : else
1338 : {
1339 : CPLDebug("GDAL_netCDF", "NOT FOUND!");
1340 : }
1341 : #endif
1342 : }
1343 : }
1344 : else
1345 : {
1346 : // Generic mapping, loop over projected values.
1347 0 : for (const auto &oValIter : oValMap)
1348 : {
1349 0 : const auto &osGDALAtt = oValIter.first;
1350 0 : const double dfValue = oValIter.second;
1351 :
1352 0 : const auto oAttIter = oAttMap.find(osGDALAtt);
1353 :
1354 0 : if (oAttIter != oAttMap.end())
1355 : {
1356 0 : oOutList.push_back(std::make_pair(oAttIter->second, dfValue));
1357 : }
1358 : /* for SRS_PP_SCALE_FACTOR write 2 mappings */
1359 0 : else if (EQUAL(osGDALAtt.c_str(), SRS_PP_SCALE_FACTOR))
1360 : {
1361 0 : oOutList.push_back(std::make_pair(
1362 0 : std::string(CF_PP_SCALE_FACTOR_MERIDIAN), dfValue));
1363 0 : oOutList.push_back(std::make_pair(
1364 0 : std::string(CF_PP_SCALE_FACTOR_ORIGIN), dfValue));
1365 : }
1366 : /* if not found insert the GDAL name */
1367 : else
1368 : {
1369 0 : oOutList.push_back(std::make_pair(osGDALAtt, dfValue));
1370 : }
1371 : }
1372 : }
1373 :
1374 122 : return oOutList;
1375 : }
1376 :
1377 : /************************************************************************/
1378 : /* exportToCF1() */
1379 : /************************************************************************/
1380 :
1381 : /**
1382 : * \brief Export a CRS to netCDF CF-1 definitions.
1383 : *
1384 : * http://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings
1385 : *
1386 : * This function is the equivalent of the C function OSRExportToCF1().
1387 : *
1388 : * @param[out] ppszGridMappingName Pointer to the suggested name for the grid
1389 : * mapping variable. ppszGridMappingName may be nullptr.
1390 : * *ppszGridMappingName should be freed with CPLFree().
1391 : *
1392 : * @param[out] ppapszKeyValues Pointer to a null-terminated list of key/value pairs,
1393 : * to write into the grid mapping variable. ppapszKeyValues may be
1394 : * nullptr. *ppapszKeyValues should be freed with CSLDestroy()
1395 : * Values may be of type string, double or a list of 2 double values (comma
1396 : * separated).
1397 : *
1398 : * @param[out] ppszUnits Pointer to the value of the "units" attribute of the
1399 : * X/Y arrays. ppszGridMappingName may be nullptr. *ppszUnits should be freed with
1400 : * CPLFree().
1401 : *
1402 : * @param[in] papszOptions Options. Currently none supported
1403 : *
1404 : * @return OGRERR_NONE on success or an error code in case of failure.
1405 : * @since 3.9
1406 : */
1407 :
1408 : OGRErr
1409 271 : OGRSpatialReference::exportToCF1(char **ppszGridMappingName,
1410 : char ***ppapszKeyValues, char **ppszUnits,
1411 : CPL_UNUSED CSLConstList papszOptions) const
1412 : {
1413 271 : if (ppszGridMappingName)
1414 143 : *ppszGridMappingName = nullptr;
1415 :
1416 271 : if (ppapszKeyValues)
1417 147 : *ppapszKeyValues = nullptr;
1418 :
1419 271 : if (ppszUnits)
1420 52 : *ppszUnits = nullptr;
1421 :
1422 271 : if (ppszGridMappingName || ppapszKeyValues)
1423 : {
1424 147 : bool bWriteWkt = true;
1425 :
1426 : struct Value
1427 : {
1428 : std::string key{};
1429 : std::string valueStr{};
1430 : std::vector<double> doubles{};
1431 : };
1432 :
1433 147 : std::vector<Value> oParams;
1434 :
1435 : const auto addParamString =
1436 882 : [&oParams](const char *key, const char *value)
1437 : {
1438 882 : Value v;
1439 441 : v.key = key;
1440 441 : v.valueStr = value;
1441 441 : oParams.emplace_back(std::move(v));
1442 441 : };
1443 :
1444 1486 : const auto addParamDouble = [&oParams](const char *key, double value)
1445 : {
1446 1486 : Value v;
1447 743 : v.key = key;
1448 743 : v.doubles.push_back(value);
1449 743 : oParams.emplace_back(std::move(v));
1450 743 : };
1451 :
1452 : const auto addParam2Double =
1453 6 : [&oParams](const char *key, double value1, double value2)
1454 : {
1455 6 : Value v;
1456 3 : v.key = key;
1457 3 : v.doubles.push_back(value1);
1458 3 : v.doubles.push_back(value2);
1459 3 : oParams.emplace_back(std::move(v));
1460 3 : };
1461 :
1462 147 : std::string osCFProjection;
1463 147 : if (IsProjected())
1464 : {
1465 : // Write CF-1.5 compliant Projected attributes.
1466 :
1467 61 : const OGR_SRSNode *poPROJCS = GetAttrNode("PROJCS");
1468 61 : if (poPROJCS == nullptr)
1469 0 : return OGRERR_FAILURE;
1470 61 : const char *pszProjName = GetAttrValue("PROJECTION");
1471 61 : if (pszProjName == nullptr)
1472 0 : return OGRERR_FAILURE;
1473 :
1474 : // Basic Projection info (grid_mapping and datum).
1475 1895 : for (int i = 0; poNetcdfSRS_PT[i].WKT_SRS != nullptr; i++)
1476 : {
1477 1895 : if (EQUAL(poNetcdfSRS_PT[i].WKT_SRS, pszProjName))
1478 : {
1479 61 : osCFProjection = poNetcdfSRS_PT[i].CF_SRS;
1480 61 : break;
1481 : }
1482 : }
1483 61 : if (osCFProjection.empty())
1484 0 : return OGRERR_FAILURE;
1485 :
1486 61 : addParamString(CF_GRD_MAPPING_NAME, osCFProjection.c_str());
1487 :
1488 : // Various projection attributes.
1489 : // PDS: keep in sync with SetProjection function
1490 122 : auto oOutList = NCDFGetProjAttribs(poPROJCS, pszProjName);
1491 :
1492 : /* Write all the values that were found */
1493 61 : double dfStdP[2] = {0, 0};
1494 61 : bool bFoundStdP1 = false;
1495 61 : bool bFoundStdP2 = false;
1496 369 : for (const auto &it : oOutList)
1497 : {
1498 308 : const char *pszParamVal = it.first.c_str();
1499 308 : double dfValue = it.second;
1500 : /* Handle the STD_PARALLEL attrib */
1501 308 : if (EQUAL(pszParamVal, CF_PP_STD_PARALLEL_1))
1502 : {
1503 3 : bFoundStdP1 = true;
1504 3 : dfStdP[0] = dfValue;
1505 : }
1506 305 : else if (EQUAL(pszParamVal, CF_PP_STD_PARALLEL_2))
1507 : {
1508 3 : bFoundStdP2 = true;
1509 3 : dfStdP[1] = dfValue;
1510 : }
1511 : else
1512 : {
1513 302 : addParamDouble(pszParamVal, dfValue);
1514 : }
1515 : }
1516 : /* Now write the STD_PARALLEL attrib */
1517 61 : if (bFoundStdP1)
1518 : {
1519 : /* one value */
1520 3 : if (!bFoundStdP2)
1521 : {
1522 0 : addParamDouble(CF_PP_STD_PARALLEL, dfStdP[0]);
1523 : }
1524 : else
1525 : {
1526 : // Two values.
1527 3 : addParam2Double(CF_PP_STD_PARALLEL, dfStdP[0], dfStdP[1]);
1528 : }
1529 : }
1530 :
1531 61 : if (EQUAL(pszProjName, SRS_PT_GEOSTATIONARY_SATELLITE))
1532 : {
1533 : const char *pszPredefProj4 =
1534 0 : GetExtension(GetRoot()->GetValue(), "PROJ4", nullptr);
1535 0 : const char *pszSweepAxisAngle =
1536 0 : (pszPredefProj4 != nullptr &&
1537 0 : strstr(pszPredefProj4, "+sweep=x"))
1538 : ? "x"
1539 : : "y";
1540 0 : addParamString(CF_PP_SWEEP_ANGLE_AXIS, pszSweepAxisAngle);
1541 : }
1542 : }
1543 86 : else if (IsDerivedGeographic())
1544 : {
1545 0 : const OGR_SRSNode *poConversion = GetAttrNode("DERIVINGCONVERSION");
1546 0 : if (poConversion == nullptr)
1547 0 : return OGRERR_FAILURE;
1548 0 : const char *pszMethod = GetAttrValue("METHOD");
1549 0 : if (pszMethod == nullptr)
1550 0 : return OGRERR_FAILURE;
1551 :
1552 0 : std::map<std::string, double> oValMap;
1553 0 : for (int iChild = 0; iChild < poConversion->GetChildCount();
1554 : iChild++)
1555 : {
1556 0 : const OGR_SRSNode *poNode = poConversion->GetChild(iChild);
1557 0 : if (!EQUAL(poNode->GetValue(), "PARAMETER") ||
1558 0 : poNode->GetChildCount() <= 2)
1559 0 : continue;
1560 0 : const char *pszParamStr = poNode->GetChild(0)->GetValue();
1561 0 : const char *pszParamVal = poNode->GetChild(1)->GetValue();
1562 0 : oValMap[pszParamStr] = CPLAtof(pszParamVal);
1563 : }
1564 :
1565 0 : constexpr const char *ROTATED_POLE_VAR_NAME = "rotated_pole";
1566 0 : if (EQUAL(pszMethod, "PROJ ob_tran o_proj=longlat"))
1567 : {
1568 : // Not enough interoperable to be written as WKT
1569 0 : bWriteWkt = false;
1570 :
1571 0 : const double dfLon0 = oValMap["lon_0"];
1572 0 : const double dfLonp = oValMap["o_lon_p"];
1573 0 : const double dfLatp = oValMap["o_lat_p"];
1574 :
1575 0 : osCFProjection = ROTATED_POLE_VAR_NAME;
1576 0 : addParamString(CF_GRD_MAPPING_NAME,
1577 : CF_PT_ROTATED_LATITUDE_LONGITUDE);
1578 0 : addParamDouble(CF_PP_GRID_NORTH_POLE_LONGITUDE, dfLon0 - 180);
1579 0 : addParamDouble(CF_PP_GRID_NORTH_POLE_LATITUDE, dfLatp);
1580 0 : addParamDouble(CF_PP_NORTH_POLE_GRID_LONGITUDE, dfLonp);
1581 : }
1582 0 : else if (EQUAL(pszMethod, "Pole rotation (netCDF CF convention)"))
1583 : {
1584 : // Not enough interoperable to be written as WKT
1585 0 : bWriteWkt = false;
1586 :
1587 : const double dfGridNorthPoleLat =
1588 0 : oValMap["Grid north pole latitude (netCDF CF convention)"];
1589 : const double dfGridNorthPoleLong =
1590 0 : oValMap["Grid north pole longitude (netCDF CF convention)"];
1591 : const double dfNorthPoleGridLong =
1592 0 : oValMap["North pole grid longitude (netCDF CF convention)"];
1593 :
1594 0 : osCFProjection = ROTATED_POLE_VAR_NAME;
1595 0 : addParamString(CF_GRD_MAPPING_NAME,
1596 : CF_PT_ROTATED_LATITUDE_LONGITUDE);
1597 0 : addParamDouble(CF_PP_GRID_NORTH_POLE_LONGITUDE,
1598 : dfGridNorthPoleLong);
1599 0 : addParamDouble(CF_PP_GRID_NORTH_POLE_LATITUDE,
1600 : dfGridNorthPoleLat);
1601 0 : addParamDouble(CF_PP_NORTH_POLE_GRID_LONGITUDE,
1602 : dfNorthPoleGridLong);
1603 : }
1604 0 : else if (EQUAL(pszMethod, "Pole rotation (GRIB convention)"))
1605 : {
1606 : // Not enough interoperable to be written as WKT
1607 0 : bWriteWkt = false;
1608 :
1609 : const double dfLatSouthernPole =
1610 0 : oValMap["Latitude of the southern pole (GRIB convention)"];
1611 : const double dfLonSouthernPole =
1612 0 : oValMap["Longitude of the southern pole (GRIB convention)"];
1613 : const double dfAxisRotation =
1614 0 : oValMap["Axis rotation (GRIB convention)"];
1615 :
1616 0 : const double dfLon0 = dfLonSouthernPole;
1617 0 : const double dfLonp = dfAxisRotation == 0 ? 0 : -dfAxisRotation;
1618 0 : const double dfLatp =
1619 0 : dfLatSouthernPole == 0 ? 0 : -dfLatSouthernPole;
1620 :
1621 0 : osCFProjection = ROTATED_POLE_VAR_NAME;
1622 0 : addParamString(CF_GRD_MAPPING_NAME,
1623 : CF_PT_ROTATED_LATITUDE_LONGITUDE);
1624 0 : addParamDouble(CF_PP_GRID_NORTH_POLE_LONGITUDE, dfLon0 - 180);
1625 0 : addParamDouble(CF_PP_GRID_NORTH_POLE_LATITUDE, dfLatp);
1626 0 : addParamDouble(CF_PP_NORTH_POLE_GRID_LONGITUDE, dfLonp);
1627 : }
1628 : else
1629 : {
1630 0 : CPLError(CE_Failure, CPLE_NotSupported,
1631 : "Unsupported method for DerivedGeographicCRS: %s",
1632 : pszMethod);
1633 0 : return OGRERR_FAILURE;
1634 : }
1635 : }
1636 : else
1637 : {
1638 : // Write CF-1.5 compliant Geographics attributes.
1639 : // Note: WKT information will not be preserved (e.g. WGS84).
1640 86 : osCFProjection = "crs";
1641 86 : addParamString(CF_GRD_MAPPING_NAME, CF_PT_LATITUDE_LONGITUDE);
1642 : }
1643 :
1644 147 : constexpr const char *CF_LNG_NAME = "long_name";
1645 147 : addParamString(CF_LNG_NAME, "CRS definition");
1646 :
1647 : // Write CF-1.5 compliant common attributes.
1648 :
1649 : // DATUM information.
1650 147 : addParamDouble(CF_PP_LONG_PRIME_MERIDIAN, GetPrimeMeridian());
1651 147 : addParamDouble(CF_PP_SEMI_MAJOR_AXIS, GetSemiMajor());
1652 147 : addParamDouble(CF_PP_INVERSE_FLATTENING, GetInvFlattening());
1653 :
1654 147 : if (bWriteWkt)
1655 : {
1656 147 : char *pszSpatialRef = nullptr;
1657 147 : exportToWkt(&pszSpatialRef);
1658 147 : if (pszSpatialRef && pszSpatialRef[0])
1659 : {
1660 147 : addParamString(NCDF_CRS_WKT, pszSpatialRef);
1661 : }
1662 147 : CPLFree(pszSpatialRef);
1663 : }
1664 :
1665 147 : if (ppszGridMappingName)
1666 143 : *ppszGridMappingName = CPLStrdup(osCFProjection.c_str());
1667 :
1668 147 : if (ppapszKeyValues)
1669 : {
1670 147 : CPLStringList aosKeyValues;
1671 1334 : for (const auto ¶m : oParams)
1672 : {
1673 1187 : if (!param.valueStr.empty())
1674 : {
1675 : aosKeyValues.AddNameValue(param.key.c_str(),
1676 441 : param.valueStr.c_str());
1677 : }
1678 : else
1679 : {
1680 1492 : std::string osVal;
1681 1495 : for (const double dfVal : param.doubles)
1682 : {
1683 749 : if (!osVal.empty())
1684 3 : osVal += ',';
1685 749 : osVal += CPLSPrintf("%.17g", dfVal);
1686 : }
1687 746 : aosKeyValues.AddNameValue(param.key.c_str(), osVal.c_str());
1688 : }
1689 : }
1690 147 : *ppapszKeyValues = aosKeyValues.StealList();
1691 : }
1692 : }
1693 :
1694 271 : if (ppszUnits)
1695 : {
1696 52 : const char *pszUnits = nullptr;
1697 52 : const char *pszUnitsToWrite = "";
1698 :
1699 52 : const double dfUnits = GetLinearUnits(&pszUnits);
1700 52 : if (fabs(dfUnits - 1.0) < 1e-15 || pszUnits == nullptr ||
1701 1 : EQUAL(pszUnits, "m") || EQUAL(pszUnits, "metre"))
1702 : {
1703 51 : pszUnitsToWrite = "m";
1704 : }
1705 1 : else if (fabs(dfUnits - 1000.0) < 1e-15)
1706 : {
1707 0 : pszUnitsToWrite = "km";
1708 : }
1709 1 : else if (fabs(dfUnits - CPLAtof(SRS_UL_US_FOOT_CONV)) < 1e-15 ||
1710 1 : EQUAL(pszUnits, SRS_UL_US_FOOT) ||
1711 0 : EQUAL(pszUnits, "US survey foot"))
1712 : {
1713 1 : pszUnitsToWrite = "US_survey_foot";
1714 : }
1715 52 : if (pszUnitsToWrite)
1716 52 : *ppszUnits = CPLStrdup(pszUnitsToWrite);
1717 : }
1718 :
1719 271 : return OGRERR_NONE;
1720 : }
1721 :
1722 : /************************************************************************/
1723 : /* OSRExportToCF1() */
1724 : /************************************************************************/
1725 : /**
1726 : * \brief Export a CRS to netCDF CF-1 definitions.
1727 : *
1728 : * This function is the same as OGRSpatialReference::exportToCF1().
1729 : */
1730 5 : OGRErr OSRExportToCF1(OGRSpatialReferenceH hSRS, char **ppszGridMappingName,
1731 : char ***ppapszKeyValues, char **ppszUnits,
1732 : CSLConstList papszOptions)
1733 :
1734 : {
1735 5 : VALIDATE_POINTER1(hSRS, "OSRExportToCF1", OGRERR_FAILURE);
1736 :
1737 5 : return OGRSpatialReference::FromHandle(hSRS)->exportToCF1(
1738 5 : ppszGridMappingName, ppapszKeyValues, ppszUnits, papszOptions);
1739 : }
|
