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