projection.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