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