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