From 6eb7135e1b1cd9f429e76ff49d8d758719cb3034 Mon Sep 17 00:00:00 2001
From: rinkk <karsten.rink@ufz.de>
Date: Thu, 5 Sep 2019 15:17:25 +0200
Subject: [PATCH] updated netcdf cli to netcdf4-cxx
---
.../Utils/FileConverter/CMakeLists.txt | 50 +---
.../Utils/FileConverter/NetCdfConverter.cpp | 274 +++++++++---------
2 files changed, 142 insertions(+), 182 deletions(-)
diff --git a/Applications/Utils/FileConverter/CMakeLists.txt b/Applications/Utils/FileConverter/CMakeLists.txt
index 230b4833461..4f326101ae1 100644
--- a/Applications/Utils/FileConverter/CMakeLists.txt
+++ b/Applications/Utils/FileConverter/CMakeLists.txt
@@ -9,7 +9,8 @@ set(TOOLS
TecPlotTools
GocadSGridReader
GocadTSurfaceReader
- Mesh2Raster)
+ Mesh2Raster
+ NetCdfConverter)
if(OGS_BUILD_GUI)
if(Shapelib_FOUND)
@@ -33,43 +34,14 @@ if(TARGET ConvertSHPToGLI)
target_link_libraries(ConvertSHPToGLI GeoLib Qt5::Xml ${Shapelib_LIBRARIES})
endif()
-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)
+if(TARGET Mesh2Shape)
+ target_link_libraries(Mesh2Shape ${Shapelib_LIBRARIES})
endif()
+if(TARGET NetCdfConverter)
+ target_link_libraries(NetCdfConverter
+ ${NETCDF_LIBRARIES_C}
+ ${NETCDF_LIBRARIES_CXX}
+ ${HDF5_LIBRARIES}
+ ${HDF5_HL_LIBRARIES})
+endif()
diff --git a/Applications/Utils/FileConverter/NetCdfConverter.cpp b/Applications/Utils/FileConverter/NetCdfConverter.cpp
index f68447b56eb..986cc6874cf 100644
--- a/Applications/Utils/FileConverter/NetCdfConverter.cpp
+++ b/Applications/Utils/FileConverter/NetCdfConverter.cpp
@@ -19,12 +19,12 @@
#include "Applications/ApplicationsLib/LogogSetup.h"
-#include <netcdfcpp.h>
+#include <netcdf>
// BaseLib
-#include "BaseLib/BuildInfo.h"
-#include "BaseLib/LogogSimpleFormatter.h"
#include "BaseLib/FileTools.h"
+#include "BaseLib/LogogSimpleFormatter.h"
+#include "InfoLib/GitInfo.h"
// GeoLib
#include "GeoLib/Raster.h"
@@ -35,10 +35,11 @@
#include "MeshLib/IO/VtkIO/VtuInterface.h"
+using namespace netCDF;
void checkExit(std::string const& str)
{
- if (str == "exit")
+ if (str == "x" || str == "exit")
exit(0);
}
@@ -87,60 +88,63 @@ std::size_t parseInput(std::string const& request_str, std::size_t max_val, bool
return std::numeric_limits<std::size_t>::max();
}
-void showArrays(NcFile const& dataset)
+NcVar getDimVar(NcFile const& dataset, NcVar const& var, std::size_t const dim)
{
- std::size_t const n_vars (dataset.num_vars());
+ NcDim const& dim_obj = var.getDim(dim);
+ return dataset.getVar(dim_obj.getName());
+}
+
+std::vector<std::string> showArrays(NcFile const& dataset)
+{
+ std::size_t const n_vars(dataset.getDimCount());
std::cout << "The NetCDF file contains the following " << n_vars << " arrays:\n\n";
std::cout << "\tIndex\tArray Name\t#Dimensions\n";
std::cout << "-------------------------------------------\n";
- for (std::size_t i=0; i<n_vars; ++i)
+ auto const& names = dataset.getVars();
+ std::vector<std::string> var_names;
+ for (auto [name, var] : names)
{
- NcVar const& var = *dataset.get_var(i);
- std::cout << "\t" << i << "\t" << var.name() << "\t(" << var.num_dims() << "D array)\n";
+ std::cout << "\t" << var_names.size() << "\t" << name << "\t(" << var.getDimCount() << "D array)\n";
+ var_names.push_back(name);
}
std::cout << "\n";
+ return var_names;
}
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());
+ std::cout << "Data array \"" << var.getName()
+ << "\" contains the following dimensions:\n";
+ std::size_t const n_dims (var.getDimCount());
for (std::size_t i=0; i<n_dims; ++i)
- std::cout << "\t" << i << "\t" << var.get_dim(i)->name() << "\t(" << var.get_dim(i)->size() << " values)\n";
+ std::cout << "\t" << i << "\t" << var.getDim(i).getName() << "\t(" << var.getDim(i).getSize() << " values)\n";
std::cout << "\n";
return true;
}
-std::size_t getDimensionBoundaries(NcFile const& dataset, NcVar const& var, std::size_t dim, double &start, double &end)
+std::pair<double, double> getBoundaries(NcVar const& var)
{
- NcVar& dim_var = *dataset.get_var(var.get_dim(dim)->name());
- if ((dim_var.num_dims()) == 1)
+ if ((var.getDimCount()) == 1)
{
- std::size_t size = dim_var.get_dim(0)->size();
-
- std::size_t length[1] = {1};
- long idx[1] = {0};
- dim_var.set_cur(idx);
- double val_at_idx[1] = {0};
- dim_var.get(val_at_idx, length);
- start = val_at_idx[0];
-
- idx[0] = size-1;
- dim_var.set_cur(idx);
- dim_var.get(val_at_idx, length);
- end = val_at_idx[0];
- return size;
+ double start, end;
+ std::size_t size = var.getDim(0).getSize();
+ var.getVar({0}, {1}, &start);
+ var.getVar({size - 1}, {1}, &end);
+ return std::make_pair(start, end);
}
- return 0;
+ return std::make_pair(0, 0);
}
-MathLib::Point3d getOrigin(std::array<std::pair<double, double>,3> bounds, std::size_t n_dims)
+MathLib::Point3d getOrigin(NcFile const& dataset, NcVar const& var)
{
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]);
+ std::size_t const n_dims = var.getDimCount();
+ for (std::size_t i=1; i<n_dims; ++i)
+ {
+ NcVar const& dim = getDimVar(dataset, var, var.getDimCount() - i);
+ auto const bounds = getBoundaries(dim);
+ origin[n_dims - i] = (bounds.first < bounds.second) ? bounds.first : bounds.second;
+ }
return origin;
}
@@ -165,21 +169,29 @@ void reverseNorthSouth(double* data, std::size_t width, std::size_t height)
delete[] cp_array;
}
-std::size_t arraySelectionLoop(NcFile const& dataset)
+bool canConvert(NcVar const& var)
+{
+ bool ret (true);
+ if (ret = (var.getDimCount() < 2))
+ ERR ("Only 2+ dimensional variables can be converted into OGS Meshes.\n");
+ return !ret;
+}
+
+std::string arraySelectionLoop(NcFile const& dataset)
{
- showArrays(dataset);
- std::size_t idx = parseInput("Enter data array index: ", dataset.num_vars(), true);
+ std::vector<std::string> const& names = showArrays(dataset);
+ std::size_t idx = parseInput("Enter data array index: ", dataset.getVarCount(), true);
- if (idx == dataset.num_vars())
+ if (idx == dataset.getVarCount() || !canConvert(dataset.getVar(names[idx])))
return arraySelectionLoop(dataset);
- return idx;
+ return names[idx];
}
bool dimensionSelectionLoop(NcVar const& var, std::array<std::size_t,4> &dim_idx_map)
{
showArraysDims(var);
- std::size_t const n_dims (var.num_dims());
+ std::size_t const n_dims(var.getDimCount());
dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
bool is_time_dep (true);
@@ -195,16 +207,17 @@ bool dimensionSelectionLoop(NcVar const& var, std::array<std::size_t,4> &dim_idx
std::stringstream str_stream(temp_str);
if (str_stream.str() == "c" || str_stream.str() == "continue")
is_time_dep = false;
- else
+ else
{
if (!(str_stream >> dim_idx_map[0]))
{
showErrorMessage(0);
+ dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
continue;
}
if (dim_idx_map[0] > n_dims-1)
{
- showErrorMessage(1, var.num_dims()-1);
+ showErrorMessage(1, var.getDimCount() - 1);
dim_idx_map[0] = std::numeric_limits<std::size_t>::max();
}
}
@@ -220,9 +233,9 @@ bool dimensionSelectionLoop(NcVar const& var, std::array<std::size_t,4> &dim_idx
dim_idx_map[i] = std::numeric_limits<std::size_t>::max();
std::string request_str ("Enter ID for dimension " + std::to_string(i) + ": ");
- std::size_t idx = parseInput(request_str, var.num_dims(), true);
+ std::size_t idx = parseInput(request_str, var.getDimCount(), true);
- if (idx == var.num_dims())
+ if (idx == var.getDimCount())
{
showArraysDims(var);
i--;
@@ -234,15 +247,16 @@ bool dimensionSelectionLoop(NcVar const& var, std::array<std::size_t,4> &dim_idx
return is_time_dep;
}
-std::pair<std::size_t, std::size_t> timestepSelectionLoop(NcFile const& dataset, NcVar const& var, std::size_t dim)
+std::pair<std::size_t, std::size_t> timestepSelectionLoop(NcFile const& dataset, NcVar const& var, std::size_t dim_idx)
{
- 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);
+ NcVar const& dim_var = getDimVar(dataset, var, dim_idx);
+ std::size_t const n_time_steps = var.getDim(dim_idx).getSize();
+ std::pair<std::size_t, std::size_t> bounds(
+ std::numeric_limits<std::size_t>::max(),
+ std::numeric_limits<std::size_t>::max());
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())
+ while (bounds.first > bounds.second || bounds.second > n_time_steps)
bounds.second = parseInput("Specify last time step to export: ", n_time_steps, false);
return bounds;
}
@@ -256,7 +270,7 @@ MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim)
while (t == MeshLib::MeshElemType::INVALID)
{
std::cout << "\nSelect element type for result, choose ";
-
+
if (dim==2) std::cout << "(t)riangle or (q)uadliteral: ";
if (dim==3) std::cout << "(p)rism or (h)exahedron: ";
std::string type ("");
@@ -291,7 +305,7 @@ MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim)
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::size_t n_time_steps(time_bounds.second - time_bounds.first + 1);
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";
@@ -299,16 +313,6 @@ bool multFilesSelectionLoop(std::pair<std::size_t, std::size_t> const& time_boun
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)
- newOrigin[i]=0;
- newOrigin[0] = time_step;
- var.set_cur(newOrigin);
- delete [] newOrigin;
-}
-
std::string getIterationString(std::size_t i, std::size_t max)
{
std::size_t const max_length (std::to_string(max).length());
@@ -316,11 +320,39 @@ std::string getIterationString(std::size_t i, std::size_t max)
return std::string(max_length - current_str.length(), '0') + current_str;
}
+double getResolution(NcFile const& dataset, NcVar const& var)
+{
+ std::size_t dim_idx = var.getDimCount() - 1;
+ auto const bounds = getBoundaries(getDimVar(dataset, var, dim_idx));
+ return fabs(bounds.second - bounds.first) / static_cast<double>(var.getDim(dim_idx).getSize());
+}
+
+GeoLib::RasterHeader createRasterHeader(
+ NcFile const& dataset, NcVar const& var,
+ std::array<std::size_t, 4> const& dim_idx_map,
+ std::vector<std::size_t> const& length, bool is_time_dep)
+{
+ MathLib::Point3d origin = getOrigin(dataset, var);
+ double const res = getResolution(dataset, var);
+ std::size_t n_dims = var.getDimCount();
+ std::size_t temp_offset = (is_time_dep) ? 1 : 0;
+ std::size_t z_length = (n_dims - temp_offset == 3) ? length[dim_idx_map.back()] : 1;
+ return {length[dim_idx_map[0 + temp_offset]], length[dim_idx_map[1 + temp_offset]], z_length, origin, res, -9999};
+}
+
int main (int argc, char* argv[])
{
ApplicationsLib::LogogSetup logog_setup;
- TCLAP::CmdLine cmd("Converts NetCDF into mesh file.", ' ', BaseLib::BuildInfo::git_describe);
+ TCLAP::CmdLine cmd(
+ "Converts NetCDF into mesh file(s).\n\n "
+ "OpenGeoSys-6 software, version " +
+ GitInfoLib::GitInfo::ogs_version +
+ ".\n"
+ "Copyright (c) 2012-2019, OpenGeoSys Community "
+ "(http://www.opengeosys.org)",
+ ' ', GitInfoLib::GitInfo::ogs_version);
+
TCLAP::ValueArg<std::string> input("i", "input-file",
"the NetCDF input file", true,
"", "the NetCDF input file");
@@ -331,9 +363,9 @@ int main (int argc, char* argv[])
cmd.add(output);
cmd.parse(argc, argv);
- NcFile dataset(input.getValue().c_str(), NcFile::ReadOnly);
+ NcFile dataset(input.getValue().c_str(), NcFile::read);
- if (!dataset.is_valid())
+ if (dataset.isNull())
{
ERR("Error opening file.");
return -1;
@@ -344,34 +376,25 @@ int main (int argc, char* argv[])
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::string const& var_name = arraySelectionLoop(dataset);
+ NcVar& var = dataset.getVar(var_name);
std::array<std::size_t, 4> dim_idx_map;
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<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 const temp_offset = (is_time_dep) ? 1 : 0;
+ std::size_t const n_dims = (var.getDimCount());
+ std::vector<std::size_t> offset(n_dims, 0);
+ std::vector<std::size_t> length(n_dims, 1);
std::size_t array_length = 1;
for (std::size_t i=temp_offset; i<n_dims; ++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]];
+ length[i] = var.getDim(i).getSize();
+ array_length *= length[i];
}
- 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;
+ GeoLib::RasterHeader header =
+ createRasterHeader(dataset, var, dim_idx_map, length, is_time_dep);
std::pair<std::size_t, std::size_t> time_bounds(0,0);
if (is_time_dep)
@@ -381,42 +404,45 @@ int main (int argc, char* argv[])
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;
+ MeshLib::MeshElemType const meshElemType = elemSelectionLoop(n_dims - temp_offset);
for (std::size_t i=time_bounds.first; i<=time_bounds.second; ++i)
{
- if (is_time_dep)
- setTimeStep(var, i);
-
- 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);
-
- 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"
+ std::cout << "Converting time step " << i << "...\n";
+ offset[0] = i;
+ std::vector<double> data_array(array_length,0);
+ var.getVar(offset, length, data_array.data());
+ std::replace_if(data_array.begin(), data_array.end(),
+ [](double& x) { return x <= -9999; }, -9999);
// reverse lines in vertical direction if the original file has its origin in the northwest corner
- if (spatial_bounds[0].first > spatial_bounds[0].second)
- reverseNorthSouth(data_array, spatial_size[1], spatial_size[0]);
+ NcVar const& dim_var = getDimVar(dataset, var, n_dims - 1);
+ auto const bounds = getBoundaries(dim_var);
+ if (bounds.first > bounds.second)
+ reverseNorthSouth(data_array.data(), length[n_dims-2], length[n_dims-1]);
+ MeshLib::UseIntensityAs const useIntensity = MeshLib::UseIntensityAs::DATAVECTOR;
if (mult_files)
{
- mesh.reset(MeshLib::RasterToMesh::convert(data_array, header, meshElemType, useIntensity, var.name()));
+ mesh.reset(MeshLib::RasterToMesh::convert(data_array.data(), header, meshElemType, useIntensity, var.getName()));
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);
}
else
{
- std::string array_name (var.name());
+ std::string array_name (var.getName());
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
+ mesh.reset(MeshLib::RasterToMesh::convert(data_array.data(), header, meshElemType, useIntensity, array_name));
+ else // copy array to mesh
{
- //copy array to mesh
+ std::unique_ptr<MeshLib::Mesh> temp (MeshLib::RasterToMesh::convert(data_array.data(), header, meshElemType, useIntensity, array_name));
+ MeshLib::PropertyVector<double> const*const vec = temp->getProperties().getPropertyVector<double>(array_name);
+ if (vec == nullptr)
+ return EXIT_FAILURE;
+ MeshLib::addPropertyToMesh<double>(*mesh, array_name, MeshLib::MeshItemType::Cell, 1, *vec);
}
if (i==time_bounds.second)
{
@@ -426,44 +452,6 @@ int main (int argc, char* argv[])
}
}
-
- /***********************
- 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;
- }
-
- if (tmp_vec->size() != mesh->getNumberOfElements())
- {
- ERR("Vector size doesn't fit.");
- continue;
- }
-
- MeshLib::Properties& props (mesh->getProperties());
-
- 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";
- }
- /***********************/
-
- delete[] length;
- delete[] data_array;
- return 0;
+ std::cout << "Conversion finished successfully.\n";
+ return EXIT_SUCCESS;
}
--
GitLab