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