diff --git a/Applications/Utils/FileConverter/CMakeLists.txt b/Applications/Utils/FileConverter/CMakeLists.txt index 761d8a6e8f8352e372aad047e94b1b16860a5298..230b4833461f110c257147a28a47ccfaeb57ea2f 100644 --- a/Applications/Utils/FileConverter/CMakeLists.txt +++ b/Applications/Utils/FileConverter/CMakeLists.txt @@ -33,7 +33,43 @@ if(TARGET ConvertSHPToGLI) target_link_libraries(ConvertSHPToGLI GeoLib Qt5::Xml ${Shapelib_LIBRARIES}) endif() -if(TARGET Mesh2Shape) - target_link_libraries(Mesh2Shape ${Shapelib_LIBRARIES}) +add_executable(OGS2VTK OGS2VTK.cpp) +set_target_properties(OGS2VTK PROPERTIES FOLDER Utilities) +target_link_libraries(OGS2VTK MeshLib) +ADD_VTK_DEPENDENCY(OGS2VTK) + +add_executable(VTK2OGS VTK2OGS.cpp) +set_target_properties(VTK2OGS PROPERTIES FOLDER Utilities) +target_link_libraries(VTK2OGS MeshLib) +ADD_VTK_DEPENDENCY(VTK2OGS) + +add_executable(VTK2TIN VTK2TIN.cpp) +set_target_properties(VTK2TIN PROPERTIES FOLDER Utilities) +target_link_libraries(VTK2TIN MeshLib) +ADD_VTK_DEPENDENCY(VTK2TIN) + +add_executable(TIN2VTK TIN2VTK.cpp) +set_target_properties(TIN2VTK PROPERTIES FOLDER Utilities) +target_link_libraries(TIN2VTK MeshLib) +ADD_VTK_DEPENDENCY(TIN2VTK) + +add_executable(NetCdfConverter NetCdfConverter.cpp) +set_target_properties(NetCdfConverter PROPERTIES FOLDER Utilities) +target_link_libraries(NetCdfConverter MeshLib) +ADD_VTK_DEPENDENCY(NetCdfConverter) + +#################### +### Installation ### +#################### +install(TARGETS generateMatPropsFromMatID GMSH2OGS OGS2VTK VTK2OGS VTK2TIN + RUNTIME DESTINATION bin COMPONENT ogs_converter) + +if(Qt5XmlPatterns_FOUND) + if(Shapelib_FOUND) + install(TARGETS ConvertSHPToGLI + RUNTIME DESTINATION bin COMPONENT ogs_converter) + endif() + install(TARGETS FEFLOW2OGS convertGEO + RUNTIME DESTINATION bin COMPONENT ogs_converter) endif() diff --git a/Applications/Utils/FileConverter/NetCdfConverter.cpp b/Applications/Utils/FileConverter/NetCdfConverter.cpp index abd41633a53a8daf4857ed070b013af27896bff7..f68447b56ebbbb23a9a56aca5006db879f8a1bdf 100644 --- a/Applications/Utils/FileConverter/NetCdfConverter.cpp +++ b/Applications/Utils/FileConverter/NetCdfConverter.cpp @@ -24,6 +24,7 @@ // BaseLib #include "BaseLib/BuildInfo.h" #include "BaseLib/LogogSimpleFormatter.h" +#include "BaseLib/FileTools.h" // GeoLib #include "GeoLib/Raster.h" @@ -35,17 +36,55 @@ #include "MeshLib/IO/VtkIO/VtuInterface.h" +void checkExit(std::string const& str) +{ + if (str == "exit") + exit(0); +} + void showErrorMessage(std::size_t error_id, std::size_t max = 0) { if (error_id == 0) { - ERR ("Index not recognised."); - std::cout << "\"info\" will display the available options again. \"exit\" will exit the programme.\n"; + ERR ("Input not valid."); } else if (error_id == 1) { ERR ("Index not valid. Valid indices are in [0,%d].", max); } + else if (error_id == 2) + { + showErrorMessage(0); + std::cout << "Type \"info\" to display the available options again. \"exit\" will exit the programme.\n"; + } +} + +std::size_t parseInput(std::string const& request_str, std::size_t max_val, bool has_info = false) +{ + while (true) + { + std::cout << request_str; + std::string str; + std::size_t val; + std::getline(std::cin, str); + checkExit(str); + if (has_info && str == "info") + return (max_val); + std::stringstream str_stream(str); + if (!(str_stream >> val)) + { + std::size_t error_val = (has_info) ? 2 : 0; + showErrorMessage(error_val); + continue; + } + if (val > max_val - 1) + { + showErrorMessage(1, max_val - 1); + continue; + } + return val; + } + return std::numeric_limits<std::size_t>::max(); } void showArrays(NcFile const& dataset) @@ -62,7 +101,7 @@ void showArrays(NcFile const& dataset) std::cout << "\n"; } -bool showArraysDims(NcFile const& dataset, NcVar const& var) +bool showArraysDims(NcVar const& var) { std::cout << "Data array \"" << var.name() << "\" contains the following dimensions:\n"; std::size_t const n_dims (var.num_dims()); @@ -74,7 +113,7 @@ bool showArraysDims(NcFile const& dataset, NcVar const& var) std::size_t getDimensionBoundaries(NcFile const& dataset, NcVar const& var, std::size_t dim, double &start, double &end) { - NcVar dim_var = *dataset.get_var(var.get_dim(dim)->name()); + NcVar& dim_var = *dataset.get_var(var.get_dim(dim)->name()); if ((dim_var.num_dims()) == 1) { std::size_t size = dim_var.get_dim(0)->size(); @@ -95,22 +134,32 @@ std::size_t getDimensionBoundaries(NcFile const& dataset, NcVar const& var, std: return 0; } +MathLib::Point3d getOrigin(std::array<std::pair<double, double>,3> bounds, std::size_t n_dims) +{ + MathLib::Point3d origin(std::array<double, 3>{ {0, 0, 0}}); + for (std::size_t i=0; i<n_dims; ++i) + origin[i] = (bounds[i].first < bounds[i].second) ? bounds[i].first : bounds[i].second; + if (n_dims > 1) + std::swap(origin[0], origin[1]); + return origin; +} + void reverseNorthSouth(double* data, std::size_t width, std::size_t height) { - double* cp_array = new double[width*height]; + std::size_t const total_length (width * height); + double* cp_array = new double[total_length]; for (std::size_t i=0; i<height; i++) { for (std::size_t j=0; j<width; j++) { - std::size_t const old_index ((width*height)-(width*(i+1))); + std::size_t const old_index (total_length - (width*(i+1))); std::size_t const new_index (width*i); cp_array[new_index+j] = data[old_index+j]; } } - std::size_t const length(height*width); - for (std::size_t i=0; i<length; i++) + for (std::size_t i=0; i<total_length; i++) data[i] = cp_array[i]; delete[] cp_array; @@ -119,39 +168,17 @@ void reverseNorthSouth(double* data, std::size_t width, std::size_t height) std::size_t arraySelectionLoop(NcFile const& dataset) { showArrays(dataset); - while (true) - { - cout << "Enter data array index: "; - std::string var_idx_str (""); - std::getline(std::cin, var_idx_str); - if (var_idx_str == "info") - { - showArrays(dataset); - continue; - } - if (var_idx_str == "exit") - exit(0); - - std::stringstream str_stream(var_idx_str); - std::size_t var_idx; - if (!(str_stream >> var_idx)) - { - showErrorMessage(0); - continue; - } - if (var_idx > dataset.num_vars()-1) - { - showErrorMessage(1, dataset.num_vars()-1); - continue; - } - return var_idx; - } - return std::numeric_limits<std::size_t>::max(); + std::size_t idx = parseInput("Enter data array index: ", dataset.num_vars(), true); + + if (idx == dataset.num_vars()) + return arraySelectionLoop(dataset); + + return idx; } -bool dimensionSelectionLoop(NcFile const& dataset, NcVar const& var, std::array<std::size_t,4> &dim_idx_map) +bool dimensionSelectionLoop(NcVar const& var, std::array<std::size_t,4> &dim_idx_map) { - showArraysDims(dataset, var); + showArraysDims(var); std::size_t const n_dims (var.num_dims()); dim_idx_map[0] = std::numeric_limits<std::size_t>::max(); bool is_time_dep (true); @@ -166,59 +193,60 @@ bool dimensionSelectionLoop(NcFile const& dataset, NcVar const& var, std::array< cout << "Enter ID for temporal dimension or \"c\" to continue: "; std::getline(std::cin, temp_str); std::stringstream str_stream(temp_str); - if (str_stream.str() == "c") + if (str_stream.str() == "c" || str_stream.str() == "continue") is_time_dep = false; else { - if (str_stream >> dim_idx_map[0]) + if (!(str_stream >> dim_idx_map[0])) + { + showErrorMessage(0); + continue; + } + if (dim_idx_map[0] > n_dims-1) { - if (dim_idx_map[0] > n_dims-1) - { - showErrorMessage(1, var.num_dims()-1); - dim_idx_map[0] = std::numeric_limits<std::size_t>::max(); - } + showErrorMessage(1, var.num_dims()-1); + dim_idx_map[0] = std::numeric_limits<std::size_t>::max(); } } } } + else + is_time_dep = false; // get spatial dimension(s) std::size_t const start_idx = (is_time_dep) ? 1 : 0; for (std::size_t i=start_idx; i<n_dims; ++i) { dim_idx_map[i] = std::numeric_limits<std::size_t>::max(); - while (dim_idx_map[i] == std::numeric_limits<std::size_t>::max()) - { - std::string dim_idx_str; - cout << "Enter ID for dimension " << i << ": "; - std::getline(std::cin, dim_idx_str); - if (dim_idx_str == "info") - { - showArraysDims(dataset, var); - continue; - } - if (dim_idx_str == "exit") - exit(0); - std::stringstream str_stream(dim_idx_str); - std::size_t dim_idx; - if (!(str_stream >> dim_idx)) - { - showErrorMessage(0); - continue; - } - if (dim_idx > var.num_dims()-1) - { - showErrorMessage(1, var.num_dims()-1); - continue; - } - dim_idx_map[i] = dim_idx; + std::string request_str ("Enter ID for dimension " + std::to_string(i) + ": "); + std::size_t idx = parseInput(request_str, var.num_dims(), true); + + if (idx == var.num_dims()) + { + showArraysDims(var); + i--; + continue; } + dim_idx_map[i] = idx; } return is_time_dep; } +std::pair<std::size_t, std::size_t> timestepSelectionLoop(NcFile const& dataset, NcVar const& var, std::size_t dim) +{ + double start, end; + std::string str(""); + std::pair<std::size_t, std::size_t> bounds(0,std::numeric_limits<std::size_t>::max()); + std::size_t n_time_steps = getDimensionBoundaries(dataset, var, dim, start, end); + std::cout << "\nThe dataset contains " << n_time_steps << " time steps.\n"; + bounds.first = parseInput("Specify first time step to export: ", n_time_steps, false); + while (bounds.first > bounds.second || bounds.second == std::numeric_limits<std::size_t>::max()) + bounds.second = parseInput("Specify last time step to export: ", n_time_steps, false); + return bounds; +} + MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim) { if (dim ==1) @@ -233,12 +261,12 @@ MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim) if (dim==3) std::cout << "(p)rism or (h)exahedron: "; std::string type (""); std::getline(std::cin, type); - + checkExit(type); if (dim==2) { - if (type != "t" || type != "q" || - type != "tri" || type != "quad" || - type != "triangle" || type != "quatliteral") + if (type != "t" && type != "q" && + type != "tri" && type != "quad" && + type != "triangle" && type != "quatliteral") continue; if (type == "t" || type == "tri" || type == "triangle") return MeshLib::MeshElemType::TRIANGLE; @@ -261,7 +289,17 @@ MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim) return t; } -void setOrigin(NcVar &var, std::size_t time_step) +bool multFilesSelectionLoop(std::pair<std::size_t, std::size_t> const& time_bounds) +{ + std::size_t n_time_steps(time_bounds.second - time_bounds.first); + std::cout << "\nThe selection includes " << n_time_steps << " time steps.\n"; + std::cout << "0. Save data in " << n_time_steps << " mesh files with one scalar array each.\n"; + std::cout << "1. Save data in one mesh file with " << n_time_steps << " scalar arrays.\n"; + std::size_t ret = parseInput("Select preferred method: ", 2, false); + return (ret == 0) ? true : false; +} + +void setTimeStep(NcVar &var, std::size_t time_step) { long* newOrigin = new long[var.num_dims()]; for (int i=0; i < var.num_dims(); ++i) @@ -271,12 +309,13 @@ void setOrigin(NcVar &var, std::size_t time_step) delete [] newOrigin; } -/* -MeshLib::Mesh convert(double* data_array, std::array<std::size_t,3> length, GeoLib::Point &origin, double resolution, MeshLib::MeshElemType meshElemType, MeshLib::UseIntensityAs useIntensity, std::string &name) +std::string getIterationString(std::size_t i, std::size_t max) { - return + std::size_t const max_length (std::to_string(max).length()); + std::string const current_str (std::to_string(i)); + return std::string(max_length - current_str.length(), '0') + current_str; } -*/ + int main (int argc, char* argv[]) { ApplicationsLib::LogogSetup logog_setup; @@ -294,126 +333,137 @@ int main (int argc, char* argv[]) NcFile dataset(input.getValue().c_str(), NcFile::ReadOnly); + if (!dataset.is_valid()) + { + ERR("Error opening file."); + return -1; + } + std::cout << "OpenGeoSys NetCDF Converter\n"; std::cout << "File " << input.getValue() << " loaded. Press ENTER to display available data arrays.\n"; std::cin.ignore(); + std::string output_name (output.getValue()); std::size_t const var_idx = arraySelectionLoop(dataset); + NcVar& var = *dataset.get_var(var_idx); std::array<std::size_t, 4> dim_idx_map; - NcVar var = *dataset.get_var(var_idx); - bool const is_time_dep = dimensionSelectionLoop(dataset, var, dim_idx_map); - - std::cin.ignore(); + bool const is_time_dep = dimensionSelectionLoop(var, dim_idx_map); std::size_t temp_offset = (is_time_dep) ? 1 : 0; std::size_t const n_dims = (var.num_dims()); std::size_t* length = new std::size_t[n_dims]; for(int i=0; i < n_dims; i++) length[i]=1; - std::array<double, 3> spatial_size = {0,0,0}; - std::array<double, 3> start = {0,0,0}; - std::array<double, 3> end = {0,0,0}; + std::array<std::size_t, 3> spatial_size = {0,0,0}; + std::array<std::pair<double, double>,3> spatial_bounds; + for (std::size_t i=0; i<3; ++i) + spatial_bounds[i] = std::pair<double, double>(0, 0); std::size_t array_length = 1; for (std::size_t i=temp_offset; i<n_dims; ++i) { - length[i] = getDimensionBoundaries(dataset, var, dim_idx_map[i], start[i-temp_offset], end[i-temp_offset]); - spatial_size[i-temp_offset] = length[i]; - array_length *= length[i]; + length[dim_idx_map[i]] = getDimensionBoundaries(dataset, var, dim_idx_map[i], spatial_bounds[i-temp_offset].first, spatial_bounds[i-temp_offset].second); + spatial_size[i-temp_offset] = length[dim_idx_map[i]]; + array_length *= length[dim_idx_map[i]]; } - setOrigin(var, 0); - - double* data_array = new double[array_length]; - for (std::size_t i=0; i < array_length; i++) - data_array[i]=0; - var.get(data_array, length); //create Array of Values - - for (std::size_t i=0; i < array_length; i++) - if (data_array[i] < -9999 ) data_array[i] = -9999; // all values < -10000, set to "no-value" - - double origin_x = (start[1] < end[1]) ? start[1] : end[1]; - double origin_y = (start[0] < end[0]) ? start[0] : end[0]; - MathLib::Point3d origin (std::array<double,3>{{origin_x, origin_y, 0}}); - double const resolution = fabs(end[0]-start[0])/(spatial_size[0]-1); - - // reverse lines in vertical direction if the original file has its origin in the northwest corner - if (start[0] > end[0]) - reverseNorthSouth(data_array, spatial_size[1], spatial_size[0]); - - GeoLib::RasterHeader header = {spatial_size[1], spatial_size[0], origin, resolution, -9999}; - - //MeshLib::MeshElemType meshElemType = elemSelectionLoop(n_dims-temp_offset); - - MeshLib::MeshElemType meshElemType = MeshLib::MeshElemType::TRIANGLE; + MathLib::Point3d origin = getOrigin(spatial_bounds, n_dims); + double const resolution = fabs(spatial_bounds[0].second - spatial_bounds[0].first) / (spatial_size[0] - 1); + GeoLib::RasterHeader header = { spatial_size[1], spatial_size[0], origin, resolution, -9999 }; + MeshLib::MeshElemType meshElemType = elemSelectionLoop(n_dims - temp_offset); MeshLib::UseIntensityAs useIntensity = MeshLib::UseIntensityAs::DATAVECTOR; - MeshLib::Mesh* mesh = MeshLib::RasterToMesh::convert(data_array, header, meshElemType, useIntensity, var.name()); + std::pair<std::size_t, std::size_t> time_bounds(0,0); + if (is_time_dep) + time_bounds = timestepSelectionLoop(dataset, var, dim_idx_map[0]); - /*********************** - std::size_t const n_time_steps = getDimensionBoundaries(dataset, var, dim_idx_map[2], start_lat, end_lat); - std::size_t const array_length (sizeLat * sizeLon); - for (std::size_t j=1; j<n_time_steps; ++j) + bool mult_files (false); + if (is_time_dep && time_bounds.first != time_bounds.second) + mult_files = multFilesSelectionLoop(time_bounds); + + double* data_array = nullptr; + std::unique_ptr<MeshLib::Mesh> mesh; + for (std::size_t i=time_bounds.first; i<=time_bounds.second; ++i) { - setOrigin(var, j); - var.get(data_array, length); //create Array of Values + if (is_time_dep) + setTimeStep(var, i); - for (std::size_t i=0; i < (sizeLat*sizeLon); i++) - if (data_array[i] < -9999 ) data_array[i] = -9999; // all values < -10000, set to "no-value" + data_array = new double[array_length]; + for (std::size_t i = 0; i < array_length; i++) + data_array[i] = 0; + var.get(data_array, length); - double const resolution = fabs(end_lat-start_lat)/(sizeLat-1); + for (std::size_t i=0; i < array_length; i++) + if (data_array[i] < -9999 ) data_array[i] = -9999; // all values < -10000, set to "no-value" // reverse lines in vertical direction if the original file has its origin in the northwest corner - if (start_lat > end_lat) - reverseNorthSouth(data_array, sizeLon, sizeLat); - GeoLib::RasterHeader header = {sizeLon, sizeLat, origin, resolution, -9999}; + if (spatial_bounds[0].first > spatial_bounds[0].second) + reverseNorthSouth(data_array, spatial_size[1], spatial_size[0]); - MeshLib::MeshElemType meshElemType = MeshLib::MeshElemType::TRIANGLE; - MeshLib::UseIntensityAs useIntensity = MeshLib::UseIntensityAs::DATAVECTOR; - std::unique_ptr<MeshLib::Mesh> tmp_mesh (MeshLib::RasterToMesh::convert(data_array, header, meshElemType, useIntensity, var.name())); - - boost::optional< MeshLib::PropertyVector<double>& > tmp_vec = tmp_mesh->getProperties().getPropertyVector<double>(var.name()); - if (!tmp_vec) + if (mult_files) { - ERR("Vector not found"); - continue; + mesh.reset(MeshLib::RasterToMesh::convert(data_array, header, meshElemType, useIntensity, var.name())); + std::string const output_file_name (BaseLib::dropFileExtension(output_name) + getIterationString(i, time_bounds.second) + ".vtu"); + MeshLib::IO::VtuInterface vtu(mesh.get()); + vtu.writeToFile(output_file_name); } - - if (tmp_vec->size() != mesh->getNumberOfElements()) + else { - ERR("Vector size doesn't fit."); - continue; + std::string array_name (var.name()); + if (time_bounds.first != time_bounds.second) + array_name.append(getIterationString(i, time_bounds.second)); + if (i==time_bounds.first) // create persistent mesh + mesh.reset(MeshLib::RasterToMesh::convert(data_array, header, meshElemType, useIntensity, array_name)); + else + { + //copy array to mesh + } + if (i==time_bounds.second) + { + MeshLib::IO::VtuInterface vtu(mesh.get()); + vtu.writeToFile(output_name); + } } + } - MeshLib::Properties& props (mesh->getProperties()); - std::string prefix (""); - if (j<100) prefix = "0"; - if (j<10) prefix = "00"; - std::string var_name_str (prefix + std::to_string(j)); + /*********************** + for (std::size_t j = time_bounds.first; j <= time_bounds.second; ++j) + { + boost::optional< MeshLib::PropertyVector<double>& > tmp_vec = tmp_mesh->getProperties().getPropertyVector<double>(var.name()); + if (!tmp_vec) + { + ERR("Vector not found"); + continue; + } - boost::optional< MeshLib::PropertyVector<double>& > vec = - props.createNewPropertyVector<double>(std::string(var.name() + var_name_str), MeshLib::MeshItemType::Cell, 1); + if (tmp_vec->size() != mesh->getNumberOfElements()) + { + ERR("Vector size doesn't fit."); + continue; + } - if (!vec) - { - ERR ("Error creating vector"); - continue; - } - vec->reserve(tmp_vec->size()); - for (std::size_t i=0; i<tmp_vec->size(); ++i) - vec->push_back((*tmp_vec)[i]); + MeshLib::Properties& props (mesh->getProperties()); - std::cout << "time step " << j << " written.\n"; + std::string iter_str (getIterationString(j, time_bounds.second); + + boost::optional< MeshLib::PropertyVector<double>& > vec = + props.createNewPropertyVector<double>(std::string(var.name() + iter_str), MeshLib::MeshItemType::Cell, 1); + + if (!vec) + { + ERR ("Error creating vector"); + continue; + } + vec->reserve(tmp_vec->size()); + for (std::size_t i=0; i<tmp_vec->size(); ++i) + vec->push_back((*tmp_vec)[i]); + + std::cout << "time step " << j << " written.\n"; } /***********************/ - MeshLib::IO::VtuInterface vtu(mesh); - vtu.writeToFile("d:/nctest.vtu"); - delete[] length; delete[] data_array; - - return 0; }