OGS
Mesh2Raster.cpp
Go to the documentation of this file.
1
10#include <tclap/CmdLine.h>
11
12#ifdef USE_PETSC
13#include <mpi.h>
14#endif
15
16#include <filesystem>
17#include <fstream>
18#include <memory>
19#include <string>
20
21#include "GeoLib/AABB.h"
22#include "InfoLib/GitInfo.h"
24#include "MeshLib/Mesh.h"
26#include "MeshLib/Node.h"
28
29int main(int argc, char* argv[])
30{
31 TCLAP::CmdLine cmd(
32 "Mesh to raster converter.\n"
33 "Rasterises a 2D mesh, pixel values are set to the elevation of a "
34 "regular grid superimposed on the mesh. If no mesh element is located "
35 "beneath a pixel it is set to NODATA.\n\n"
36 "OpenGeoSys-6 software, version " +
38 ".\n"
39 "Copyright (c) 2012-2024, OpenGeoSys Community "
40 "(http://www.opengeosys.org)",
42 TCLAP::ValueArg<double> cell_arg("c", "cellsize",
43 "edge length of raster cells in result",
44 false, 1, "real");
45 cmd.add(cell_arg);
46 TCLAP::ValueArg<std::string> output_arg(
47 "o", "output-file", "Raster output file (*.asc)", true, "", "string");
48 cmd.add(output_arg);
49 TCLAP::ValueArg<std::string> input_arg("i", "input-file",
50 "Mesh input file (*.vtu, *.msh)",
51 true, "", "string");
52 cmd.add(input_arg);
53 cmd.parse(argc, argv);
54
55#ifdef USE_PETSC
56 MPI_Init(&argc, &argv);
57#endif
58
59 INFO("Rasterising mesh...");
60 std::unique_ptr<MeshLib::Mesh> const mesh(
61 MeshLib::IO::readMeshFromFile(input_arg.getValue()));
62 if (mesh == nullptr)
63 {
64 ERR("Error reading mesh file.");
65#ifdef USE_PETSC
66 MPI_Finalize();
67#endif
68 return 1;
69 }
70 if (mesh->getDimension() != 2)
71 {
72 ERR("The programme requires a mesh containing two-dimensional elements "
73 "(i.e. triangles or quadrilaterals.");
74#ifdef USE_PETSC
75 MPI_Finalize();
76#endif
77 return 2;
78 }
79
80 auto const edgeLengths = minMaxEdgeLength(mesh->getElements());
81 double const cellsize =
82 (cell_arg.isSet()) ? cell_arg.getValue() : edgeLengths.first;
83 INFO("Cellsize set to {:f}", cellsize);
84
85 double const max_edge = edgeLengths.second + cellsize;
86
87 std::vector<MeshLib::Node*> const& nodes_vec(mesh->getNodes());
88 GeoLib::AABB const bounding_box(nodes_vec.begin(), nodes_vec.end());
89 auto const& [min, max] = bounding_box.getMinMaxPoints();
90 auto const n_cols =
91 static_cast<std::size_t>(std::ceil((max[0] - min[0]) / cellsize));
92 auto const n_rows =
93 static_cast<std::size_t>(std::ceil((max[1] - min[1]) / cellsize));
94 double const half_cell = cellsize / 2.0;
95
96 // raster header
97 std::string output_name = output_arg.getValue();
98 std::string const ext = ".asc";
99 if (std::filesystem::path(output_name).extension() != ext)
100 {
101 WARN("Adding extension '*.asc' to output file name '{:s}'.",
102 output_name);
103 output_name += ext;
104 }
105 std::ofstream out(output_name);
106 out << "ncols " << n_cols << "\n";
107 out << "nrows " << n_rows << "\n";
108 out << std::fixed << "xllcorner " << (min[0] - half_cell) << "\n";
109 out << std::fixed << "yllcorner " << (min[1] - half_cell) << "\n";
110 out << std::fixed << "cellsize " << cellsize << "\n";
111 out << "NODATA_value "
112 << "-9999\n";
113 INFO("Writing raster with {:d} x {:d} pixels.", n_cols, n_rows);
114
115 MeshLib::MeshElementGrid const grid(*mesh);
116
117 for (std::size_t row = 0; row < n_rows; ++row)
118 {
119 double const y = max[1] - row * cellsize;
120 for (std::size_t column = 0; column < n_cols; ++column)
121 {
122 // pixel values
123 double const x = min[0] + column * cellsize;
124 MeshLib::Node const node(x, y, 0);
125 MathLib::Point3d min_vol{{x - max_edge, y - max_edge,
126 -std::numeric_limits<double>::max()}};
127 MathLib::Point3d max_vol{{x + max_edge, y + max_edge,
128 std::numeric_limits<double>::max()}};
129 std::vector<const MeshLib::Element*> const& elems =
130 grid.getElementsInVolume(min_vol, max_vol);
131 auto const* element =
133 node);
134 // centre of the pixel is located within a mesh element
135 if (element != nullptr)
136 {
138 node)
139 << " ";
140 }
141 else
142 {
143 std::array<double, 4> const x_off{
144 {-half_cell, half_cell, -half_cell, half_cell}};
145 std::array<double, 4> const y_off{
146 {-half_cell, -half_cell, half_cell, half_cell}};
147 double sum(0);
148 std::size_t nonzero_count(0);
149 // test all of the pixel's corners if there are any within an
150 // element
151 for (std::size_t i = 0; i < 4; ++i)
152 {
153 MeshLib::Node const corner_node(x + x_off[i], y + y_off[i],
154 0);
155 auto const* corner_element =
157 elems, corner_node);
158 if (corner_element != nullptr)
159 {
161 *corner_element, node);
162 nonzero_count++;
163 }
164 }
165 if (nonzero_count > 0)
166 {
167 // calculate pixel value as average of corner values
168 out << sum / nonzero_count << " ";
169 }
170 else
171 {
172 // if none of the corners give a value, set pixel to NODATA
173 out << "-9999 ";
174 }
175 }
176 }
177 out << "\n";
178 }
179 out.close();
180 INFO("Result written to {:s}", output_name);
181#ifdef USE_PETSC
182 MPI_Finalize();
183#endif
184 return 0;
185}
Definition of the AABB class.
Git information.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
int main(int argc, char *argv[])
Definition of the Mesh class.
Definition of the Node class.
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
MinMaxPoints getMinMaxPoints() const
Definition AABB.h:174
std::vector< MeshLib::Element const * > getElementsInVolume(POINT const &min, POINT const &max) const
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
double getElevation(MeshLib::Element const &element, MathLib::Point3d const &node)
MeshLib::Element const * getProjectedElement(std::vector< const MeshLib::Element * > const &elements, MathLib::Point3d const &node)
Definition of readMeshFromFile function.