OGS
Mesh2Raster.cpp File Reference

Detailed Description

Definition in file Mesh2Raster.cpp.

#include <tclap/CmdLine.h>
#include <filesystem>
#include <fstream>
#include <memory>
#include <string>
#include "BaseLib/Logging.h"
#include "BaseLib/MPI.h"
#include "BaseLib/TCLAPArguments.h"
#include "GeoLib/AABB.h"
#include "InfoLib/GitInfo.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshSearch/MeshElementGrid.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshEditing/ProjectPointOnMesh.h"
Include dependency graph for Mesh2Raster.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 28 of file Mesh2Raster.cpp.

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);
53 auto log_level_arg = BaseLib::makeLogLevelArg();
54 cmd.add(log_level_arg);
55 cmd.parse(argc, argv);
56
57 BaseLib::MPI::Setup mpi_setup(argc, argv);
58 BaseLib::initOGSLogger(log_level_arg.getValue());
59
60 INFO("Rasterising mesh...");
61 std::unique_ptr<MeshLib::Mesh> const mesh(
62 MeshLib::IO::readMeshFromFile(input_arg.getValue()));
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
75 auto const edgeLengths = minMaxEdgeLength(mesh->getElements());
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 // raster header
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
110 MeshLib::MeshElementGrid const grid(*mesh);
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 // pixel values
118 double const x = min[0] + column * cellsize;
119 MeshLib::Node const node(x, y, 0);
120 MathLib::Point3d min_vol{{x - max_edge, y - max_edge,
121 -std::numeric_limits<double>::max()}};
122 MathLib::Point3d max_vol{{x + max_edge, y + max_edge,
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 // centre of the pixel is located within a mesh element
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 // test all of the pixel's corners if there are any within an
145 // element
146 for (std::size_t i = 0; i < 4; ++i)
147 {
148 MeshLib::Node const corner_node(x + x_off[i], y + y_off[i],
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 // calculate pixel value as average of corner values
163 out << sum / nonzero_count << " ";
164 }
165 else
166 {
167 // if none of the corners give a value, set pixel to NODATA
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)
Definition Logging.h:36
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:42
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:64
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.
Definition Mesh.cpp:190
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)

References ERR(), MeshLib::MeshElementGrid::getElementsInVolume(), MeshToolsLib::ProjectPointOnMesh::getElevation(), GeoLib::AABB::getMinMaxPoints(), MeshToolsLib::ProjectPointOnMesh::getProjectedElement(), INFO(), BaseLib::initOGSLogger(), BaseLib::makeLogLevelArg(), GitInfoLib::GitInfo::ogs_version, MeshLib::IO::readMeshFromFile(), and WARN().