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