tprojection.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tprojection.cc (14160B)
---
1 /* Copyright (C) 2015, 2016, 2017, 2018, 2019, 2020 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> // fabs
22
23 #include "projection.hh"
24 #include "VariableMetadata.hh"
25 #include "error_handling.hh"
26 #include "io/File.hh"
27 #include "io/io_helpers.hh"
28 #include "pism/util/IceGrid.hh"
29 #include "pism/util/iceModelVec.hh"
30
31 #include "pism/pism_config.hh"
32
33 #if (Pism_USE_PROJ==1)
34 #include "pism/util/Proj.hh"
35 #endif
36
37 namespace pism {
38
39 MappingInfo::MappingInfo(const std::string &mapping_name, units::System::Ptr unit_system)
40 : mapping(mapping_name, unit_system) {
41 // empty
42 }
43
44 //! @brief Return CF-Convention "mapping" variable corresponding to an EPSG code specified in a
45 //! PROJ string.
46 VariableMetadata epsg_to_cf(units::System::Ptr system, const std::string &proj_string) {
47 VariableMetadata mapping("mapping", system);
48
49 int auth_len = 5; // length of "epsg:"
50 std::string::size_type position = std::string::npos;
51 for (auto &auth : {"epsg:", "EPSG:"}) {
52 position = proj_string.find(auth);
53 if (position != std::string::npos) {
54 break;
55 }
56 }
57
58 if (position == std::string::npos) {
59 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
60 "could not parse EPSG code '%s'", proj_string.c_str());
61 }
62
63 long int epsg = 0;
64 {
65 std::string epsg_code = proj_string.substr(position + auth_len);
66 const char *str = epsg_code.c_str();
67 char *endptr = NULL;
68 epsg = strtol(str, &endptr, 10);
69 if (endptr == str) {
70 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
71 "failed to parse EPSG code at '%s' in PROJ string '%s'",
72 epsg_code.c_str(), proj_string.c_str());
73 }
74 }
75
76 switch (epsg) {
77 case 3413:
78 mapping.set_number("latitude_of_projection_origin", 90.0);
79 mapping.set_number("scale_factor_at_projection_origin", 1.0);
80 mapping.set_number("straight_vertical_longitude_from_pole", -45.0);
81 mapping.set_number("standard_parallel", 70.0);
82 mapping.set_number("false_northing", 0.0);
83 mapping.set_string("grid_mapping_name", "polar_stereographic");
84 mapping.set_number("false_easting", 0.0);
85 break;
86 case 3031:
87 mapping.set_number("latitude_of_projection_origin", -90.0);
88 mapping.set_number("scale_factor_at_projection_origin", 1.0);
89 mapping.set_number("straight_vertical_longitude_from_pole", 0.0);
90 mapping.set_number("standard_parallel", -71.0);
91 mapping.set_number("false_northing", 0.0);
92 mapping.set_string("grid_mapping_name", "polar_stereographic");
93 mapping.set_number("false_easting", 0.0);
94 break;
95 case 3057:
96 mapping.set_string("grid_mapping_name", "lambert_conformal_conic") ;
97 mapping.set_number("longitude_of_central_meridian", -19.) ;
98 mapping.set_number("false_easting", 500000.) ;
99 mapping.set_number("false_northing", 500000.) ;
100 mapping.set_number("latitude_of_projection_origin", 65.) ;
101 mapping.set_numbers("standard_parallel", {64.25, 65.75}) ;
102 mapping.set_string("long_name", "CRS definition") ;
103 mapping.set_number("longitude_of_prime_meridian", 0.) ;
104 mapping.set_number("semi_major_axis", 6378137.) ;
105 mapping.set_number("inverse_flattening", 298.257222101) ;
106 break;
107 case 5936:
108 mapping.set_number("latitude_of_projection_origin", 90.0);
109 mapping.set_number("scale_factor_at_projection_origin", 1.0);
110 mapping.set_number("straight_vertical_longitude_from_pole", -150.0);
111 mapping.set_number("standard_parallel", 90.0);
112 mapping.set_number("false_northing", 2000000.0);
113 mapping.set_string("grid_mapping_name", "polar_stereographic");
114 mapping.set_number("false_easting", 2000000.0);
115 break;
116 case 26710:
117 mapping.set_number("longitude_of_central_meridian", -123.0);
118 mapping.set_number("false_easting", 500000.0);
119 mapping.set_number("false_northing", 0.0);
120 mapping.set_string("grid_mapping_name", "transverse_mercator");
121 mapping.set_number("inverse_flattening", 294.978698213898);
122 mapping.set_number("latitude_of_projection_origin", 0.0);
123 mapping.set_number("scale_factor_at_central_meridian", 0.9996);
124 mapping.set_number("semi_major_axis", 6378206.4);
125 mapping.set_string("unit", "metre");
126 break;
127 default:
128 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "unknown EPSG code '%d' in PROJ string '%s'",
129 (int)epsg, proj_string.c_str());
130 }
131
132 return mapping;
133 }
134
135 void check_consistency_epsg(const MappingInfo &info) {
136
137 VariableMetadata epsg_mapping = epsg_to_cf(info.mapping.unit_system(), info.proj);
138
139 bool mapping_is_empty = not info.mapping.has_attributes();
140 bool epsg_mapping_is_empty = not epsg_mapping.has_attributes();
141
142 if (mapping_is_empty and epsg_mapping_is_empty) {
143 // empty mapping variables are equivalent
144 return;
145 } else {
146 // Check if the "info.mapping" variable in the input file matches the EPSG code.
147 // Check strings.
148 for (auto s : epsg_mapping.get_all_strings()) {
149 if (not info.mapping.has_attribute(s.first)) {
150 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "inconsistent metadata:\n"
151 "PROJ string \"%s\" requires %s = \"%s\",\n"
152 "but the mapping variable has no %s.",
153 info.proj.c_str(),
154 s.first.c_str(), s.second.c_str(),
155 s.first.c_str());
156 }
157
158 std::string string = info.mapping.get_string(s.first);
159
160 if (not (string == s.second)) {
161 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "inconsistent metadata:\n"
162 "%s requires %s = \"%s\",\n"
163 "but the mapping variable has %s = \"%s\".",
164 info.proj.c_str(),
165 s.first.c_str(), s.second.c_str(),
166 s.first.c_str(),
167 string.c_str());
168 }
169 }
170
171 // Check doubles
172 for (auto d : epsg_mapping.get_all_doubles()) {
173 if (not info.mapping.has_attribute(d.first)) {
174 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "inconsistent metadata:\n"
175 "%s requires %s = %f,\n"
176 "but the mapping variable has no %s.",
177 info.proj.c_str(),
178 d.first.c_str(), d.second[0],
179 d.first.c_str());
180 }
181
182 double value = info.mapping.get_number(d.first);
183
184 if (fabs(value - d.second[0]) > 1e-12) {
185 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "inconsistent metadata:\n"
186 "%s requires %s = %f,\n"
187 "but the mapping variable has %s = %f.",
188 info.proj.c_str(),
189 d.first.c_str(), d.second[0],
190 d.first.c_str(),
191 value);
192 }
193 }
194 }
195 }
196
197 MappingInfo get_projection_info(const File &input_file, const std::string &mapping_name,
198 units::System::Ptr unit_system) {
199 MappingInfo result(mapping_name, unit_system);
200
201 result.proj = input_file.read_text_attribute("PISM_GLOBAL", "proj");
202
203 bool proj_is_epsg = false;
204 for (auto &auth : {"epsg:", "EPSG:"}) {
205 if (result.proj.find(auth) != std::string::npos) {
206 proj_is_epsg = true;
207 break;
208 }
209 }
210
211 if (input_file.find_variable(mapping_name)) {
212 // input file has a mapping variable
213
214 io::read_attributes(input_file, mapping_name, result.mapping);
215
216 if (proj_is_epsg) {
217 // check consistency
218 try {
219 check_consistency_epsg(result);
220 } catch (RuntimeError &e) {
221 e.add_context("getting projection info from %s", input_file.filename().c_str());
222 throw;
223 }
224 } else {
225 // use mapping read from input_file (can't check consistency here)
226 }
227 } else {
228 // no mapping variable in the input file
229
230 if (proj_is_epsg) {
231 result.mapping = epsg_to_cf(unit_system, result.proj);
232 } else {
233 // leave mapping empty
234 }
235 }
236 return result;
237 }
238
239 enum LonLat {LONGITUDE, LATITUDE};
240
241 #if (Pism_USE_PROJ==1)
242
243 //! Computes the area of a triangle using vector cross product.
244 static double triangle_area(double *A, double *B, double *C) {
245 double V1[3], V2[3];
246 for (int j = 0; j < 3; ++j) {
247 V1[j] = B[j] - A[j];
248 V2[j] = C[j] - A[j];
249 }
250
251 return 0.5*sqrt(PetscSqr(V1[1]*V2[2] - V2[1]*V1[2]) +
252 PetscSqr(V1[0]*V2[2] - V2[0]*V1[2]) +
253 PetscSqr(V1[0]*V2[1] - V2[0]*V1[1]));
254 }
255
256 void compute_cell_areas(const std::string &projection, IceModelVec2S &result) {
257 IceGrid::ConstPtr grid = result.grid();
258
259 Proj pism_to_geocent(projection, "+proj=geocent +datum=WGS84 +ellps=WGS84");
260
261 // Cell layout:
262 // (nw) (ne)
263 // +-----------+
264 // | |
265 // | |
266 // | o |
267 // | |
268 // | |
269 // +-----------+
270 // (sw) (se)
271
272 double dx2 = 0.5 * grid->dx(), dy2 = 0.5 * grid->dy();
273
274 IceModelVec::AccessList list(result);
275
276 for (Points p(*grid); p; p.next()) {
277 const int i = p.i(), j = p.j();
278
279 const double
280 x = grid->x(i),
281 y = grid->y(j);
282 double
283 x_nw = x - dx2, y_nw = y + dy2,
284 x_ne = x + dx2, y_ne = y + dy2,
285 x_se = x + dx2, y_se = y - dy2,
286 x_sw = x - dx2, y_sw = y - dy2;
287
288 PJ_COORD in, out;
289
290 in.xy = {x_nw, y_nw};
291 out = proj_trans(*pism_to_geocent, PJ_FWD, in);
292 double nw[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
293
294 in.xy = {x_ne, y_ne};
295 out = proj_trans(*pism_to_geocent, PJ_FWD, in);
296 double ne[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
297
298 in.xy = {x_se, y_se};
299 out = proj_trans(*pism_to_geocent, PJ_FWD, in);
300 double se[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
301
302 in.xy = {x_sw, y_sw};
303 out = proj_trans(*pism_to_geocent, PJ_FWD, in);
304 double sw[3] = {out.xyz.x, out.xyz.y, out.xyz.z};
305
306 result(i, j) = triangle_area(sw, se, ne) + triangle_area(ne, nw, sw);
307 }
308 }
309
310 static void compute_lon_lat(const std::string &projection,
311 LonLat which, IceModelVec2S &result) {
312
313 Proj crs(projection, "EPSG:4326");
314
315 // Cell layout:
316 // (nw) (ne)
317 // +-----------+
318 // | |
319 // | |
320 // | o |
321 // | |
322 // | |
323 // +-----------+
324 // (sw) (se)
325
326 IceGrid::ConstPtr grid = result.grid();
327
328 IceModelVec::AccessList list{&result};
329
330 for (Points p(*grid); p; p.next()) {
331 const int i = p.i(), j = p.j();
332
333 PJ_COORD in, out;
334
335 in.xy = {grid->x(i), grid->y(j)};
336 out = proj_trans(*crs, PJ_FWD, in);
337
338 if (which == LONGITUDE) {
339 result(i, j) = out.lp.phi;
340 } else {
341 result(i, j) = out.lp.lam;
342 }
343 }
344 }
345
346 static void compute_lon_lat_bounds(const std::string &projection,
347 LonLat which,
348 IceModelVec3D &result) {
349
350 Proj crs(projection, "EPSG:4326");
351
352 IceGrid::ConstPtr grid = result.grid();
353
354 double dx2 = 0.5 * grid->dx(), dy2 = 0.5 * grid->dy();
355 double x_offsets[] = {-dx2, dx2, dx2, -dx2};
356 double y_offsets[] = {-dy2, -dy2, dy2, dy2};
357
358 IceModelVec::AccessList list{&result};
359
360 for (Points p(*grid); p; p.next()) {
361 const int i = p.i(), j = p.j();
362
363 double x0 = grid->x(i), y0 = grid->y(j);
364
365 double *values = result.get_column(i, j);
366
367 for (int k = 0; k < 4; ++k) {
368
369 PJ_COORD in, out;
370
371 in.xy = {x0 + x_offsets[k], y0 + y_offsets[k]};
372
373 // compute lon,lat coordinates:
374 out = proj_trans(*crs, PJ_FWD, in);
375
376 if (which == LATITUDE) {
377 values[k] = out.lp.lam;
378 } else {
379 values[k] = out.lp.phi;
380 }
381 }
382 }
383 }
384
385 #else
386
387 void compute_cell_areas(const std::string &projection, IceModelVec2S &result) {
388 (void) projection;
389
390 IceGrid::ConstPtr grid = result.grid();
391 result.set(grid->dx() * grid->dy());
392 }
393
394 static void compute_lon_lat(const std::string &projection, LonLat which,
395 IceModelVec2S &result) {
396 (void) projection;
397 (void) which;
398 (void) result;
399
400 throw RuntimeError(PISM_ERROR_LOCATION, "Cannot compute longitude and latitude."
401 " Please rebuild PISM with PROJ.");
402 }
403
404 static void compute_lon_lat_bounds(const std::string &projection,
405 LonLat which,
406 IceModelVec3D &result) {
407 (void) projection;
408 (void) which;
409 (void) result;
410
411 throw RuntimeError(PISM_ERROR_LOCATION, "Cannot compute longitude and latitude bounds."
412 " Please rebuild PISM with PROJ.");
413 }
414
415 #endif
416
417 void compute_longitude(const std::string &projection, IceModelVec2S &result) {
418 compute_lon_lat(projection, LONGITUDE, result);
419 }
420 void compute_latitude(const std::string &projection, IceModelVec2S &result) {
421 compute_lon_lat(projection, LATITUDE, result);
422 }
423
424 void compute_lon_bounds(const std::string &projection, IceModelVec3D &result) {
425 compute_lon_lat_bounds(projection, LONGITUDE, result);
426 }
427
428 void compute_lat_bounds(const std::string &projection, IceModelVec3D &result) {
429 compute_lon_lat_bounds(projection, LATITUDE, result);
430 }
431
432 } // end of namespace pism