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"
35#include "pism/pism_config.hh"
36#include "pism/util/pism_utilities.hh"
39#include "pism/util/Proj.hh"
45 : cf_mapping(mapping_variable_name, unit_system) {
50 : cf_mapping(mapping_variable), proj_string(proj_string_) {
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) {
65 if (position == std::string::npos) {
71 auto epsg_string = proj_string.substr(position + auth_len);
73 long int result = strtol(epsg_string.c_str(), &endptr, 10);
74 if (endptr == epsg_string.c_str()) {
82 e.
add_context(
"parsing the EPSG code in the PROJ string '%s'",
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};
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};
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} ;
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};
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";
148 (
int)epsg, proj_string.c_str());
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");
162 if (mapping_is_empty and epsg_mapping_is_empty) {
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(),
180 std::string
string = cf_mapping[s.first];
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());
193 const double eps = 1e-12;
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],
204 double value = cf_mapping.
get_number(d.first);
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);
220 std::vector<std::string> dimensions = file.
dimensions(variable_name);
222 if (dimensions.empty()) {
225 dimensions =
split(dimension_list,
' ');
228 if (dimensions.empty()) {
230 "variable '%s' in '%s' does not define a grid",
231 variable_name.c_str(), file.
name().c_str());
234 std::string result =
split(file.
name(),
'/').back();
236 for (
const auto &d : dimensions) {
245 if (piecewise_constant) {
246 result +=
":piecewise_constant";
256 for (
const auto &auth : {
"epsg:",
"EPSG:" }) {
257 if (
string.find(auth) != std::string::npos) {
280 std::string proj_string;
282 if (not mapping_variable_name.empty()) {
287 if (proj_string.empty()) {
293 if (proj_string.empty()) {
298 if (proj_string.empty()) {
320 auto mapping_variable_name = input_file.
read_text_attribute(variable_name,
"grid_mapping");
345 e.
add_context(
"getting projection info from variable '%s' in '%s'", variable_name.c_str(),
346 input_file.
name().c_str());
372#if (Pism_USE_PROJ==1)
375static double triangle_area(
const double *A,
const double *B,
const double *C) {
377 for (
int j = 0; j < 3; ++j) {
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));
389 auto grid = result.grid();
391 Proj pism_to_geocent(projection,
"+proj=geocent +datum=WGS84 +ellps=WGS84");
404 double dx2 = 0.5 *
grid->dx(), dy2 = 0.5 *
grid->dy();
406 array::AccessScope list(result);
408 for (
auto p =
grid->points(); p; p.next()) {
409 const int i = p.i(), j = p.j();
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;
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};
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};
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};
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};
443 LonLat which, array::Scalar &result) {
445 Proj crs(projection,
"EPSG:4326");
458 auto grid = result.grid();
460 array::AccessScope list{&result};
462 for (
auto p =
grid->points(); p; p.next()) {
463 const int i = p.i(), j = p.j();
468 out = proj_trans(crs, PJ_FWD, in);
471 result(i, j) = out.lp.phi;
473 result(i, j) = out.lp.lam;
480 array::Array3D &result) {
482 Proj crs(projection,
"EPSG:4326");
484 auto grid = result.grid();
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};
490 array::AccessScope list{&result};
492 for (
auto p =
grid->points(); p; p.next()) {
493 const int i = p.i(), j = p.j();
495 double x0 =
grid->x(i), y0 =
grid->y(j);
497 double *values = result.get_column(i, j);
499 for (
int k = 0;
k < 4; ++
k) {
503 in.xy = {x0 + x_offsets[
k], y0 + y_offsets[
k]};
506 out = proj_trans(crs, PJ_FWD, in);
509 values[
k] = out.lp.lam;
511 values[
k] = out.lp.phi;
522 auto grid = result.
grid();
523 result.
set(grid->dx() * grid->dy());
533 " Please rebuild PISM with PROJ.");
544 " Please rebuild PISM with PROJ.");
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"];
575 std::vector<double> standard_parallel = mapping[
"standard_parallel"];
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,
582 standard_parallel[0]);
583 if (standard_parallel.size() == 2) {
584 result +=
pism::printf(
" +lat_2=%f", standard_parallel[1]);
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);
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);
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"];
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,
636 standard_parallel[0]);
637 if (standard_parallel.size() == 2) {
638 result +=
pism::printf(
" +lat_2=%f", standard_parallel[1]);
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"];
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);
658 result +=
pism::printf(
" +lat_ts=%f", (
double)mapping[
"standard_parallel"]);
660 result +=
pism::printf(
" +k_0=%f", (
double)mapping[
"scale_factor_at_projection_origin"]);
668 return "+proj=lonlat +type=crs";
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"];
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);
687 result +=
pism::printf(
" +lat_ts=%f", (
double)mapping[
"standard_parallel"]);
689 result +=
pism::printf(
" +k_0=%f", (
double)mapping[
"scale_factor_at_projection_origin"]);
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"];
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,
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"];
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,
736 result +=
pism::printf(
" +lat_ts=%f", (
double)mapping[
"standard_parallel"]);
738 result +=
pism::printf(
" +k_0=%f", (
double)mapping[
"scale_factor_at_projection_origin"]);
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;
767 north_pole_grid_longitude = mapping[
"north_pole_grid_longitude"];
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);
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"];
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,
797 scale_factor_at_projection_origin);
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"];
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,
823 std::string name = mapping[
"grid_mapping_name"];
831 auto flag = [&mapping](
const char* name,
const char* proj_name) {
833 return pism::printf(
" +%s=%f", proj_name, (
double)mapping[name]);
835 return std::string{};
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");
847 result += flag(
"longitude_of_prime_meridian",
"pm");
866 std::string wkt = mapping[
"crs_wkt"];
867 if (not wkt.empty()) {
873 std::map<std::string, cf_to_proj_fn> mappings = {
889 std::string mapping_name = mapping[
"grid_mapping_name"];
891 if (mappings.find(mapping_name) == mappings.end()) {
893 "conversion CF -> PROJ for a '%s' grid is not implemented",
894 mapping_name.c_str());
897 return mappings[mapping_name](mapping) +
common_flags(mapping);
903 if (mapping.has_attributes()) {
AxisType dimension_type(const std::string &name, units::System::Ptr unit_system) const
Get the "type" of a dimension.
bool variable_exists(const std::string &short_name) const
Checks if a variable exists.
void define_variable(const std::string &name, io::Type nctype, const std::vector< std::string > &dims) const
Define a variable.
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.
std::vector< std::string > dimensions(const std::string &variable_name) const
std::string read_text_attribute(const std::string &var_name, const std::string &att_name) const
Get a text attribute.
High-level PISM I/O class.
MappingInfo(const std::string &mapping_variable_name, units::System::Ptr unit_system)
VariableMetadata cf_mapping
grid mapping description following CF conventions
std::string proj_string
a projection definition string in a format recognized by PROJ 6.x+
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
A virtual class collecting methods common to ice and bedrock 3D fields.
std::shared_ptr< const Grid > grid() const
void set(double c)
Result: v[j] <- c for all j.
std::shared_ptr< System > Ptr
#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.
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)
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)
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.
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)
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)