29{
30 TCLAP::CmdLine cmd(
31 "Mesh to raster converter.\n"
32 "Rasterises a 2D mesh, pixel values are set to the elevation of a "
33 "regular grid superimposed on the mesh. If no mesh element is located "
34 "beneath a pixel it is set to NODATA.\n\n"
35 "OpenGeoSys-6 software, version " +
37 ".\n"
38 "Copyright (c) 2012-2025, OpenGeoSys Community "
39 "(http://www.opengeosys.org)",
41 TCLAP::ValueArg<double> cell_arg(
42 "c", "cellsize", "edge length of raster cells in result. (min = 0)",
43 false, 1, "CELLSIZE");
44 cmd.add(cell_arg);
45 TCLAP::ValueArg<std::string> output_arg("o", "output-file",
46 "Output (.asc). Raster output file",
47 true, "", "OUTPUT_FILE");
48 cmd.add(output_arg);
49 TCLAP::ValueArg<std::string> input_arg(
50 "i", "input-file", "Input (.vtu | .msh). Mesh input file", true, "",
51 "string");
52 cmd.add(input_arg);
54 cmd.add(log_level_arg);
55 cmd.parse(argc, argv);
56
59
60 INFO(
"Rasterising mesh...");
61 std::unique_ptr<MeshLib::Mesh> const mesh(
63 if (mesh == nullptr)
64 {
65 ERR(
"Error reading mesh file.");
66 return 1;
67 }
68 if (mesh->getDimension() != 2)
69 {
70 ERR(
"The programme requires a mesh containing two-dimensional elements "
71 "(i.e. triangles or quadrilaterals.");
72 return 2;
73 }
74
76 double const cellsize =
77 (cell_arg.isSet()) ? cell_arg.getValue() : edgeLengths.first;
78 INFO(
"Cellsize set to {:f}", cellsize);
79
80 double const max_edge = edgeLengths.second + cellsize;
81
82 std::vector<MeshLib::Node*> const& nodes_vec(mesh->getNodes());
83 GeoLib::AABB const bounding_box(nodes_vec.begin(), nodes_vec.end());
84 auto const& [min, max] = bounding_box.getMinMaxPoints();
85 auto const n_cols =
86 static_cast<std::size_t>(std::ceil((max[0] - min[0]) / cellsize));
87 auto const n_rows =
88 static_cast<std::size_t>(std::ceil((max[1] - min[1]) / cellsize));
89 double const half_cell = cellsize / 2.0;
90
91
92 std::string output_name = output_arg.getValue();
93 std::string const ext = ".asc";
94 if (std::filesystem::path(output_name).extension() != ext)
95 {
96 WARN(
"Adding extension '*.asc' to output file name '{:s}'.",
97 output_name);
98 output_name += ext;
99 }
100 std::ofstream out(output_name);
101 out << "ncols " << n_cols << "\n";
102 out << "nrows " << n_rows << "\n";
103 out << std::fixed << "xllcorner " << (min[0] - half_cell) << "\n";
104 out << std::fixed << "yllcorner " << (min[1] - half_cell) << "\n";
105 out << std::fixed << "cellsize " << cellsize << "\n";
106 out << "NODATA_value "
107 << "-9999\n";
108 INFO(
"Writing raster with {:d} x {:d} pixels.", n_cols, n_rows);
109
111
112 for (std::size_t row = 0; row < n_rows; ++row)
113 {
114 double const y = max[1] - row * cellsize;
115 for (std::size_t column = 0; column < n_cols; ++column)
116 {
117
118 double const x = min[0] + column * cellsize;
121 -std::numeric_limits<double>::max()}};
123 std::numeric_limits<double>::max()}};
124 std::vector<const MeshLib::Element*> const& elems =
125 grid.getElementsInVolume(min_vol, max_vol);
126 auto const* element =
128 node);
129
130 if (element != nullptr)
131 {
133 node)
134 << " ";
135 }
136 else
137 {
138 std::array<double, 4> const x_off{
139 {-half_cell, half_cell, -half_cell, half_cell}};
140 std::array<double, 4> const y_off{
141 {-half_cell, -half_cell, half_cell, half_cell}};
142 double sum(0);
143 std::size_t nonzero_count(0);
144
145
146 for (std::size_t i = 0; i < 4; ++i)
147 {
149 0);
150 auto const* corner_element =
152 elems, corner_node);
153 if (corner_element != nullptr)
154 {
156 *corner_element, node);
157 nonzero_count++;
158 }
159 }
160 if (nonzero_count > 0)
161 {
162
163 out << sum / nonzero_count << " ";
164 }
165 else
166 {
167
168 out << "-9999 ";
169 }
170 }
171 }
172 out << "\n";
173 }
174 out.close();
175 INFO(
"Result written to {:s}", output_name);
176 return 0;
177}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
std::pair< double, double > minMaxEdgeLength(std::vector< Element * > const &elements)
Returns the minimum and maximum edge length for given elements.