29int main(
int argc,
char* argv[])
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 " +
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",
46 TCLAP::ValueArg<std::string> output_arg(
47 "o",
"output-file",
"Raster output file (*.asc)",
true,
"",
"string");
49 TCLAP::ValueArg<std::string> input_arg(
"i",
"input-file",
50 "Mesh input file (*.vtu, *.msh)",
53 cmd.parse(argc, argv);
56 MPI_Init(&argc, &argv);
59 INFO(
"Rasterising mesh...");
60 std::unique_ptr<MeshLib::Mesh>
const mesh(
64 ERR(
"Error reading mesh file.");
70 if (mesh->getDimension() != 2)
72 ERR(
"The programme requires a mesh containing two-dimensional elements "
73 "(i.e. triangles or quadrilaterals.");
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);
85 double const max_edge = edgeLengths.second + cellsize;
87 std::vector<MeshLib::Node*>
const& nodes_vec(mesh->getNodes());
88 GeoLib::AABB const bounding_box(nodes_vec.begin(), nodes_vec.end());
91 static_cast<std::size_t
>(std::ceil((max[0] - min[0]) / cellsize));
93 static_cast<std::size_t
>(std::ceil((max[1] - min[1]) / cellsize));
94 double const half_cell = cellsize / 2.0;
97 std::string output_name = output_arg.getValue();
98 std::string
const ext =
".asc";
99 if (std::filesystem::path(output_name).extension() != ext)
101 WARN(
"Adding extension '*.asc' to output file name '{:s}'.",
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 "
113 INFO(
"Writing raster with {:d} x {:d} pixels.", n_cols, n_rows);
117 for (std::size_t row = 0; row < n_rows; ++row)
119 double const y = max[1] - row * cellsize;
120 for (std::size_t column = 0; column < n_cols; ++column)
123 double const x = min[0] + column * cellsize;
126 -std::numeric_limits<double>::max()}};
128 std::numeric_limits<double>::max()}};
129 std::vector<const MeshLib::Element*>
const& elems =
131 auto const* element =
135 if (element !=
nullptr)
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}};
148 std::size_t nonzero_count(0);
151 for (std::size_t i = 0; i < 4; ++i)
155 auto const* corner_element =
158 if (corner_element !=
nullptr)
161 *corner_element, node);
165 if (nonzero_count > 0)
168 out << sum / nonzero_count <<
" ";
180 INFO(
"Result written to {:s}", output_name);