11 #include <tclap/CmdLine.h>
32 using namespace netCDF;
47 if (str ==
"x" || str ==
"exit")
52 std::size_t
const max = 0)
56 ERR(
"Input not valid.");
58 else if (error_id == 1)
60 ERR(
"Index not valid. Valid indices are in [0,{:d}].", max);
62 else if (error_id == 2)
64 ERR(
"Input not valid.");
65 std::cout <<
"Type \"info\" to display the available options again. "
66 "\"exit\" will exit the programme.\n";
70 static std::size_t
parseInput(std::string
const& request_str,
71 std::size_t
const max_val,
72 bool const has_info =
false)
76 std::cout << request_str;
79 std::getline(std::cin, str);
81 if (has_info && str ==
"info")
83 std::stringstream str_stream(str);
84 if (!(str_stream >> val))
86 std::size_t
const error_val = (has_info) ? 2 : 0;
90 if (val > max_val - 1)
97 return std::numeric_limits<std::size_t>::max();
100 static NcVar
getDimVar(NcFile
const& dataset, NcVar
const& var,
101 std::size_t
const dim)
103 NcDim
const& dim_obj = var.getDim(dim);
104 return dataset.getVar(dim_obj.getName());
108 std::size_t
const dim)
110 return std::make_pair(0.0,
static_cast<double>(var.getDim(dim).getSize()));
113 static std::vector<std::string>
getArrays(NcFile
const& dataset)
115 auto const& names = dataset.getVars();
116 std::vector<std::string> var_names;
117 for (
auto [
name, var] : names)
120 var_names.push_back(
name);
127 std::size_t
const n_vars(dataset.getDimCount());
128 std::cout <<
"The NetCDF file contains the following " << n_vars
130 std::cout <<
"\tIndex\tArray Name\t#Dimensions\n";
131 std::cout <<
"-------------------------------------------\n";
132 auto const& names = dataset.getVars();
133 std::size_t count = 0;
134 for (
auto [
name, var] : names)
136 std::cout <<
"\t" << count++ <<
"\t" <<
name <<
"\t("
137 << var.getDimCount() <<
"D array)\n";
144 std::cout <<
"Data array \"" << var.getName()
145 <<
"\" contains the following dimensions:\n";
146 std::size_t
const n_dims(var.getDimCount());
147 for (std::size_t i = 0; i < n_dims; ++i)
148 std::cout <<
"\t" << i <<
"\t" << var.getDim(i).getName() <<
"\t("
149 << var.getDim(i).getSize() <<
" values)\n";
155 if (var.getDimCount() == 1)
158 std::size_t
const size = var.getDim(0).getSize();
159 var.getVar({0}, {1}, &start);
160 var.getVar({size - 1}, {1}, &end);
161 return std::make_pair(start, end);
163 return std::make_pair(0, 0);
167 std::vector<std::size_t>
const& dim_idx_map,
168 std::size_t
const time_offset)
171 std::size_t
const n_dims = var.getDimCount();
172 for (std::size_t i = time_offset; i < n_dims; ++i)
174 NcVar
const& dim =
getDimVar(dataset, var, dim_idx_map[i]);
175 auto const bounds = (dim.isNull()) ?
getDimLength(var, dim_idx_map[i])
177 origin[i - time_offset] =
178 (bounds.first < bounds.second) ? bounds.first : bounds.second;
183 static void flipRaster(std::vector<double>& data, std::size_t
const layers,
184 std::size_t
const width, std::size_t
const height)
186 std::size_t
const length(data.size());
187 std::vector<double> tmp_vec;
188 tmp_vec.reserve(length);
189 for (std::size_t k = 0; k < layers; k++)
191 std::size_t
const layer_end = (k + 1) * height * width;
192 for (std::size_t i = 0; i < height; i++)
194 std::size_t
const line_idx(layer_end - (width * (i + 1)));
195 for (std::size_t j = 0; j < width; j++)
197 tmp_vec.push_back(data[line_idx + j]);
201 std::copy(tmp_vec.cbegin(), tmp_vec.cend(), data.begin());
206 bool ret(var.getDimCount() < 2);
208 ERR(
"Only 2+ dimensional variables can be converted into OGS "
215 std::vector<std::string>
const& names =
getArrays(dataset);
217 std::size_t
const idx =
218 parseInput(
"Enter data array index: ", dataset.getVarCount(),
true);
220 if (
static_cast<int>(idx) == dataset.getVarCount() ||
228 std::vector<std::size_t>& dim_idx_map)
231 std::size_t
const n_dims(var.getDimCount());
232 dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
233 bool is_time_dep(
true);
238 std::string temp_str(
"");
239 cout <<
"Is the parameter time-dependent?\n";
240 while (dim_idx_map[0] == std::numeric_limits<std::size_t>::max() &&
243 cout <<
"Enter ID for temporal dimension or \"c\" to continue: ";
244 std::getline(std::cin, temp_str);
245 std::stringstream str_stream(temp_str);
246 if (str_stream.str() ==
"c" || str_stream.str() ==
"continue")
250 if (!(str_stream >> dim_idx_map[0]))
253 dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
256 if (dim_idx_map[0] > n_dims - 1)
259 dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
268 std::size_t
const start_idx = (is_time_dep) ? 1 : 0;
269 std::array<std::string, 4>
const dim_comment{
270 "(x / longitude)",
"(y / latitude)",
"(z / height / depth)",
271 "[Error: 4-dimensional non-temporal arrays are not supported]"};
272 for (std::size_t i = start_idx; i < n_dims; ++i)
274 dim_idx_map[i] = std::numeric_limits<std::size_t>::max();
276 std::string
const request_str(
"Enter ID for dimension " +
277 std::to_string(i) +
" " +
278 dim_comment[i - start_idx] +
": ");
279 std::size_t
const idx =
280 parseInput(request_str, var.getDimCount(),
true);
282 if (
static_cast<int>(idx) == var.getDimCount())
288 dim_idx_map[i] = idx;
295 NcVar
const& var, std::size_t
const dim_idx)
297 std::size_t
const n_time_steps = var.getDim(dim_idx).getSize();
298 std::size_t
const max_val = std::numeric_limits<std::size_t>::max();
299 std::pair<std::size_t, std::size_t> bounds(max_val, max_val);
300 std::cout <<
"\nThe dataset contains " << n_time_steps <<
" time steps.\n";
301 while (bounds.first == max_val)
304 "Specify first time step to export: ", n_time_steps,
false);
306 while (bounds.first > bounds.second || bounds.second > n_time_steps)
309 "Specify last time step to export: ", n_time_steps,
false);
322 std::cout <<
"\nSelect element type for result, choose ";
325 std::cout <<
"(t)riangle or (q)uadliteral: ";
327 std::cout <<
"(p)rism or (h)exahedron: ";
328 std::string type(
"");
329 std::getline(std::cin, type);
333 if (type !=
"t" && type !=
"q" && type !=
"tri" && type !=
"quad" &&
334 type !=
"triangle" && type !=
"quatliteral")
336 if (type ==
"t" || type ==
"tri" || type ==
"triangle")
343 if (type !=
"p" && type !=
"h" && type !=
"prism" &&
344 type !=
"hex" && type !=
"hexahedron")
346 if (type ==
"p" || type ==
"prism")
355 std::pair<std::size_t, std::size_t>
const& time_bounds)
360 std::size_t
const n_time_steps(time_bounds.second - time_bounds.first +
362 std::cout <<
"\nThe selection includes " << n_time_steps
364 std::cout <<
"0. Save data in " << n_time_steps
365 <<
" mesh files with one scalar array each.\n";
366 std::cout <<
"1. Save data in one mesh file with " << n_time_steps
367 <<
" scalar arrays.\n";
368 std::cout <<
"2. Save data as " << n_time_steps <<
" ASC images.\n";
370 std::size_t
const ret =
371 parseInput(
"Select preferred method: ", 3,
false);
385 std::size_t
const max_length(std::to_string(max).length());
386 std::string
const current_str(std::to_string(i));
387 return std::string(max_length - current_str.length(),
'0') + current_str;
392 std::size_t
const dim_idx = var.getDimCount() - 1;
393 NcVar
const dim_var(
getDimVar(dataset, var, dim_idx));
394 auto const bounds = (dim_var.isNull()) ?
getDimLength(var, dim_idx)
396 std::size_t
const dim_size = var.getDim(dim_idx).getSize();
399 OGS_FATAL(
"Dimension '{:s}' has size 0. Aborting...",
400 var.getDim(dim_idx).getName());
402 return std::fabs(bounds.second - bounds.first) /
403 static_cast<double>(dim_size);
407 NcFile
const& dataset, NcVar
const& var,
408 std::vector<std::size_t>
const& dim_idx_map,
409 std::vector<std::size_t>
const& length, std::size_t
const time_offset)
412 getOrigin(dataset, var, dim_idx_map, time_offset);
414 std::size_t n_dims = var.getDimCount();
415 std::size_t z_length =
416 (n_dims - time_offset == 3) ? length[dim_idx_map.back()] : 1;
417 return {length[dim_idx_map[0 + time_offset]],
418 length[dim_idx_map[1 + time_offset]],
425 static std::vector<std::size_t>
getLength(NcVar
const& var,
426 std::size_t
const time_offset)
428 std::size_t
const n_dims = (var.getDimCount());
429 std::vector<std::size_t> length(n_dims, 1);
430 for (std::size_t i = time_offset; i < n_dims; ++i)
432 length[i] = var.getDim(i).getSize();
437 static std::vector<double>
getData(NcVar
const& var,
438 std::size_t
const total_length,
439 std::size_t
const time_step,
440 std::vector<std::size_t>
const& length)
442 std::size_t
const n_dims(var.getDimCount());
443 std::vector<std::size_t> offset(n_dims, 0);
444 offset[0] = time_step;
445 std::vector<double> data_vec(total_length, 0);
446 var.getVar(offset, length, data_vec.data());
449 data_vec.begin(), data_vec.end(),
450 [](
double const& x) { return x == no_data_input; },
no_data_output);
456 std::vector<std::size_t>& dim_idx_map,
457 TCLAP::ValueArg<std::size_t>& arg_dim_time,
458 TCLAP::ValueArg<std::size_t>& arg_dim1,
459 TCLAP::ValueArg<std::size_t>& arg_dim2,
460 TCLAP::ValueArg<std::size_t>& arg_dim3)
462 std::size_t dim_param_count = 0;
463 if (arg_dim_time.isSet())
465 if (arg_dim1.isSet())
467 if (arg_dim2.isSet())
469 if (arg_dim3.isSet())
472 std::size_t
const n_dims = var.getDimCount();
473 if (dim_param_count != n_dims)
475 ERR(
"Number of parameters set does not fit number of parameters for "
476 "specified variable.");
480 if (arg_dim_time.getValue() >= n_dims || arg_dim1.getValue() >= n_dims ||
481 arg_dim2.getValue() >= n_dims || arg_dim3.getValue() >= n_dims)
483 ERR(
"Maximum allowed dimension for variable \"{:s}\" is {:d}.",
484 var.getName(), n_dims - 1);
488 if (arg_dim_time.isSet())
489 dim_idx_map[0] = arg_dim_time.getValue();
490 std::size_t
const temp_offset = (arg_dim_time.isSet()) ? 1 : 0;
491 dim_idx_map[0 + temp_offset] = arg_dim1.getValue();
492 dim_idx_map[1 + temp_offset] = arg_dim2.getValue();
493 if (n_dims == (3 + temp_offset))
494 dim_idx_map[2 + temp_offset] = arg_dim3.getValue();
501 TCLAP::ValueArg<std::size_t>& arg_time_start,
502 TCLAP::ValueArg<std::size_t>& arg_time_end)
505 if (arg_time_start.getValue() > bounds.second)
507 ERR(
"Start time step larger than total number of time steps. Resetting "
509 arg_time_start.reset();
512 if (!arg_time_end.isSet())
513 return {arg_time_start.getValue(), arg_time_start.getValue()};
515 if (arg_time_end.getValue() > bounds.second)
517 ERR(
"End time step larger than total number of time steps. Resetting "
518 "to starting time step");
519 return {arg_time_start.getValue(), arg_time_start.getValue()};
522 if (arg_time_end.getValue() < arg_time_start.getValue())
524 ERR(
"End time step larger than starting time step. Swapping values");
525 return {arg_time_end.getValue(), arg_time_start.getValue()};
528 return {arg_time_start.getValue(), arg_time_end.getValue()};
532 TCLAP::ValueArg<std::string>& arg_elem_type)
534 if (arg_elem_type.getValue() ==
"tri")
536 if (arg_elem_type.getValue() ==
"quad")
538 if (arg_elem_type.getValue() ==
"prism")
540 if (arg_elem_type.getValue() ==
"hex")
546 static bool convert(NcFile
const& dataset, NcVar
const& var,
547 std::string
const& output_name,
548 std::vector<std::size_t>
const& dim_idx_map,
549 std::size_t
const time_offset,
550 std::pair<std::size_t, std::size_t>
const& time_bounds,
554 std::unique_ptr<MeshLib::Mesh> mesh;
555 std::vector<std::size_t>
const length =
getLength(var, time_offset);
556 std::size_t
const array_length = std::accumulate(
557 length.cbegin(), length.cend(), 1, std::multiplies<std::size_t>());
558 for (std::size_t i = time_bounds.first; i <= time_bounds.second; ++i)
560 std::string
const step_str =
561 (time_bounds.first != time_bounds.second)
562 ? std::string(
" time step " + std::to_string(i))
564 std::cout <<
"Converting" << step_str <<
"...\n";
565 std::vector<double> data_vec =
getData(var, array_length, i, length);
569 std::size_t
const n_dims = length.size();
570 NcVar
const dim_var(
getDimVar(dataset, var, n_dims - 2));
571 auto const bounds = (dim_var.isNull()) ?
getDimLength(var, n_dims - 2)
573 if (bounds.first > bounds.second)
575 std::size_t n_layers =
576 (length.size() - time_offset == 3) ? length[n_dims - 3] : 1;
577 flipRaster(data_vec, n_layers, length[n_dims - 1],
588 elem_type, useIntensity,
590 std::string
const output_file_name(
598 std::string array_name(var.getName());
599 if (time_bounds.first != time_bounds.second)
601 if (i == time_bounds.first)
603 elem_type, useIntensity,
607 std::unique_ptr<MeshLib::Mesh>
const temp(
609 elem_type, useIntensity,
612 temp->getProperties().getPropertyVector<
double>(array_name);
615 MeshLib::addPropertyToMesh<double>(
618 if (i == time_bounds.second)
621 std::string
const output_file_name =
624 : output_name +
".vtu";
631 data_vec.data() + header.
n_cols *
634 std::string
const output_file_name(
644 int main(
int argc,
char* argv[])
647 "Converts NetCDF data into mesh file(s).\n\n "
648 "OpenGeoSys-6 software, version " +
651 "Copyright (c) 2012-2021, OpenGeoSys Community "
652 "(http://www.opengeosys.org)",
655 TCLAP::ValueArg<int> arg_nodata(
657 "explicitly specifies the no data value used in the dataset (usually "
658 "it is not necessary to set this)",
662 std::vector<std::string> allowed_elems{
"tri",
"quad",
"prism",
"hex"};
663 TCLAP::ValuesConstraint<std::string> allowed_elem_vals(allowed_elems);
664 TCLAP::ValueArg<std::string> arg_elem_type(
665 "e",
"elem-type",
"the element type used in the resulting OGS mesh",
666 false,
"", &allowed_elem_vals);
667 cmd.add(arg_elem_type);
669 TCLAP::SwitchArg arg_images(
671 "if set, all time steps will be written as ESRI image files (*.asc)");
674 TCLAP::SwitchArg arg_multi_files(
676 "if set, each time step will be written to a separate mesh file");
677 cmd.add(arg_multi_files);
679 TCLAP::SwitchArg arg_single_file(
681 "if set, all time steps will be written to a single mesh file (with "
682 "one scalar array per time step)");
683 cmd.add(arg_single_file);
685 TCLAP::ValueArg<std::size_t> arg_time_end(
687 "last time step to be extracted (only for time-dependent variables!)",
688 false, 0,
"integer specifying index of time step");
689 cmd.add(arg_time_end);
691 TCLAP::ValueArg<std::size_t> arg_time_start(
692 "",
"timestep-first",
693 "first time step to be extracted (only for time-dependent variables!)",
694 false, 0,
"integer specifying index of time step");
695 cmd.add(arg_time_start);
697 std::vector<std::size_t> allowed_dims{0, 1, 2, 3};
698 TCLAP::ValuesConstraint<std::size_t> allowed_dim_vals(allowed_dims);
699 TCLAP::ValueArg<std::size_t> arg_dim3(
701 "index of third dimension (z/height/depth) for the selected variable",
702 false, 0, &allowed_dim_vals);
705 TCLAP::ValueArg<std::size_t> arg_dim2(
707 "index of second dimension (y/latitude) for the selected variable",
708 false, 0, &allowed_dim_vals);
711 TCLAP::ValueArg<std::size_t> arg_dim1(
713 "index of first dimension (x/longitude) for the selected variable",
714 false, 0, &allowed_dim_vals);
717 TCLAP::ValueArg<std::size_t> arg_dim_time(
719 "index of the time-dependent dimension for the selected variable",
720 false, 0, &allowed_dim_vals);
721 cmd.add(arg_dim_time);
723 TCLAP::ValueArg<std::string> arg_varname(
724 "v",
"var",
"variable included in the the netCDF file",
false,
"",
725 "string containing the variable name");
726 cmd.add(arg_varname);
728 TCLAP::ValueArg<std::string> arg_output(
729 "o",
"output",
"the OGS mesh output file",
true,
"",
730 "string containing the path and file name");
732 TCLAP::ValueArg<std::string> arg_input(
733 "i",
"input",
"the netCDF input file",
true,
"",
734 "string containing the path and file name");
736 cmd.parse(argc, argv);
738 NcFile dataset(arg_input.getValue().c_str(), NcFile::read);
740 if (dataset.isNull())
742 ERR(
"Error opening file.");
746 std::size_t
const mutex =
747 arg_single_file.isSet() + arg_multi_files.isSet() + arg_images.isSet();
750 ERR(
"Only one output format can be specified (single-file, multi-file, "
755 std::cout <<
"OpenGeoSys NetCDF Converter\n";
756 if (!arg_varname.isSet())
758 std::cout <<
"File " << arg_input.getValue()
759 <<
" loaded. Press ENTER to display available data arrays.\n";
763 std::string
const& output_name(arg_output.getValue());
764 std::string
const& var_name = (arg_varname.isSet())
765 ? arg_varname.getValue()
767 NcVar
const& var = dataset.getVar(var_name);
770 ERR(
"Variable \"{:s}\" not found in file.", arg_varname.getValue());
774 std::vector<std::size_t> dim_idx_map(var.getDimCount(), 0);
775 bool is_time_dep(
false);
776 if (arg_dim1.isSet() && arg_dim2.isSet())
778 is_time_dep = arg_dim_time.isSet();
779 if (!
assignDimParams(var, dim_idx_map, arg_dim_time, arg_dim1, arg_dim2,
788 std::pair<std::size_t, std::size_t> time_bounds(0, 0);
791 (arg_time_start.isSet())
793 arg_time_start, arg_time_end)
797 if (arg_images.isSet())
801 else if (arg_multi_files.isSet())
805 else if (arg_single_file.isSet() || !is_time_dep ||
806 time_bounds.first == time_bounds.second)
815 std::size_t
const temp_offset = (is_time_dep) ? 1 : 0;
816 std::size_t
const n_dims = (var.getDimCount());
821 elem_type = (arg_elem_type.isSet())
828 if (arg_nodata.isSet())
833 if (!
convert(dataset, var, output_name, dim_idx_map, is_time_dep,
834 time_bounds, output, elem_type))
837 std::cout <<
"Conversion finished successfully.\n";
Definition of the AsciiRasterInterface class.
void ERR(char const *fmt, Args const &... args)
Definition of the Mesh class.
static void showArrays(NcFile const &dataset)
static std::pair< std::size_t, std::size_t > timestepSelectionLoop(NcVar const &var, std::size_t const dim_idx)
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 std::pair< double, double > getBoundaries(NcVar const &var)
static void showErrorMessage(std::size_t const error_id, std::size_t const max=0)
static std::vector< std::size_t > getLength(NcVar const &var, std::size_t const time_offset)
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::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 std::vector< std::string > getArrays(NcFile const &dataset)
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 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 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 std::pair< double, double > getDimLength(NcVar const &var, std::size_t const dim)
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.
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, MeshElemType elem_type, 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
void copy(PETScVector const &x, PETScVector &y)
UseIntensityAs
Selection of possible interpretations for intensities.
MeshElemType
Types of mesh elements supported by OpenGeoSys. Values are from VTKCellType enum.