36 std::vector<MeshLib::Node*>
const& nodes,
37 double const& max_dist)
39 double sqr_shortest_dist(max_dist);
40 double elevation(p[2]);
43 double sqr_dist = (p[0] - (*node)[0]) * (p[0] - (*node)[0]) +
44 (p[1] - (*node)[1]) * (p[1] - (*node)[1]);
45 if (sqr_dist < sqr_shortest_dist)
47 sqr_shortest_dist = sqr_dist;
48 elevation = (*node)[2];
54int main(
int argc,
char* argv[])
57 "Changes the elevation of 2D mesh nodes based on either raster data or "
58 "another 2D mesh. In addition, a low pass filter can be applied to "
59 "node elevation based connected nodes.\n\n"
60 "OpenGeoSys-6 software, version " +
63 "Copyright (c) 2012-2025, OpenGeoSys Community "
64 "(http://www.opengeosys.org)",
67 cmd.add(log_level_arg);
68 TCLAP::SwitchArg lowpass_arg(
70 "Applies a lowpass filter to elevation over connected nodes.",
false);
72 TCLAP::ValueArg<double> map_static_arg(
74 "Static elevation to map the input file to. This can be combined with "
75 "mapping based on rasters or other meshes to deal with locations where "
76 "no corresponding data exists.",
77 false, 0,
"MAP_STATIC");
78 cmd.add(map_static_arg);
79 TCLAP::ValueArg<double> max_dist_arg(
81 "Maximum distance to search for mesh nodes if there is no "
82 "corresponding data for input mesh nodes on the mesh it should be "
83 "mapped on. (min = 0)",
84 false, 1,
"MAX_DISTANCE");
85 cmd.add(max_dist_arg);
86 TCLAP::ValueArg<std::string> map_mesh_arg(
87 "m",
"mesh",
"Input (.vtu). 2D mesh file to map the input file on",
88 false,
"",
"INPUT_FILE");
89 cmd.add(map_mesh_arg);
90 TCLAP::ValueArg<std::string> map_raster_arg(
92 "Input (.asc | .grd | .xyz) Raster file to map the input file on",
93 false,
"",
"INPUT_FILE");
94 cmd.add(map_raster_arg);
95 TCLAP::ValueArg<std::string> output_arg(
96 "o",
"output",
"Output (.vtu) mesh file",
true,
"",
"OUTPUT_FILE");
98 TCLAP::ValueArg<std::string> input_arg(
99 "i",
"input",
"Input (.vtu | .msh) mesh file",
true,
"",
"INPUT_FILE");
101 cmd.parse(argc, argv);
106 std::unique_ptr<MeshLib::Mesh> mesh(
110 ERR(
"Error reading mesh file.");
114 if (!(map_static_arg.isSet() || map_raster_arg.isSet() ||
115 map_mesh_arg.isSet()))
117 ERR(
"Nothing to do. Please choose mapping based on a raster or mesh "
118 "file, or to a static value.");
122 if (map_raster_arg.isSet() && map_mesh_arg.isSet())
124 ERR(
"Please select mapping based on *either* a mesh or a raster file.");
129 if (map_static_arg.isSet())
132 *mesh, map_static_arg.getValue());
136 if (map_raster_arg.isSet())
138 std::string
const raster_path = map_raster_arg.getValue();
139 std::ifstream file_stream(raster_path, std::ifstream::in);
140 if (!file_stream.good())
142 ERR(
"Opening raster file {} failed.", raster_path);
147 std::unique_ptr<GeoLib::Raster>
const raster(
153 if (map_mesh_arg.isSet())
155 std::unique_ptr<MeshLib::Mesh> ground_truth(
157 if (ground_truth ==
nullptr)
159 ERR(
"Error reading mesh file.");
163 std::vector<MeshLib::Node*>
const& nodes = mesh->getNodes();
165 auto const edgeLengths = minMaxEdgeLength(mesh->getElements());
166 double const max_edge = edgeLengths.second;
167 double const max_dist(pow(max_dist_arg.getValue(), 2));
172 (*node)[1] - max_edge,
173 -std::numeric_limits<double>::max()}};
175 (*node)[1] + max_edge,
176 std::numeric_limits<double>::max()}};
177 std::vector<const MeshLib::Element*>
const& elems =
179 auto const* element =
182 if (element !=
nullptr)
190 *node, ground_truth->getNodes(), max_dist);
198 if (lowpass_arg.isSet())
201 const std::size_t nNodes(mesh->getNumberOfNodes());
202 std::vector<MeshLib::Node*> nodes(mesh->getNodes());
204 std::vector<double> elevation(nNodes);
205 for (std::size_t i = 0; i < nNodes; i++)
207 elevation[i] = (*nodes[i])[2];
210 auto const& connections =
212 for (std::size_t i = 0; i < nNodes; i++)
214 auto const& conn_nodes(connections[nodes[i]->getID()]);
215 const unsigned nConnNodes(conn_nodes.size());
216 elevation[i] = (2 * (*nodes[i])[2]);
217 for (std::size_t j = 0; j < nConnNodes; ++j)
219 elevation[i] += (*conn_nodes[j])[2];
221 elevation[i] /= (nConnNodes + 2);
224 for (std::size_t i = 0; i < nNodes; i++)
226 (*nodes[i])[2] = elevation[i];
235 INFO(
"Result successfully written.");