27 "Mesh to raster converter.\n"
28 "Rasterises a 2D mesh, pixel values are set to the elevation of a "
29 "regular grid superimposed on the mesh. If no mesh element is located "
30 "beneath a pixel it is set to NODATA.\n\n"
31 "OpenGeoSys-6 software, version " +
34 "Copyright (c) 2012-2021, OpenGeoSys Community "
35 "(http://www.opengeosys.org)",
37 TCLAP::ValueArg<double> cell_arg(
"c",
"cellsize",
38 "edge length of raster cells in result",
41 TCLAP::ValueArg<std::string> output_arg(
42 "o",
"output-file",
"Raster output file (*.asc)",
true,
"",
"string");
44 TCLAP::ValueArg<std::string> input_arg(
"i",
"input-file",
45 "Mesh input file (*.vtu, *.msh)",
48 cmd.parse(argc, argv);
50 INFO(
"Rasterising mesh...");
51 std::unique_ptr<MeshLib::Mesh>
const mesh(
55 ERR(
"Error reading mesh file.");
58 if (mesh->getDimension() != 2)
60 ERR(
"The programme requires a mesh containing two-dimensional elements "
61 "(i.e. triangles or quadrilaterals.");
65 double const cellsize =
66 (cell_arg.isSet()) ? cell_arg.getValue() : mesh->getMinEdgeLength();
67 INFO(
"Cellsize set to {:f}", cellsize);
69 std::vector<MeshLib::Node*>
const& nodes_vec(mesh->getNodes());
70 GeoLib::AABB const bounding_box(nodes_vec.begin(), nodes_vec.end());
71 auto const& min(bounding_box.getMinPoint());
72 auto const& max(bounding_box.getMaxPoint());
74 static_cast<std::size_t
>(std::ceil((max[0] - min[0]) / cellsize));
76 static_cast<std::size_t
>(std::ceil((max[1] - min[1]) / cellsize));
77 double const half_cell = cellsize / 2.0;
80 std::ofstream out(output_arg.getValue());
81 out <<
"ncols " << n_cols + 1 <<
"\n";
82 out <<
"nrows " << n_rows + 1 <<
"\n";
83 out << std::fixed <<
"xllcorner " << (min[0] - half_cell) <<
"\n";
84 out << std::fixed <<
"yllcorner " << (min[1] - half_cell) <<
"\n";
85 out << std::fixed <<
"cellsize " << cellsize <<
"\n";
86 out <<
"NODATA_value "
88 INFO(
"Writing raster with {:d} x {:d} pixels.", n_cols, n_rows);
91 double const max_edge(mesh->getMaxEdgeLength() + cellsize);
93 for (std::size_t row = 0; row <= n_rows; ++row)
95 double const y = max[1] - row * cellsize;
96 for (std::size_t column = 0; column <= n_cols; ++column)
99 double const x = min[0] + column * cellsize;
102 -std::numeric_limits<double>::max()}};
104 std::numeric_limits<double>::max()}};
105 std::vector<const MeshLib::Element*>
const& elems =
106 grid.getElementsInVolume(min_vol, max_vol);
107 auto const* element =
110 if (element !=
nullptr)
117 std::array<double, 4>
const x_off{
118 {-half_cell, half_cell, -half_cell, half_cell}};
119 std::array<double, 4>
const y_off{
120 {-half_cell, -half_cell, half_cell, half_cell}};
122 std::size_t nonzero_count(0);
125 for (std::size_t i = 0; i < 4; ++i)
129 auto const* corner_element =
132 if (corner_element !=
nullptr)
135 *corner_element, corner_node);
139 if (nonzero_count > 0)
142 out << sum / nonzero_count <<
" ";
154 INFO(
"Result written to {:s}", output_arg.getValue());
void INFO(char const *fmt, Args const &... args)
void ERR(char const *fmt, Args const &... args)
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name)
double getElevation(Element const &element, Node const &node)
Element const * getProjectedElement(std::vector< const Element * > const &elements, Node const &node)