Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
projection.cc
Go to the documentation of this file.
1/* Copyright (C) 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20#include <cstdlib> // strtol
21#include <cmath> // std::pow, std::sqrt, std::fabs
22#include <string>
23#include <vector>
24#include <tuple>
25
26#include "pism/util/projection.hh"
27#include "pism/util/VariableMetadata.hh"
28#include "pism/util/error_handling.hh"
29#include "pism/util/io/File.hh"
30#include "pism/util/io/io_helpers.hh"
31#include "pism/util/Grid.hh"
32#include "pism/util/array/Scalar.hh"
33#include "pism/util/array/Array3D.hh"
34
35#include "pism/pism_config.hh"
36#include "pism/util/pism_utilities.hh"
37
38#if (Pism_USE_PROJ==1)
39#include "pism/util/Proj.hh"
40#endif
41
42namespace pism {
43
44MappingInfo::MappingInfo(const std::string &mapping_variable_name, units::System::Ptr unit_system)
45 : cf_mapping(mapping_variable_name, unit_system) {
46 // empty
47}
48
49MappingInfo::MappingInfo(const VariableMetadata &mapping_variable, const std::string &proj_string_)
50 : cf_mapping(mapping_variable), proj_string(proj_string_) {
51 // empty
52}
53
54int parse_epsg(const std::string &proj_string) {
55 try {
56 int auth_len = 5; // length of "epsg:"
57 std::string::size_type position = std::string::npos;
58 for (const auto &auth : { "epsg:", "EPSG:" }) {
59 position = proj_string.find(auth);
60 if (position != std::string::npos) {
61 break;
62 }
63 }
64
65 if (position == std::string::npos) {
66 throw RuntimeError(PISM_ERROR_LOCATION, "could not find 'epsg:' or 'EPSG:'");
67 }
68
69
70 {
71 auto epsg_string = proj_string.substr(position + auth_len);
72 char *endptr = NULL;
73 long int result = strtol(epsg_string.c_str(), &endptr, 10);
74 if (endptr == epsg_string.c_str()) {
75 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't parse %s (expected an integer)",
76 epsg_string.c_str());
77 }
78 return (int)result;
79 }
80
81 } catch (RuntimeError &e) {
82 e.add_context("parsing the EPSG code in the PROJ string '%s'",
83 proj_string.c_str());
84 throw;
85 }
86}
87
88//! @brief Return CF-Convention "mapping" variable corresponding to an EPSG code specified in a
89//! PROJ string.
90VariableMetadata epsg_to_cf(units::System::Ptr system, const std::string &proj_string) {
91 VariableMetadata mapping("mapping", system);
92
93 int epsg = parse_epsg(proj_string);
94
95 switch (epsg) {
96 case 3413:
97 mapping["latitude_of_projection_origin"] = {90.0};
98 mapping["scale_factor_at_projection_origin"] = {1.0};
99 mapping["straight_vertical_longitude_from_pole"] = {-45.0};
100 mapping["standard_parallel"] = {70.0};
101 mapping["false_northing"] = {0.0};
102 mapping["grid_mapping_name"] = "polar_stereographic";
103 mapping["false_easting"] = {0.0};
104 break;
105 case 3031:
106 mapping["latitude_of_projection_origin"] = {-90.0};
107 mapping["scale_factor_at_projection_origin"] = {1.0};
108 mapping["straight_vertical_longitude_from_pole"] = {0.0};
109 mapping["standard_parallel"] = {-71.0};
110 mapping["false_northing"] = {0.0};
111 mapping["grid_mapping_name"] = "polar_stereographic";
112 mapping["false_easting"] = {0.0};
113 break;
114 case 3057:
115 mapping["grid_mapping_name"] = "lambert_conformal_conic" ;
116 mapping["longitude_of_central_meridian"] = {-19.} ;
117 mapping["false_easting"] = {500000.} ;
118 mapping["false_northing"] = {500000.} ;
119 mapping["latitude_of_projection_origin"] = {65.} ;
120 mapping["standard_parallel"] = {64.25, 65.75} ;
121 mapping["long_name"] = "CRS definition" ;
122 mapping["longitude_of_prime_meridian"] = {0.} ;
123 mapping["semi_major_axis"] = {6378137.} ;
124 mapping["inverse_flattening"] = {298.257222101} ;
125 break;
126 case 5936:
127 mapping["latitude_of_projection_origin"] = {90.0};
128 mapping["scale_factor_at_projection_origin"] = {1.0};
129 mapping["straight_vertical_longitude_from_pole"] = {-150.0};
130 mapping["standard_parallel"] = {90.0};
131 mapping["false_northing"] = {2000000.0};
132 mapping["grid_mapping_name"] = "polar_stereographic";
133 mapping["false_easting"] = {2000000.0};
134 break;
135 case 26710:
136 mapping["longitude_of_central_meridian"] = {-123.0};
137 mapping["false_easting"] = {500000.0};
138 mapping["false_northing"] = {0.0};
139 mapping["grid_mapping_name"] = "transverse_mercator";
140 mapping["inverse_flattening"] = {294.978698213898};
141 mapping["latitude_of_projection_origin"] = {0.0};
142 mapping["scale_factor_at_central_meridian"] = {0.9996};
143 mapping["semi_major_axis"] = {6378206.4};
144 mapping["unit"] = "metre";
145 break;
146 default:
147 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "unknown EPSG code '%d' in PROJ string '%s'",
148 (int)epsg, proj_string.c_str());
149 }
150
151 return mapping;
152}
153
154void check_consistency_epsg(const VariableMetadata &cf_mapping, const std::string &proj_string) {
155
156 VariableMetadata epsg_mapping = epsg_to_cf(cf_mapping.unit_system(), proj_string);
157
158 // All CF grid mapping variables have to have the "grid_mapping_name" attribute.
159 bool mapping_is_empty = not cf_mapping.has_attribute("grid_mapping_name");
160 bool epsg_mapping_is_empty = not epsg_mapping.has_attribute("grid_mapping_name");
161
162 if (mapping_is_empty and epsg_mapping_is_empty) {
163 // empty mapping variables are equivalent
164 return;
165 }
166
167 // Check if the `cf_mapping` variable matches the EPSG code in `proj_string`.
168 //
169 // Check strings.
170 for (const auto &s : epsg_mapping.all_strings()) {
171 if (not cf_mapping.has_attribute(s.first)) {
173 "inconsistent metadata:\n"
174 "PROJ string \"%s\" requires %s = \"%s\",\n"
175 "but the mapping variable has no attribute \"%s\".",
176 proj_string.c_str(), s.first.c_str(), s.second.c_str(),
177 s.first.c_str());
178 }
179
180 std::string string = cf_mapping[s.first];
181
182 if (not(string == s.second)) {
184 "inconsistent metadata:\n"
185 "%s requires %s = \"%s\",\n"
186 "but the mapping variable has %s = \"%s\".",
187 proj_string.c_str(), s.first.c_str(), s.second.c_str(),
188 s.first.c_str(), string.c_str());
189 }
190 }
191
192 // Check doubles
193 const double eps = 1e-12;
194 for (auto d : epsg_mapping.all_doubles()) {
195 if (not cf_mapping.has_attribute(d.first)) {
197 "inconsistent metadata:\n"
198 "%s requires %s = %f,\n"
199 "but the mapping variable has no attribute \"%s\".",
200 proj_string.c_str(), d.first.c_str(), d.second[0],
201 d.first.c_str());
202 }
203
204 double value = cf_mapping.get_number(d.first);
205
206 if (std::fabs(value - d.second[0]) > eps) {
208 "inconsistent metadata:\n"
209 "%s requires %s = %f,\n"
210 "but the mapping variable has %s = %f.",
211 proj_string.c_str(), d.first.c_str(), d.second[0],
212 d.first.c_str(), value);
213 }
214 }
215}
216
217std::string grid_name(const pism::File &file, const std::string &variable_name,
218 pism::units::System::Ptr sys, bool piecewise_constant) {
219
220 std::vector<std::string> dimensions = file.dimensions(variable_name);
221
222 if (dimensions.empty()) {
223 // variable `variable_name` has no dimensions, so we check if it is a domain variable
224 auto dimension_list = file.read_text_attribute(variable_name, "dimensions");
225 dimensions = split(dimension_list, ' ');
226 }
227
228 if (dimensions.empty()) {
230 "variable '%s' in '%s' does not define a grid",
231 variable_name.c_str(), file.name().c_str());
232 }
233
234 std::string result = split(file.name(), '/').back();
235
236 for (const auto &d : dimensions) {
237 auto type = file.dimension_type(d, sys);
238
239 if (type == pism::X_AXIS or type == pism::Y_AXIS) {
240 result += ":";
241 result += d;
242 }
243 }
244
245 if (piecewise_constant) {
246 result += ":piecewise_constant";
247 }
248
249 return result;
250}
251
252/*!
253 * Return true if `string` contains "epsg:" or "EPSG:", otherwise false.
254 */
255static bool contains_epsg(const std::string &string) {
256 for (const auto &auth : { "epsg:", "EPSG:" }) {
257 if (string.find(auth) != std::string::npos) {
258 return true;
259 }
260 }
261 return false;
262}
263
264/*!
265 * Get PROJ parameters from the grid mapping variable `mapping_variable_name` in
266 * `input_file` *or* from global attributes.
267 *
268 * Return the string containing PROJ parameters.
269 *
270 * Tries
271 *
272 * - reading the "crs_wkt" attribute of the grid mapping variable (CF conventions)
273 * - reading the "proj_params" attribute of the grid mapping variable (convention used by CDO)
274 *
275 * If failed to discover PROJ parameters, tries global attributes "proj" and "proj4", in
276 * this order (PISM's old convention).
277 */
278static std::string get_proj_parameters(const File &input_file, const std::string &mapping_variable_name) {
279
280 std::string proj_string;
281
282 if (not mapping_variable_name.empty()) {
283 // try the crs_wkt attribute
284 proj_string = input_file.read_text_attribute(mapping_variable_name, "crs_wkt");
285
286 // if failed, try the "proj_params" attribute used by CDO:
287 if (proj_string.empty()) {
288 proj_string = input_file.read_text_attribute(mapping_variable_name, "proj_params");
289 }
290 }
291
292 // if failed, try the global attribute "proj"
293 if (proj_string.empty()) {
294 proj_string = input_file.read_text_attribute("PISM_GLOBAL", "proj");
295 }
296
297 // if failed, try the global attribute "proj4"
298 if (proj_string.empty()) {
299 proj_string = input_file.read_text_attribute("PISM_GLOBAL", "proj4");
300 }
301
302 if (contains_epsg(proj_string)) {
303 // re-create the PROJ string to make sure it does not contain "+init="
304 proj_string = pism::printf("EPSG:%d", parse_epsg(proj_string));
305 }
306
307 return proj_string;
308}
309
310
311/*!
312 * Get projection info for a variable `variable_name` from `input_file`.
313 *
314 * Obtains the string containing PROJ parameters as described in `get_proj_parameters()`.
315 * If the grid mapping variable has
316 */
317MappingInfo MappingInfo::FromFile(const File &input_file, const std::string &variable_name,
318 units::System::Ptr unit_system) {
319
320 auto mapping_variable_name = input_file.read_text_attribute(variable_name, "grid_mapping");
321
322 auto proj_string = get_proj_parameters(input_file, mapping_variable_name);
323
324 auto proj_is_epsg = contains_epsg(proj_string);
325
326 // Initialize (and possibly validate) the CF-style grid mapping variable by reading
327 // metadata from the input file:
328 VariableMetadata cf_mapping(mapping_variable_name, unit_system);
329 if (input_file.variable_exists(mapping_variable_name)) {
330 // input file has a mapping variable
331
332 cf_mapping = io::read_attributes(input_file, mapping_variable_name, unit_system);
333
334 // From the CF Conventions document, section 5.6:
335 //
336 // The one attribute that all grid mapping variables must have is grid_mapping_name, which
337 // takes a string value that contains the mapping’s name.
338
339 if (cf_mapping.has_attribute("grid_mapping_name") and proj_is_epsg) {
340 // The grid mapping variable contains the CF Conventions-style grid mapping
341 // definition: check consistency
342 try {
344 } catch (RuntimeError &e) {
345 e.add_context("getting projection info from variable '%s' in '%s'", variable_name.c_str(),
346 input_file.name().c_str());
347 throw;
348 }
349 }
350 }
351
352 if (cf_mapping.has_attribute("grid_mapping_name")) {
353 if (proj_string.empty()) {
354 // CF-style mapping was initialized, by the PROJ string was not: set proj_string
355 // from cf_mapping
357 }
358 } else {
359 if (proj_is_epsg) {
360 // cf_mapping was not initialized by the code above and proj_string contains an EPSG
361 // code: we may be able to convert to a CF-style mapping variable and set cf_mapping
362 // from proj_string
363 cf_mapping = epsg_to_cf(unit_system, proj_string);
364 }
365 }
366
367 return {cf_mapping, proj_string};
368}
369
371
372#if (Pism_USE_PROJ==1)
373
374//! Computes the area of a triangle using vector cross product.
375static double triangle_area(const double *A, const double *B, const double *C) {
376 double V1[3], V2[3];
377 for (int j = 0; j < 3; ++j) {
378 V1[j] = B[j] - A[j];
379 V2[j] = C[j] - A[j];
380 }
381 using std::pow;
382 using std::sqrt;
383 return 0.5*sqrt(pow(V1[1]*V2[2] - V2[1]*V1[2], 2) +
384 pow(V1[0]*V2[2] - V2[0]*V1[2], 2) +
385 pow(V1[0]*V2[1] - V2[0]*V1[1], 2));
386}
387
388void compute_cell_areas(const std::string &projection, array::Scalar &result) {
389 auto grid = result.grid();
390
391 Proj pism_to_geocent(projection, "+proj=geocent +datum=WGS84 +ellps=WGS84");
392
393// Cell layout:
394// (nw) (ne)
395// +-----------+
396// | |
397// | |
398// | o |
399// | |
400// | |
401// +-----------+
402// (sw) (se)
403
404 double dx2 = 0.5 * grid->dx(), dy2 = 0.5 * grid->dy();
405
406 array::AccessScope list(result);
407
408 for (auto p = grid->points(); p; p.next()) {
409 const int i = p.i(), j = p.j();
410
411 const double
412 x = grid->x(i),
413 y = grid->y(j);
414 double
415 x_nw = x - dx2, y_nw = y + dy2,
416 x_ne = x + dx2, y_ne = y + dy2,
417 x_se = x + dx2, y_se = y - dy2,
418 x_sw = x - dx2, y_sw = y - dy2;
419
420 PJ_COORD in, out;
421
422 in.xy = {x_nw, y_nw};
423 out = proj_trans(pism_to_geocent, PJ_FWD, in);
424 double nw[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
425
426 in.xy = {x_ne, y_ne};
427 out = proj_trans(pism_to_geocent, PJ_FWD, in);
428 double ne[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
429
430 in.xy = {x_se, y_se};
431 out = proj_trans(pism_to_geocent, PJ_FWD, in);
432 double se[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
433
434 in.xy = {x_sw, y_sw};
435 out = proj_trans(pism_to_geocent, PJ_FWD, in);
436 double sw[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
437
438 result(i, j) = triangle_area(sw, se, ne) + triangle_area(ne, nw, sw);
439 }
440}
441
442static void compute_lon_lat(const std::string &projection,
443 LonLat which, array::Scalar &result) {
444
445 Proj crs(projection, "EPSG:4326");
446
447// Cell layout:
448// (nw) (ne)
449// +-----------+
450// | |
451// | |
452// | o |
453// | |
454// | |
455// +-----------+
456// (sw) (se)
457
458 auto grid = result.grid();
459
460 array::AccessScope list{&result};
461
462 for (auto p = grid->points(); p; p.next()) {
463 const int i = p.i(), j = p.j();
464
465 PJ_COORD in, out;
466
467 in.xy = {grid->x(i), grid->y(j)};
468 out = proj_trans(crs, PJ_FWD, in);
469
470 if (which == LONGITUDE) {
471 result(i, j) = out.lp.phi;
472 } else {
473 result(i, j) = out.lp.lam;
474 }
475 }
476}
477
478static void compute_lon_lat_bounds(const std::string &projection,
479 LonLat which,
480 array::Array3D &result) {
481
482 Proj crs(projection, "EPSG:4326");
483
484 auto grid = result.grid();
485
486 double dx2 = 0.5 * grid->dx(), dy2 = 0.5 * grid->dy();
487 double x_offsets[] = {-dx2, dx2, dx2, -dx2};
488 double y_offsets[] = {-dy2, -dy2, dy2, dy2};
489
490 array::AccessScope list{&result};
491
492 for (auto p = grid->points(); p; p.next()) {
493 const int i = p.i(), j = p.j();
494
495 double x0 = grid->x(i), y0 = grid->y(j);
496
497 double *values = result.get_column(i, j);
498
499 for (int k = 0; k < 4; ++k) {
500
501 PJ_COORD in, out;
502
503 in.xy = {x0 + x_offsets[k], y0 + y_offsets[k]};
504
505 // compute lon,lat coordinates:
506 out = proj_trans(crs, PJ_FWD, in);
507
508 if (which == LATITUDE) {
509 values[k] = out.lp.lam;
510 } else {
511 values[k] = out.lp.phi;
512 }
513 }
514 }
515}
516
517#else
518
519void compute_cell_areas(const std::string &projection, array::Scalar &result) {
520 (void) projection;
521
522 auto grid = result.grid();
523 result.set(grid->dx() * grid->dy());
524}
525
526static void compute_lon_lat(const std::string &projection, LonLat which,
527 array::Scalar &result) {
528 (void) projection;
529 (void) which;
530 (void) result;
531
532 throw RuntimeError(PISM_ERROR_LOCATION, "Cannot compute longitude and latitude."
533 " Please rebuild PISM with PROJ.");
534}
535
536static void compute_lon_lat_bounds(const std::string &projection,
537 LonLat which,
538 array::Array3D &result) {
539 (void) projection;
540 (void) which;
541 (void) result;
542
543 throw RuntimeError(PISM_ERROR_LOCATION, "Cannot compute longitude and latitude bounds."
544 " Please rebuild PISM with PROJ.");
545}
546
547#endif
548
549void compute_longitude(const std::string &projection, array::Scalar &result) {
550 compute_lon_lat(projection, LONGITUDE, result);
551}
552void compute_latitude(const std::string &projection, array::Scalar &result) {
553 compute_lon_lat(projection, LATITUDE, result);
554}
555
556void compute_lon_bounds(const std::string &projection, array::Array3D &result) {
557 compute_lon_lat_bounds(projection, LONGITUDE, result);
558}
559
560void compute_lat_bounds(const std::string &projection, array::Array3D &result) {
561 compute_lon_lat_bounds(projection, LATITUDE, result);
562}
563
564static std::string albers_conical_equal_area_to_proj(const VariableMetadata &mapping) {
565 // standard_parallel - There may be 1 or 2 values.
566 // longitude_of_central_meridian
567 // latitude_of_projection_origin
568 // false_easting
569 // false_northing
570 double longitude_of_projection_origin = mapping["longitude_of_projection_origin"];
571 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
572 double false_easting = mapping["false_easting"];
573 double false_northing = mapping["false_northing"];
574
575 std::vector<double> standard_parallel = mapping["standard_parallel"];
576
577 auto result = pism::printf("+proj=aea +type=crs +lon_0=%f +lat_0=%f +x_0=%f +y_0=%f +lat_1=%f",
578 longitude_of_projection_origin,
579 latitude_of_projection_origin,
580 false_easting,
581 false_northing,
582 standard_parallel[0]);
583 if (standard_parallel.size() == 2) {
584 result += pism::printf(" +lat_2=%f", standard_parallel[1]);
585 }
586
587 return result;
588}
589
590static std::string azimuthal_equidistant_to_proj(const VariableMetadata &mapping) {
591 // Parameters:
592 // longitude_of_projection_origin
593 // latitude_of_projection_origin
594 // false_easting
595 // false_northing
596 double longitude_of_projection_origin = mapping["longitude_of_projection_origin"];
597 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
598 double false_easting = mapping["false_easting"];
599 double false_northing = mapping["false_northing"];
600 return pism::printf("+proj=aeqd +type=crs +lon_0=%f +lat_0=%f +x_0=%f +y_0=%f",
601 longitude_of_projection_origin,
602 latitude_of_projection_origin, false_easting, false_northing);
603}
604
605static std::string lambert_azimuthal_equal_area_to_proj(const VariableMetadata &mapping) {
606 // longitude_of_projection_origin
607 // latitude_of_projection_origin
608 // false_easting
609 // false_northing
610 double longitude_of_projection_origin = mapping["longitude_of_projection_origin"];
611 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
612 double false_easting = mapping["false_easting"];
613 double false_northing = mapping["false_northing"];
614 return pism::printf("+proj=laea +type=crs +lon_0=%f +lat_0=%f +x_0=%f +y_0=%f",
615 longitude_of_projection_origin,
616 latitude_of_projection_origin, false_easting, false_northing);
617}
618
619static std::string lambert_conformal_conic_to_proj(const VariableMetadata &mapping) {
620 // standard_parallel - There may be 1 or 2 values.
621 // longitude_of_central_meridian
622 // latitude_of_projection_origin
623 // false_easting
624 // false_northing
625 double longitude_of_projection_origin = mapping["longitude_of_projection_origin"];
626 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
627 double false_easting = mapping["false_easting"];
628 double false_northing = mapping["false_northing"];
629 std::vector<double> standard_parallel = mapping["standard_parallel"];
630
631 auto result = pism::printf("+proj=lcc +type=crs +lon_0=%f +lat_0=%f +x_0=%f +y_0=%f +lat_1=%f",
632 longitude_of_projection_origin,
633 latitude_of_projection_origin,
634 false_easting,
635 false_northing,
636 standard_parallel[0]);
637 if (standard_parallel.size() == 2) {
638 result += pism::printf(" +lat_2=%f", standard_parallel[1]);
639 }
640 return result;
641}
642
644 // Parameters:
645 // longitude_of_central_meridian
646 // Either standard_parallel or scale_factor_at_projection_origin
647 // false_easting
648 // false_northing
649
650 double longitude_of_central_meridian = mapping["longitude_of_central_meridian"];
651 double false_easting = mapping["false_easting"];
652 double false_northing = mapping["false_northing"];
653
654 auto result = pism::printf("+proj=cea +type=crs +lon_0=%f +x_0=%f +y_0=%f",
655 longitude_of_central_meridian, false_easting, false_northing);
656
657 if (mapping.has_attribute("standard_parallel")) {
658 result += pism::printf(" +lat_ts=%f", (double)mapping["standard_parallel"]);
659 } else {
660 result += pism::printf(" +k_0=%f", (double)mapping["scale_factor_at_projection_origin"]);
661 }
662
663 return result;
664}
665
666static std::string latitude_longitude_to_proj(const VariableMetadata &/* unused */) {
667 // No parameters
668 return "+proj=lonlat +type=crs";
669}
670
671static std::string mercator_to_proj(const VariableMetadata &mapping) {
672 // Parameters:
673 // longitude_of_projection_origin
674 // Either standard_parallel (EPSG 9805) or scale_factor_at_projection_origin (EPSG 9804)
675 // false_easting
676 // false_northing
677
678 double longitude_of_projection_origin = mapping["longitude_of_projection_origin"];
679 double false_easting = mapping["false_easting"];
680 double false_northing = mapping["false_northing"];
681
682 auto result = pism::printf("+proj=merc +type=crs +lon_0=%f +x_0=%f +y_0=%f",
683 longitude_of_projection_origin,
684 false_easting, false_northing);
685
686 if (mapping.has_attribute("standard_parallel")) {
687 result += pism::printf(" +lat_ts=%f", (double)mapping["standard_parallel"]);
688 } else {
689 result += pism::printf(" +k_0=%f", (double)mapping["scale_factor_at_projection_origin"]);
690 }
691
692 return result;
693}
694
695static std::string orthographic_to_proj(const VariableMetadata &mapping) {
696 // Parameters:
697 // longitude_of_projection_origin
698 // latitude_of_projection_origin
699 // false_easting
700 // false_northing
701 double longitude_of_projection_origin = mapping["longitude_of_projection_origin"];
702 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
703 double false_easting = mapping["false_easting"];
704 double false_northing = mapping["false_northing"];
705
706 // clang-format off
707 return pism::printf("+proj=ortho +type=crs +lon_0=%f +lat_0=%f +x_0=%f +y_0=%f",
708 longitude_of_projection_origin,
709 latitude_of_projection_origin,
710 false_easting,
711 false_northing);
712 // clang-format on
713}
714
715static std::string polar_stereographic_to_proj(const VariableMetadata &mapping) {
716 /* Parameters:
717 straight_vertical_longitude_from_pole
718 latitude_of_projection_origin - Either +90. or -90.
719 Either standard_parallel or scale_factor_at_projection_origin
720 false_easting
721 false_northing
722 */
723
724 double straight_vertical_longitude_from_pole = mapping["straight_vertical_longitude_from_pole"];
725 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
726 double false_easting = mapping["false_easting"];
727 double false_northing = mapping["false_northing"];
728
729 auto result = pism::printf("+proj=stere +type=crs +lat_0=%f +lon_0=%f +x_0=%f +y_0=%f",
730 latitude_of_projection_origin,
731 straight_vertical_longitude_from_pole,
732 false_easting,
733 false_northing);
734
735 if (mapping.has_attribute("standard_parallel")) {
736 result += pism::printf(" +lat_ts=%f", (double)mapping["standard_parallel"]);
737 } else {
738 result += pism::printf(" +k_0=%f", (double)mapping["scale_factor_at_projection_origin"]);
739 }
740
741 return result;
742}
743
744static std::string rotated_latitude_longitude_to_proj(const VariableMetadata &mapping) {
745 // grid_north_pole_latitude
746 // grid_north_pole_longitude
747 // north_pole_grid_longitude - This parameter is optional (default is 0).
748
749 /*
750 * From https://github.com/OSGeo/PROJ/blob/356e255022cea47bf242bc9e6345a1e87ab1031d/src/iso19111/operation/conversion.cpp#L2419
751 *
752 * Several conventions for the pole rotation method exists.
753 * The parameters provided in this method are remapped to the PROJ ob_tran
754 * operation with:
755 * <pre>
756 * +proj=ob_tran +o_proj=longlat +o_lon_p=northPoleGridLongitude
757 * +o_lat_p=gridNorthPoleLatitude
758 * +lon_0=180+gridNorthPoleLongitude
759 * </pre>
760 */
761
762 double grid_north_pole_latitude = mapping["grid_north_pole_latitude"];
763 double grid_north_pole_longitude = mapping["grid_north_pole_longitude"];
764 double north_pole_grid_longitude = 0.0;
765
766 if (mapping.has_attribute("north_pole_grid_longitude")) {
767 north_pole_grid_longitude = mapping["north_pole_grid_longitude"];
768 }
769
770 return pism::printf("+proj=ob_tran +type=crs +o_proj=longlat +o_lon_p=%f +o_lat_p=%f +lon_0=%f",
771 north_pole_grid_longitude,
772 grid_north_pole_latitude,
773 180.0 + grid_north_pole_longitude);
774}
775
776static std::string stereographic_to_proj(const VariableMetadata &mapping) {
777 /* Parameters:
778 longitude_of_projection_origin
779 latitude_of_projection_origin
780 scale_factor_at_projection_origin
781 false_easting
782 false_northing
783*/
784
785 double straight_vertical_longitude_from_pole = mapping["straight_vertical_longitude_from_pole"];
786 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
787 double false_easting = mapping["false_easting"];
788 double false_northing = mapping["false_northing"];
789 double scale_factor_at_projection_origin = mapping["scale_factor_at_projection_origin"];
790
791 // clang-format off
792 return pism::printf("+proj=stere +type=crs +lat_0=%f +lon_0=%f +x_0=%f +y_0=%f +k_0=%f",
793 latitude_of_projection_origin,
794 straight_vertical_longitude_from_pole,
795 false_easting,
796 false_northing,
797 scale_factor_at_projection_origin);
798 // clang-format on
799}
800
801static std::string transverse_mercator_to_proj(const VariableMetadata &mapping) {
802 double scale_factor_at_central_meridian = mapping["scale_factor_at_central_meridian"];
803 double longitude_of_central_meridian = mapping["longitude_of_central_meridian"];
804 double latitude_of_projection_origin = mapping["latitude_of_projection_origin"];
805 double false_easting = mapping["false_easting"];
806 double false_northing = mapping["false_northing"];
807
808 return pism::printf("+proj=tmerc +type=crs +lon_0=%f +lat_0=%f +k=%f +x_0=%f +y_0=%f",
809 longitude_of_central_meridian,
810 latitude_of_projection_origin,
811 scale_factor_at_central_meridian,
812 false_easting,
813 false_northing);
814}
815
816static std::string vertical_perspective_to_proj(const VariableMetadata &mapping) {
817 // Parameters:
818 // latitude_of_projection_origin
819 // longitude_of_projection_origin
820 // perspective_point_height
821 // false_easting
822 // false_northing
823 std::string name = mapping["grid_mapping_name"];
824 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "grid mapping '%s' is not supported",
825 name.c_str());
826 return "";
827}
828
829static std::string common_flags(const VariableMetadata &mapping) {
830
831 auto flag = [&mapping](const char* name, const char* proj_name) {
832 if (mapping.has_attribute(name)) {
833 return pism::printf(" +%s=%f", proj_name, (double)mapping[name]);
834 }
835 return std::string{};
836 };
837
838 std::string result;
839
840 // ellipsoid
841 result += flag("earth_radius", "R");
842 result += flag("inverse_flattening", "rf");
843 result += flag("semi_major_axis", "a");
844 result += flag("semi_minor_axis", "b");
845
846 // prime meridian
847 result += flag("longitude_of_prime_meridian", "pm");
848
849 return result;
850}
851
852/*!
853 * Convert a CF-style grid mapping definition to a PROJ string.
854 *
855 * See http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/apf.html
856 *
857 * See also https://proj.org/en/9.4/operations/projections/index.html for a list of
858 * projections and their parameters.
859 *
860 * PROJ strings returned by this function are inspired by
861 * https://scitools.org.uk/cartopy/docs/latest/index.html
862 */
863std::string cf_to_proj(const VariableMetadata &mapping) {
864
865 // all projections: use the crs_wkt attribute if it is present
866 std::string wkt = mapping["crs_wkt"];
867 if (not wkt.empty()) {
868 return wkt;
869 }
870
871 typedef std::string(*cf_to_proj_fn)(const VariableMetadata&);
872
873 std::map<std::string, cf_to_proj_fn> mappings = {
874 { "albers_conical_equal_area", albers_conical_equal_area_to_proj },
875 { "azimuthal_equidistant", azimuthal_equidistant_to_proj },
876 { "lambert_azimuthal_equal_area", lambert_azimuthal_equal_area_to_proj },
877 { "lambert_conformal_conic", lambert_conformal_conic_to_proj },
878 { "lambert_cylindrical_equal_area", lambert_cylindrical_equal_area_to_proj },
879 { "latitude_longitude", latitude_longitude_to_proj },
880 { "mercator", mercator_to_proj },
881 { "orthographic", orthographic_to_proj },
882 { "polar_stereographic", polar_stereographic_to_proj },
883 { "rotated_latitude_longitude", rotated_latitude_longitude_to_proj },
884 { "stereographic", stereographic_to_proj },
885 { "transverse_mercator", transverse_mercator_to_proj },
886 { "vertical_perspective", vertical_perspective_to_proj },
887 };
888
889 std::string mapping_name = mapping["grid_mapping_name"];
890
891 if (mappings.find(mapping_name) == mappings.end()) {
893 "conversion CF -> PROJ for a '%s' grid is not implemented",
894 mapping_name.c_str());
895 }
896
897 return mappings[mapping_name](mapping) + common_flags(mapping);
898}
899
900void write_mapping(const File &file, const pism::MappingInfo &info) {
901
902 const auto &mapping = info.cf_mapping;
903 if (mapping.has_attributes()) {
904 auto name = mapping.get_name();
905 if (not file.variable_exists(name)) {
906 file.define_variable(name, io::PISM_DOUBLE, {});
907 }
909
910 // Write the PROJ string to mapping:proj_params (for CDO).
911 if (not info.proj_string.empty()) {
912 file.write_attribute(name, "proj_params", info.proj_string);
913 }
914 }
915}
916
917} // end of namespace pism
AxisType dimension_type(const std::string &name, units::System::Ptr unit_system) const
Get the "type" of a dimension.
Definition File.cc:460
bool variable_exists(const std::string &short_name) const
Checks if a variable exists.
Definition File.cc:378
void define_variable(const std::string &name, io::Type nctype, const std::vector< std::string > &dims) const
Define a variable.
Definition File.cc:543
std::string name() const
Definition File.cc:274
void write_attribute(const std::string &var_name, const std::string &att_name, io::Type nctype, const std::vector< double > &values) const
Write a multiple-valued double attribute.
Definition File.cc:590
std::vector< std::string > dimensions(const std::string &variable_name) const
Definition File.cc:390
std::string read_text_attribute(const std::string &var_name, const std::string &att_name) const
Get a text attribute.
Definition File.cc:645
High-level PISM I/O class.
Definition File.hh:55
MappingInfo(const std::string &mapping_variable_name, units::System::Ptr unit_system)
Definition projection.cc:44
VariableMetadata cf_mapping
grid mapping description following CF conventions
Definition projection.hh:69
std::string proj_string
a projection definition string in a format recognized by PROJ 6.x+
Definition projection.hh:72
static MappingInfo FromFile(const File &input_file, const std::string &variable_name, units::System::Ptr unit_system)
Get projection info from a file.
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
double get_number(const std::string &name) const
Get a single-valued scalar attribute.
const std::map< std::string, std::string > & all_strings() const
units::System::Ptr unit_system() const
bool has_attribute(const std::string &name) const
std::string get_name() const
const std::map< std::string, std::vector< double > > & all_doubles() const
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition Array3D.hh:33
std::shared_ptr< const Grid > grid() const
Definition Array.cc:131
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
std::shared_ptr< System > Ptr
Definition Units.hh:47
#define PISM_ERROR_LOCATION
VariableMetadata read_attributes(const File &file, const std::string &variable_name, std::shared_ptr< units::System > unit_system)
void write_attributes(const File &file, const VariableMetadata &variable, io::Type nctype)
Write variable attributes to a NetCDF file.
@ PISM_DOUBLE
Definition IO_Flags.hh:52
static std::string lambert_azimuthal_equal_area_to_proj(const VariableMetadata &mapping)
void check_consistency_epsg(const VariableMetadata &cf_mapping, const std::string &proj_string)
Check consistency of the "mapping" variable with the EPSG code in the proj string.
static bool contains_epsg(const std::string &string)
static std::string get_proj_parameters(const File &input_file, const std::string &mapping_variable_name)
@ X_AXIS
Definition IO_Flags.hh:33
@ Y_AXIS
Definition IO_Flags.hh:33
void write_mapping(const File &file, const pism::MappingInfo &info)
static std::string mercator_to_proj(const VariableMetadata &mapping)
static std::string stereographic_to_proj(const VariableMetadata &mapping)
static std::string polar_stereographic_to_proj(const VariableMetadata &mapping)
static void compute_lon_lat_bounds(const std::string &projection, LonLat which, array::Array3D &result)
void compute_latitude(const std::string &projection, array::Scalar &result)
int parse_epsg(const std::string &proj_string)
Definition projection.cc:54
static std::string lambert_conformal_conic_to_proj(const VariableMetadata &mapping)
std::string printf(const char *format,...)
void compute_longitude(const std::string &projection, array::Scalar &result)
static std::string azimuthal_equidistant_to_proj(const VariableMetadata &mapping)
VariableMetadata epsg_to_cf(units::System::Ptr system, const std::string &proj_string)
Return CF-Convention "mapping" variable corresponding to an EPSG code specified in a PROJ string.
Definition projection.cc:90
static const double k
Definition exactTestP.cc:42
std::string grid_name(const pism::File &file, const std::string &variable_name, pism::units::System::Ptr sys, bool piecewise_constant)
static std::string vertical_perspective_to_proj(const VariableMetadata &mapping)
@ LONGITUDE
@ LATITUDE
static std::string transverse_mercator_to_proj(const VariableMetadata &mapping)
static double triangle_area(const Point &a, const Point &b, const Point &c)
static void compute_lon_lat(const std::string &projection, LonLat which, array::Scalar &result)
static std::string albers_conical_equal_area_to_proj(const VariableMetadata &mapping)
void compute_lat_bounds(const std::string &projection, array::Array3D &result)
static std::string rotated_latitude_longitude_to_proj(const VariableMetadata &mapping)
static std::string lambert_cylindrical_equal_area_to_proj(const VariableMetadata &mapping)
void compute_cell_areas(const std::string &projection, array::Scalar &result)
void compute_lon_bounds(const std::string &projection, array::Array3D &result)
std::string cf_to_proj(const VariableMetadata &mapping)
static std::string orthographic_to_proj(const VariableMetadata &mapping)
static std::string latitude_longitude_to_proj(const VariableMetadata &)
static std::string common_flags(const VariableMetadata &mapping)
std::vector< std::string > split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a vector of strings.
std::shared_ptr< Grid > grid(std::shared_ptr< Context > ctx, char testname)
Definition pism.cc:176