OGS
NetCdfConverter.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// STL
17#include <cctype>
18#include <iostream>
19#include <limits>
20#include <memory>
21#include <netcdf>
22#include <numeric>
23#include <sstream>
24#include <string>
25#include <utility>
26
27#include "BaseLib/FileTools.h"
28#include "BaseLib/Logging.h"
30#include "GeoLib/Raster.h"
31#include "InfoLib/GitInfo.h"
33#include "MeshLib/Mesh.h"
36
37using namespace netCDF;
38
39static const double no_data_output = -9999;
40static double no_data_input = -9999;
41
42enum class OutputType
43{
44 INVALID,
45 IMAGES,
48};
49
50static void checkExit(std::string const& str)
51{
52 if (str == "x" || str == "exit")
53 exit(0);
54}
55
56static void showErrorMessage(std::size_t const error_id,
57 std::size_t const max = 0)
58{
59 if (error_id == 0)
60 {
61 ERR("Input not valid.");
62 }
63 else if (error_id == 1)
64 {
65 ERR("Index not valid. Valid indices are in [0,{:d}].", max);
66 }
67 else if (error_id == 2)
68 {
69 ERR("Input not valid.");
70 std::cout << "Type \"info\" to display the available options again. "
71 "\"exit\" will exit the programme.\n";
72 }
73}
74
75static std::size_t parseInput(std::string const& request_str,
76 std::size_t const max_val,
77 bool const has_info = false)
78{
79 while (true)
80 {
81 std::cout << request_str;
82 std::string str;
83 std::size_t val;
84 std::getline(std::cin, str);
85 checkExit(str);
86 if (has_info && str == "info")
87 return (max_val);
88 std::stringstream str_stream(str);
89 if (!(str_stream >> val))
90 {
91 std::size_t const error_val = (has_info) ? 2 : 0;
92 showErrorMessage(error_val);
93 continue;
94 }
95 if (val > max_val - 1)
96 {
97 showErrorMessage(1, max_val - 1);
98 continue;
99 }
100 return val;
101 }
102 return std::numeric_limits<std::size_t>::max();
103}
104
105static NcVar getDimVar(NcFile const& dataset, NcVar const& var,
106 std::size_t const dim)
107{
108 NcDim const& dim_obj = var.getDim(dim);
109 return dataset.getVar(dim_obj.getName());
110}
111
112static std::pair<double, double> getDimLength(NcVar const& var,
113 std::size_t const dim)
114{
115 return std::make_pair(0.0, static_cast<double>(var.getDim(dim).getSize()));
116}
117
118static std::vector<std::string> getArrays(NcFile const& dataset)
119{
120 auto const& names = dataset.getVars();
121 std::vector<std::string> var_names;
122 for (auto [name, var] : names)
123 {
124 (void)var;
125 var_names.push_back(name);
126 }
127 return var_names;
128}
129
130static void showArrays(NcFile const& dataset)
131{
132 std::size_t const n_vars(dataset.getDimCount());
133 std::cout << "The NetCDF file contains the following " << n_vars
134 << " arrays:\n\n";
135 std::cout << "\tIndex\tArray Name\t#Dimensions\n";
136 std::cout << "-------------------------------------------\n";
137 auto const& names = dataset.getVars();
138 std::size_t count = 0;
139 for (auto [name, var] : names)
140 {
141 std::cout << "\t" << count++ << "\t" << name << "\t("
142 << var.getDimCount() << "D array)\n";
143 }
144 std::cout << "\n";
145}
146
147static void showArraysDims(NcVar const& var)
148{
149 std::cout << "Data array \"" << var.getName()
150 << "\" contains the following dimensions:\n";
151 std::size_t const n_dims(var.getDimCount());
152 for (std::size_t i = 0; i < n_dims; ++i)
153 std::cout << "\t" << i << "\t" << var.getDim(i).getName() << "\t("
154 << var.getDim(i).getSize() << " values)\n";
155 std::cout << "\n";
156}
157
158static std::pair<double, double> getBoundaries(NcVar const& var)
159{
160 if (var.getDimCount() == 1)
161 {
162 double start, end;
163 std::size_t const size = var.getDim(0).getSize();
164 var.getVar({0}, {1}, &start);
165 var.getVar({size - 1}, {1}, &end);
166 return std::make_pair(start, end);
167 }
168 return std::make_pair(0, 0);
169}
170
171static MathLib::Point3d getOrigin(NcFile const& dataset, NcVar const& var,
172 std::vector<std::size_t> const& dim_idx_map,
173 std::size_t const time_offset)
174{
176 std::size_t const n_dims = var.getDimCount();
177 for (std::size_t i = time_offset; i < n_dims; ++i)
178 {
179 NcVar const& dim = getDimVar(dataset, var, dim_idx_map[i]);
180 auto const bounds = (dim.isNull()) ? getDimLength(var, dim_idx_map[i])
181 : getBoundaries(dim);
182 origin[i - time_offset] =
183 (bounds.first < bounds.second) ? bounds.first : bounds.second;
184 }
185 return origin;
186}
187
188static void flipRaster(std::vector<double>& data, std::size_t const layers,
189 std::size_t const width, std::size_t const height)
190{
191 std::size_t const length(data.size());
192 std::vector<double> tmp_vec;
193 tmp_vec.reserve(length);
194 for (std::size_t k = 0; k < layers; k++)
195 {
196 std::size_t const layer_end = (k + 1) * height * width;
197 for (std::size_t i = 0; i < height; i++)
198 {
199 std::size_t const line_idx(layer_end - (width * (i + 1)));
200 for (std::size_t j = 0; j < width; j++)
201 {
202 tmp_vec.push_back(data[line_idx + j]);
203 }
204 }
205 }
206 std::copy(tmp_vec.cbegin(), tmp_vec.cend(), data.begin());
207}
208
209static bool canConvert(NcVar const& var)
210{
211 bool ret(var.getDimCount() < 2);
212 if (ret)
213 ERR("Only 2+ dimensional variables can be converted into OGS "
214 "Meshes.\n");
215 return !ret;
216}
217
218static std::string arraySelectionLoop(NcFile const& dataset)
219{
220 std::vector<std::string> const& names = getArrays(dataset);
221 showArrays(dataset);
222 std::size_t const idx =
223 parseInput("Enter data array index: ", dataset.getVarCount(), true);
224
225 if (static_cast<int>(idx) == dataset.getVarCount() ||
226 !canConvert(dataset.getVar(names[idx])))
227 return arraySelectionLoop(dataset);
228
229 return names[idx];
230}
231
232static bool dimensionSelectionLoop(NcVar const& var,
233 std::vector<std::size_t>& dim_idx_map)
234{
235 showArraysDims(var);
236 std::size_t const n_dims(var.getDimCount());
237 dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
238 bool is_time_dep(true);
239
240 // get temporal dimension
241 if (n_dims > 1)
242 {
243 std::string temp_str("");
244 cout << "Is the parameter time-dependent?\n";
245 while (dim_idx_map[0] == std::numeric_limits<std::size_t>::max() &&
246 is_time_dep == true)
247 {
248 cout << "Enter ID for temporal dimension or \"c\" to continue: ";
249 std::getline(std::cin, temp_str);
250 std::stringstream str_stream(temp_str);
251 if (str_stream.str() == "c" || str_stream.str() == "continue")
252 is_time_dep = false;
253 else
254 {
255 if (!(str_stream >> dim_idx_map[0]))
256 {
258 dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
259 continue;
260 }
261 if (dim_idx_map[0] > n_dims - 1)
262 {
263 showErrorMessage(1, var.getDimCount() - 1);
264 dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
265 }
266 }
267 }
268 }
269 else
270 is_time_dep = false;
271
272 // get spatial dimension(s)
273 std::size_t const start_idx = (is_time_dep) ? 1 : 0;
274 std::array<std::string, 4> const dim_comment{
275 "(x / longitude)", "(y / latitude)", "(z / height / depth)",
276 "[Error: 4-dimensional non-temporal arrays are not supported]"};
277 for (std::size_t i = start_idx; i < n_dims; ++i)
278 {
279 dim_idx_map[i] = std::numeric_limits<std::size_t>::max();
280
281 std::string const request_str("Enter ID for dimension " +
282 std::to_string(i) + " " +
283 dim_comment[i - start_idx] + ": ");
284 std::size_t const idx =
285 parseInput(request_str, var.getDimCount(), true);
286
287 if (static_cast<int>(idx) == var.getDimCount())
288 {
289 showArraysDims(var);
290 i--;
291 continue;
292 }
293 dim_idx_map[i] = idx;
294 }
295
296 return is_time_dep;
297}
298
299static std::pair<std::size_t, std::size_t> timestepSelectionLoop(
300 NcVar const& var, std::size_t const dim_idx)
301{
302 std::size_t const n_time_steps = var.getDim(dim_idx).getSize();
303 std::size_t const max_val = std::numeric_limits<std::size_t>::max();
304 std::pair<std::size_t, std::size_t> bounds(max_val, max_val);
305 std::cout << "\nThe dataset contains " << n_time_steps << " time steps.\n";
306 while (bounds.first == max_val)
307 {
308 bounds.first = parseInput(
309 "Specify first time step to export: ", n_time_steps, false);
310 }
311 while (bounds.first > bounds.second || bounds.second > n_time_steps)
312 {
313 bounds.second = parseInput(
314 "Specify last time step to export: ", n_time_steps, false);
315 }
316 return bounds;
317}
318
319static MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim)
320{
321 if (dim == 1)
323
326 {
327 std::cout << "\nSelect element type for result, choose ";
328
329 if (dim == 2)
330 std::cout << "(t)riangle or (q)uadliteral: ";
331 if (dim == 3)
332 std::cout << "(p)rism or (h)exahedron: ";
333 std::string type("");
334 std::getline(std::cin, type);
335 checkExit(type);
336 if (dim == 2)
337 {
338 if (type != "t" && type != "q" && type != "tri" && type != "quad" &&
339 type != "triangle" && type != "quatliteral")
340 continue;
341 if (type == "t" || type == "tri" || type == "triangle")
344 }
345
346 if (dim == 3)
347 {
348 if (type != "p" && type != "h" && type != "prism" &&
349 type != "hex" && type != "hexahedron")
350 continue;
351 if (type == "p" || type == "prism")
354 }
355 }
356 return t;
357}
358
360 std::pair<std::size_t, std::size_t> const& time_bounds)
361{
363 while (t == OutputType::INVALID)
364 {
365 std::size_t const n_time_steps(time_bounds.second - time_bounds.first +
366 1);
367 std::cout << "\nThe selection includes " << n_time_steps
368 << " time steps.\n";
369 std::cout << "0. Save data in " << n_time_steps
370 << " mesh files with one scalar array each.\n";
371 std::cout << "1. Save data in one mesh file with " << n_time_steps
372 << " scalar arrays.\n";
373 std::cout << "2. Save data as " << n_time_steps << " ASC images.\n";
374
375 std::size_t const ret =
376 parseInput("Select preferred method: ", 3, false);
377
378 if (ret == 0)
380 else if (ret == 1)
382 else if (ret == 2)
384 }
385 return t;
386}
387
388static std::string getIterationString(std::size_t i, std::size_t max)
389{
390 std::size_t const max_length(std::to_string(max).length());
391 std::string const current_str(std::to_string(i));
392 return std::string(max_length - current_str.length(), '0') + current_str;
393}
394
395static double getResolution(NcFile const& dataset, NcVar const& var)
396{
397 std::size_t const dim_idx = var.getDimCount() - 1;
398 NcVar const dim_var(getDimVar(dataset, var, dim_idx));
399 auto const bounds = (dim_var.isNull()) ? getDimLength(var, dim_idx)
400 : getBoundaries(dim_var);
401 std::size_t const dim_size = var.getDim(dim_idx).getSize();
402 if (dim_size == 0)
403 {
404 OGS_FATAL("Dimension '{:s}' has size 0. Aborting...",
405 var.getDim(dim_idx).getName());
406 }
407 return std::fabs(bounds.second - bounds.first) /
408 static_cast<double>(dim_size);
409}
410
412 NcFile const& dataset, NcVar const& var,
413 std::vector<std::size_t> const& dim_idx_map,
414 std::vector<std::size_t> const& length, std::size_t const time_offset)
415{
416 MathLib::Point3d const origin =
417 getOrigin(dataset, var, dim_idx_map, time_offset);
418 double const res = getResolution(dataset, var);
419 std::size_t n_dims = var.getDimCount();
420 std::size_t z_length =
421 (n_dims - time_offset == 3) ? length[dim_idx_map.back()] : 1;
422 return {length[dim_idx_map[0 + time_offset]],
423 length[dim_idx_map[1 + time_offset]],
424 z_length,
425 origin,
426 res,
428}
429
430static std::vector<std::size_t> getLength(NcVar const& var,
431 std::size_t const time_offset)
432{
433 std::size_t const n_dims = (var.getDimCount());
434 std::vector<std::size_t> length(n_dims, 1);
435 for (std::size_t i = time_offset; i < n_dims; ++i)
436 {
437 length[i] = var.getDim(i).getSize();
438 }
439 return length;
440}
441
442static std::vector<double> getData(NcVar const& var,
443 std::size_t const total_length,
444 std::size_t const time_step,
445 std::vector<std::size_t> const& length)
446{
447 std::size_t const n_dims(var.getDimCount());
448 std::vector<std::size_t> offset(n_dims, 0);
449 offset[0] = time_step;
450 std::vector<double> data_vec(total_length, 0);
451 var.getVar(offset, length, data_vec.data());
452
453 std::replace_if(
454 data_vec.begin(), data_vec.end(),
455 [](double const& x) { return x == no_data_input; }, no_data_output);
456
457 return data_vec;
458}
459
460static bool assignDimParams(NcVar const& var,
461 std::vector<std::size_t>& dim_idx_map,
462 TCLAP::ValueArg<std::size_t>& arg_dim_time,
463 TCLAP::ValueArg<std::size_t>& arg_dim1,
464 TCLAP::ValueArg<std::size_t>& arg_dim2,
465 TCLAP::ValueArg<std::size_t>& arg_dim3)
466{
467 std::size_t dim_param_count = 0;
468 if (arg_dim_time.isSet())
469 dim_param_count++;
470 if (arg_dim1.isSet())
471 dim_param_count++;
472 if (arg_dim2.isSet())
473 dim_param_count++;
474 if (arg_dim3.isSet())
475 dim_param_count++;
476
477 std::size_t const n_dims = var.getDimCount();
478 if (dim_param_count != n_dims)
479 {
480 ERR("Number of parameters set does not fit number of parameters for "
481 "specified variable.");
482 return false;
483 }
484
485 if (arg_dim_time.getValue() >= n_dims || arg_dim1.getValue() >= n_dims ||
486 arg_dim2.getValue() >= n_dims || arg_dim3.getValue() >= n_dims)
487 {
488 ERR("Maximum allowed dimension for variable \"{:s}\" is {:d}.",
489 var.getName(), n_dims - 1);
490 return false;
491 }
492
493 if (arg_dim_time.isSet())
494 dim_idx_map[0] = arg_dim_time.getValue();
495 std::size_t const temp_offset = (arg_dim_time.isSet()) ? 1 : 0;
496 dim_idx_map[0 + temp_offset] = arg_dim1.getValue();
497 dim_idx_map[1 + temp_offset] = arg_dim2.getValue();
498 if (n_dims == (3 + temp_offset))
499 dim_idx_map[2 + temp_offset] = arg_dim3.getValue();
500
501 return true;
502}
503
504static std::pair<std::size_t, std::size_t> assignTimeBounds(
505 NcVar const& var,
506 TCLAP::ValueArg<std::size_t>& arg_time_start,
507 TCLAP::ValueArg<std::size_t>& arg_time_end)
508{
509 auto const bounds = getBoundaries(var);
510 if (arg_time_start.getValue() > bounds.second)
511 {
512 ERR("Start time step larger than total number of time steps. Resetting "
513 "to 0.");
514 arg_time_start.reset();
515 }
516
517 if (!arg_time_end.isSet())
518 return {arg_time_start.getValue(), arg_time_start.getValue()};
519
520 if (arg_time_end.getValue() > bounds.second)
521 {
522 ERR("End time step larger than total number of time steps. Resetting "
523 "to starting time step");
524 return {arg_time_start.getValue(), arg_time_start.getValue()};
525 }
526
527 if (arg_time_end.getValue() < arg_time_start.getValue())
528 {
529 ERR("End time step larger than starting time step. Swapping values");
530 return {arg_time_end.getValue(), arg_time_start.getValue()};
531 }
532
533 return {arg_time_start.getValue(), arg_time_end.getValue()};
534}
535
537 TCLAP::ValueArg<std::string>& arg_elem_type)
538{
539 if (arg_elem_type.getValue() == "tri")
541 if (arg_elem_type.getValue() == "quad")
543 if (arg_elem_type.getValue() == "prism")
545 if (arg_elem_type.getValue() == "hex")
547 // this should never happen
549}
550
551static bool convert(NcFile const& dataset, NcVar const& var,
552 std::string const& output_name,
553 std::vector<std::size_t> const& dim_idx_map,
554 std::size_t const time_offset,
555 std::pair<std::size_t, std::size_t> const& time_bounds,
556 OutputType const output,
557 MeshLib::MeshElemType const elem_type)
558{
559 std::unique_ptr<MeshLib::Mesh> mesh;
560 std::vector<std::size_t> const length = getLength(var, time_offset);
561 std::size_t const array_length = std::accumulate(
562 length.cbegin(), length.cend(), 1, std::multiplies<std::size_t>());
563 for (std::size_t i = time_bounds.first; i <= time_bounds.second; ++i)
564 {
565 std::string const step_str =
566 (time_bounds.first != time_bounds.second)
567 ? std::string(" time step " + std::to_string(i))
568 : "";
569 std::cout << "Converting" << step_str << "...\n";
570 std::vector<double> data_vec = getData(var, array_length, i, length);
571
572 // reverse lines in vertical direction if file has its origin in the
573 // northwest corner
574 std::size_t const n_dims = length.size();
575 NcVar const dim_var(getDimVar(dataset, var, n_dims - 2));
576 auto const bounds = (dim_var.isNull()) ? getDimLength(var, n_dims - 2)
577 : getBoundaries(dim_var);
578 if (bounds.first > bounds.second)
579 {
580 std::size_t n_layers =
581 (length.size() - time_offset == 3) ? length[n_dims - 3] : 1;
582 flipRaster(data_vec, n_layers, length[n_dims - 1],
583 length[n_dims - 2]);
584 }
585
586 GeoLib::RasterHeader const header =
587 createRasterHeader(dataset, var, dim_idx_map, length, time_offset);
588 MeshLib::UseIntensityAs const useIntensity =
590 if (output == OutputType::MULTIMESH)
591 {
592 mesh = MeshToolsLib::RasterToMesh::convert(data_vec.data(), header,
593 elem_type, useIntensity,
594 var.getName());
595 std::string const output_file_name(
596 BaseLib::dropFileExtension(output_name) +
597 getIterationString(i, time_bounds.second) + ".vtu");
598 MeshLib::IO::VtuInterface vtu(mesh.get());
599 vtu.writeToFile(output_file_name);
600 }
601 else if (output == OutputType::SINGLEMESH)
602 {
603 std::string array_name(var.getName());
604 if (time_bounds.first != time_bounds.second)
605 array_name.append(getIterationString(i, time_bounds.second));
606 if (i == time_bounds.first) // create persistent mesh
608 data_vec.data(), header, elem_type, useIntensity,
609 array_name);
610 else // copy array to mesh
611 {
612 std::unique_ptr<MeshLib::Mesh> const temp(
613 MeshToolsLib::RasterToMesh::convert(data_vec.data(), header,
614 elem_type, useIntensity,
615 array_name));
616 MeshLib::PropertyVector<double> const* const vec =
617 temp->getProperties().getPropertyVector<double>(array_name);
618 if (vec == nullptr)
619 return false;
620 MeshLib::addPropertyToMesh<double>(
621 *mesh, array_name, MeshLib::MeshItemType::Cell, 1, *vec);
622 }
623 if (i == time_bounds.second)
624 {
625 MeshLib::IO::VtuInterface vtu(mesh.get());
626 std::string const output_file_name =
627 (BaseLib::getFileExtension(output_name) == ".vtu")
628 ? output_name
629 : output_name + ".vtu";
630 vtu.writeToFile(output_file_name);
631 }
632 }
633 else // OutputType::IMAGES
634 {
635 GeoLib::Raster const raster(header, data_vec.data(),
636 data_vec.data() + header.n_cols *
637 header.n_rows *
638 header.n_depth);
639 std::string const output_file_name(
640 BaseLib::dropFileExtension(output_name) +
641 getIterationString(i, time_bounds.second) + ".asc");
643 output_file_name);
644 }
645 }
646 return true;
647}
648
649int main(int argc, char* argv[])
650{
651 TCLAP::CmdLine cmd(
652 "Converts NetCDF data into mesh file(s).\n\n "
653 "OpenGeoSys-6 software, version " +
655 ".\n"
656 "Copyright (c) 2012-2024, OpenGeoSys Community "
657 "(http://www.opengeosys.org)",
659
660 TCLAP::ValueArg<int> arg_nodata(
661 "n", "nodata",
662 "explicitly specifies the no data value used in the dataset (usually "
663 "it is not necessary to set this)",
664 false, no_data_input, "integer specifying no data value");
665 cmd.add(arg_nodata);
666
667 std::vector<std::string> allowed_elems{"tri", "quad", "prism", "hex"};
668 TCLAP::ValuesConstraint<std::string> allowed_elem_vals(allowed_elems);
669 TCLAP::ValueArg<std::string> arg_elem_type(
670 "e", "elem-type", "the element type used in the resulting OGS mesh",
671 false, "", &allowed_elem_vals);
672 cmd.add(arg_elem_type);
673
674 TCLAP::SwitchArg arg_images(
675 "", "images",
676 "if set, all time steps will be written as ESRI image files (*.asc)");
677 cmd.add(arg_images);
678
679 TCLAP::SwitchArg arg_multi_files(
680 "", "multi-file",
681 "if set, each time step will be written to a separate mesh file");
682 cmd.add(arg_multi_files);
683
684 TCLAP::SwitchArg arg_single_file(
685 "", "single-file",
686 "if set, all time steps will be written to a single mesh file (with "
687 "one scalar array per time step)");
688 cmd.add(arg_single_file);
689
690 TCLAP::ValueArg<std::size_t> arg_time_end(
691 "", "timestep-last",
692 "last time step to be extracted (only for time-dependent variables!)",
693 false, 0, "integer specifying index of time step");
694 cmd.add(arg_time_end);
695
696 TCLAP::ValueArg<std::size_t> arg_time_start(
697 "", "timestep-first",
698 "first time step to be extracted (only for time-dependent variables!)",
699 false, 0, "integer specifying index of time step");
700 cmd.add(arg_time_start);
701
702 std::vector<std::size_t> allowed_dims{0, 1, 2, 3};
703 TCLAP::ValuesConstraint<std::size_t> allowed_dim_vals(allowed_dims);
704 TCLAP::ValueArg<std::size_t> arg_dim3(
705 "", "dim3",
706 "index of third dimension (z/height/depth) for the selected variable",
707 false, 0, &allowed_dim_vals);
708 cmd.add(arg_dim3);
709
710 TCLAP::ValueArg<std::size_t> arg_dim2(
711 "", "dim2",
712 "index of second dimension (y/latitude) for the selected variable",
713 false, 0, &allowed_dim_vals);
714 cmd.add(arg_dim2);
715
716 TCLAP::ValueArg<std::size_t> arg_dim1(
717 "", "dim1",
718 "index of first dimension (x/longitude) for the selected variable",
719 false, 0, &allowed_dim_vals);
720 cmd.add(arg_dim1);
721
722 TCLAP::ValueArg<std::size_t> arg_dim_time(
723 "t", "time",
724 "index of the time-dependent dimension for the selected variable",
725 false, 0, &allowed_dim_vals);
726 cmd.add(arg_dim_time);
727
728 TCLAP::ValueArg<std::string> arg_varname(
729 "v", "var", "variable included in the the netCDF file", false, "",
730 "string containing the variable name");
731 cmd.add(arg_varname);
732
733 TCLAP::ValueArg<std::string> arg_output(
734 "o", "output", "the OGS mesh output file", true, "",
735 "string containing the path and file name");
736 cmd.add(arg_output);
737 TCLAP::ValueArg<std::string> arg_input(
738 "i", "input", "the netCDF input file", true, "",
739 "string containing the path and file name");
740 cmd.add(arg_input);
741 cmd.parse(argc, argv);
742
743#ifdef USE_PETSC
744 MPI_Init(&argc, &argv);
745#endif
746
747 NcFile dataset(arg_input.getValue().c_str(), NcFile::read);
748
749 if (dataset.isNull())
750 {
751 ERR("Error opening file.");
752#ifdef USE_PETSC
753 MPI_Finalize();
754#endif
755 return -1;
756 }
757
758 std::size_t const mutex =
759 arg_single_file.isSet() + arg_multi_files.isSet() + arg_images.isSet();
760 if (mutex > 1)
761 {
762 ERR("Only one output format can be specified (single-file, multi-file, "
763 "or images)");
764#ifdef USE_PETSC
765 MPI_Finalize();
766#endif
767 return EXIT_FAILURE;
768 }
769
770 std::cout << "OpenGeoSys NetCDF Converter\n";
771 if (!arg_varname.isSet())
772 {
773 std::cout << "File " << arg_input.getValue()
774 << " loaded. Press ENTER to display available data arrays.\n";
775 std::cin.ignore();
776 }
777
778 std::string const& output_name(arg_output.getValue());
779 std::string const& var_name = (arg_varname.isSet())
780 ? arg_varname.getValue()
781 : arraySelectionLoop(dataset);
782 NcVar const& var = dataset.getVar(var_name);
783 if (var.isNull())
784 {
785 ERR("Variable \"{:s}\" not found in file.", arg_varname.getValue());
786#ifdef USE_PETSC
787 MPI_Finalize();
788#endif
789 return EXIT_FAILURE;
790 }
791
792 std::vector<std::size_t> dim_idx_map(var.getDimCount(), 0);
793 bool is_time_dep(false);
794 if (arg_dim1.isSet() && arg_dim2.isSet())
795 {
796 is_time_dep = arg_dim_time.isSet();
797 if (!assignDimParams(var, dim_idx_map, arg_dim_time, arg_dim1, arg_dim2,
798 arg_dim3))
799 {
800#ifdef USE_PETSC
801 MPI_Finalize();
802#endif
803 return EXIT_FAILURE;
804 }
805 }
806 else
807 {
808 is_time_dep = dimensionSelectionLoop(var, dim_idx_map);
809 }
810
811 std::pair<std::size_t, std::size_t> time_bounds(0, 0);
812 if (is_time_dep)
813 time_bounds =
814 (arg_time_start.isSet())
815 ? assignTimeBounds(getDimVar(dataset, var, dim_idx_map[0]),
816 arg_time_start, arg_time_end)
817 : timestepSelectionLoop(var, dim_idx_map[0]);
818
820 if (arg_images.isSet())
821 {
822 output = OutputType::IMAGES;
823 }
824 else if (arg_multi_files.isSet())
825 {
826 output = OutputType::MULTIMESH;
827 }
828 else if (arg_single_file.isSet() || !is_time_dep ||
829 time_bounds.first == time_bounds.second)
830 {
831 output = OutputType::SINGLEMESH;
832 }
833 else
834 {
835 output = multFilesSelectionLoop(time_bounds);
836 }
837
838 std::size_t const temp_offset = (is_time_dep) ? 1 : 0;
839 std::size_t const n_dims = (var.getDimCount());
841
842 if (output != OutputType::IMAGES)
843 {
844 elem_type = (arg_elem_type.isSet())
845 ? assignElemType(arg_elem_type)
846 : elemSelectionLoop(n_dims - temp_offset);
847 if (elem_type == MeshLib::MeshElemType::INVALID)
848 elemSelectionLoop(n_dims - temp_offset);
849 }
850
851 if (arg_nodata.isSet())
852 {
853 no_data_input = arg_nodata.getValue();
854 }
855
856 if (!convert(dataset, var, output_name, dim_idx_map, is_time_dep,
857 time_bounds, output, elem_type))
858 {
859#ifdef USE_PETSC
860 MPI_Finalize();
861#endif
862 return EXIT_FAILURE;
863 }
864
865 std::cout << "Conversion finished successfully.\n";
866#ifdef USE_PETSC
867 MPI_Finalize();
868#endif
869 return EXIT_SUCCESS;
870}
Definition of the AsciiRasterInterface class.
#define OGS_FATAL(...)
Definition Error.h:26
Filename manipulation routines.
Git information.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Definition of the Mesh class.
static void showArrays(NcFile const &dataset)
static double no_data_input
static OutputType multFilesSelectionLoop(std::pair< std::size_t, std::size_t > const &time_bounds)
int main(int argc, char *argv[])
static std::string arraySelectionLoop(NcFile const &dataset)
static MathLib::Point3d getOrigin(NcFile const &dataset, NcVar const &var, std::vector< std::size_t > const &dim_idx_map, std::size_t const time_offset)
static GeoLib::RasterHeader createRasterHeader(NcFile const &dataset, NcVar const &var, std::vector< std::size_t > const &dim_idx_map, std::vector< std::size_t > const &length, std::size_t const time_offset)
static void showErrorMessage(std::size_t const error_id, std::size_t const max=0)
static std::pair< double, double > getDimLength(NcVar const &var, std::size_t const dim)
static std::pair< std::size_t, std::size_t > assignTimeBounds(NcVar const &var, TCLAP::ValueArg< std::size_t > &arg_time_start, TCLAP::ValueArg< std::size_t > &arg_time_end)
static std::pair< double, double > getBoundaries(NcVar const &var)
static std::size_t parseInput(std::string const &request_str, std::size_t const max_val, bool const has_info=false)
static std::string getIterationString(std::size_t i, std::size_t max)
static MeshLib::MeshElemType assignElemType(TCLAP::ValueArg< std::string > &arg_elem_type)
static MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim)
static bool assignDimParams(NcVar const &var, std::vector< std::size_t > &dim_idx_map, TCLAP::ValueArg< std::size_t > &arg_dim_time, TCLAP::ValueArg< std::size_t > &arg_dim1, TCLAP::ValueArg< std::size_t > &arg_dim2, TCLAP::ValueArg< std::size_t > &arg_dim3)
static std::pair< std::size_t, std::size_t > timestepSelectionLoop(NcVar const &var, std::size_t const dim_idx)
static const double no_data_output
static void checkExit(std::string const &str)
static void flipRaster(std::vector< double > &data, std::size_t const layers, std::size_t const width, std::size_t const height)
static std::vector< std::size_t > getLength(NcVar const &var, std::size_t const time_offset)
static std::vector< std::string > getArrays(NcFile const &dataset)
static NcVar getDimVar(NcFile const &dataset, NcVar const &var, std::size_t const dim)
static bool dimensionSelectionLoop(NcVar const &var, std::vector< std::size_t > &dim_idx_map)
static double getResolution(NcFile const &dataset, NcVar const &var)
static std::vector< double > getData(NcVar const &var, std::size_t const total_length, std::size_t const time_step, std::vector< std::size_t > const &length)
static bool convert(NcFile const &dataset, NcVar const &var, std::string const &output_name, std::vector< std::size_t > const &dim_idx_map, std::size_t const time_offset, std::pair< std::size_t, std::size_t > const &time_bounds, OutputType const output, MeshLib::MeshElemType const elem_type)
static bool canConvert(NcVar const &var)
static void showArraysDims(NcVar const &var)
Definition of the GeoLib::Raster class.
Implementation of the VtuInterface class.
static void writeRasterAsASC(GeoLib::Raster const &raster, std::string const &file_name)
Writes an Esri asc-file.
Class Raster is used for managing raster data.
Definition Raster.h:49
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
bool writeToFile(std::filesystem::path const &file_path)
static std::unique_ptr< MeshLib::Mesh > convert(GeoLib::Raster const &raster, MeshLib::MeshElemType elem_type, MeshLib::UseIntensityAs intensity_type, std::string const &array_name="Colour")
std::string getFileExtension(const std::string &path)
std::string dropFileExtension(std::string const &filename)
GITINFOLIB_EXPORT const std::string ogs_version
const Point3d ORIGIN
Definition Point3d.h:98
UseIntensityAs
Selection of possible interpretations for intensities.
Definition MeshEnums.h:83
MeshElemType
Types of mesh elements supported by OpenGeoSys. Values are from VTKCellType enum.
Definition MeshEnums.h:27
Contains the relevant information when storing a geoscientific raster data.
Definition Raster.h:28
std::size_t n_depth
Definition Raster.h:31
std::size_t n_cols
Definition Raster.h:29
std::size_t n_rows
Definition Raster.h:30