diff --git a/Applications/Utils/FileConverter/NetCdfConverter.cpp b/Applications/Utils/FileConverter/NetCdfConverter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..abd41633a53a8daf4857ed070b013af27896bff7 --- /dev/null +++ b/Applications/Utils/FileConverter/NetCdfConverter.cpp @@ -0,0 +1,419 @@ +/** + * @copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/LICENSE.txt + */ + +// STL +#include <cctype> +#include <iostream> +#include <limits> +#include <memory> +#include <sstream> +#include <string> +#include <utility> + +#include <tclap/CmdLine.h> + +#include "Applications/ApplicationsLib/LogogSetup.h" + +#include <netcdfcpp.h> + +// BaseLib +#include "BaseLib/BuildInfo.h" +#include "BaseLib/LogogSimpleFormatter.h" + +// GeoLib +#include "GeoLib/Raster.h" + +// MeshLib +#include "MeshLib/Mesh.h" +#include "MeshLib/MeshGenerators/RasterToMesh.h" + +#include "MeshLib/IO/VtkIO/VtuInterface.h" + + +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"; + } + else if (error_id == 1) + { + ERR ("Index not valid. Valid indices are in [0,%d].", max); + } +} + +void showArrays(NcFile const& dataset) +{ + std::size_t const n_vars (dataset.num_vars()); + 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) + { + NcVar const& var = *dataset.get_var(i); + std::cout << "\t" << i << "\t" << var.name() << "\t(" << var.num_dims() << "D array)\n"; + } + std::cout << "\n"; +} + +bool showArraysDims(NcFile const& dataset, NcVar const& var) +{ + std::cout << "Data array \"" << var.name() << "\" contains the following dimensions:\n"; + std::size_t const n_dims (var.num_dims()); + 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 << "\n"; + return true; +} + +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()); + if ((dim_var.num_dims()) == 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; + } + return 0; +} + +void reverseNorthSouth(double* data, std::size_t width, std::size_t height) +{ + double* cp_array = new double[width*height]; + + 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 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++) + data[i] = cp_array[i]; + + delete[] cp_array; +} + +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(); +} + +bool dimensionSelectionLoop(NcFile const& dataset, NcVar const& var, std::array<std::size_t,4> &dim_idx_map) +{ + showArraysDims(dataset, 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); + + // get temporal dimension + if (n_dims > 1) + { + std::string temp_str(""); + cout << "Is the parameter time-dependent?\n"; + while (dim_idx_map[0] == std::numeric_limits<std::size_t>::max() && is_time_dep == true) + { + 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") + is_time_dep = false; + else + { + if (str_stream >> dim_idx_map[0]) + { + 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(); + } + } + } + } + } + + // 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; + } + } + + return is_time_dep; +} + +MeshLib::MeshElemType elemSelectionLoop(std::size_t const dim) +{ + if (dim ==1) + return MeshLib::MeshElemType::LINE; + + MeshLib::MeshElemType t = MeshLib::MeshElemType::INVALID; + 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 (""); + std::getline(std::cin, type); + + if (dim==2) + { + 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; + else + return MeshLib::MeshElemType::QUAD; + } + + if (dim==3) + { + if (type != "p" || type != "h" || + type != "prism" || type != "hex" || + type != "hexahedron") + continue; + if (type == "p" || type == "prism") + return MeshLib::MeshElemType::PRISM; + else + return MeshLib::MeshElemType::HEXAHEDRON; + } + } + return t; +} + +void setOrigin(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; +} + +/* +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) +{ + return +} +*/ +int main (int argc, char* argv[]) +{ + ApplicationsLib::LogogSetup logog_setup; + + TCLAP::CmdLine cmd("Converts NetCDF into mesh file.", ' ', BaseLib::BuildInfo::git_describe); + TCLAP::ValueArg<std::string> input("i", "input-file", + "the NetCDF input file", true, + "", "the NetCDF input file"); + cmd.add(input); + TCLAP::ValueArg<std::string> output("o", "output-file", + "the OGS mesh output file", true, + "", "the OGS mesh output file"); + cmd.add(output); + cmd.parse(argc, argv); + + NcFile dataset(input.getValue().c_str(), NcFile::ReadOnly); + + std::cout << "OpenGeoSys NetCDF Converter\n"; + std::cout << "File " << input.getValue() << " loaded. Press ENTER to display available data arrays.\n"; + std::cin.ignore(); + + std::size_t const var_idx = arraySelectionLoop(dataset); + + 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(); + + 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::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]; + } + + 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; + MeshLib::UseIntensityAs useIntensity = MeshLib::UseIntensityAs::DATAVECTOR; + + MeshLib::Mesh* mesh = MeshLib::RasterToMesh::convert(data_array, header, meshElemType, useIntensity, var.name()); + + /*********************** + 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) + { + setOrigin(var, j); + var.get(data_array, length); //create Array of Values + + 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" + + double const resolution = fabs(end_lat-start_lat)/(sizeLat-1); + + // 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}; + + 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) + { + 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 prefix (""); + if (j<100) prefix = "0"; + if (j<10) prefix = "00"; + std::string var_name_str (prefix + std::to_string(j)); + + boost::optional< MeshLib::PropertyVector<double>& > vec = + props.createNewPropertyVector<double>(std::string(var.name() + var_name_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; +}