Merge branch 'RasterParameter' into 'master'

Implementation of parameter type Raster

See merge request ogs/ogs!4788

Applications/ApplicationsLib/ProjectData.cpp

@@ -31,6 +31,9 @@
 #include "BaseLib/Logging.h"
 #include "BaseLib/StringTools.h"
 #include "GeoLib/GEOObjects.h"
+#include "GeoLib/IO/AsciiRasterInterface.h"
+#include "GeoLib/IO/NetCDFRasterReader.h"
+#include "GeoLib/Raster.h"
 #include "InfoLib/CMakeInfo.h"
 #include "MaterialLib/MPL/CreateMedium.h"
 #include "MathLib/Curve/CreatePiecewiseLinearCurve.h"
@@ -271,6 +274,52 @@ std::vector<std::unique_ptr<MeshLib::Mesh>> readMeshes(
     return meshes;
 }
 
+std::vector<GeoLib::NamedRaster> readRasters(
+    BaseLib::ConfigTree const& config, std::string const& raster_directory,
+    GeoLib::MinMaxPoints const& min_max_points)
+{
+    INFO("readRasters ... ");
+    std::vector<GeoLib::NamedRaster> named_rasters;
+
+    //! \ogs_file_param{prj__rasters}
+    auto optional_rasters = config.getConfigSubtreeOptional("rasters");
+    if (optional_rasters)
+    {
+        INFO("Reading rasters.");
+        //! \ogs_file_param{prj__rasters__raster}
+        auto const configs = optional_rasters->getConfigSubtreeList("raster");
+        std::transform(
+            configs.begin(), configs.end(), std::back_inserter(named_rasters),
+            [&raster_directory, &min_max_points](auto const& raster_config)
+            {
+                return GeoLib::IO::readRaster(raster_config, raster_directory,
+                                              min_max_points);
+            });
+    }
+    INFO("readRasters done");
+    return named_rasters;
+}
+
+// for debugging raster reading implementation
+// void writeRasters(std::vector<GeoLib::NamedRaster> const& named_rasters,
+//                  std::string const& output_directory)
+//{
+//    for (auto const& named_raster : named_rasters)
+//    {
+// #if defined(USE_PETSC)
+//        int my_mpi_rank;
+//        MPI_Comm_rank(MPI_COMM_WORLD, &my_mpi_rank);
+// #endif
+//        FileIO::AsciiRasterInterface::writeRasterAsASC(
+//            named_raster.raster, output_directory + "/" +
+//                                     named_raster.raster_name +
+// #if defined(USE_PETSC)
+//                                     "_" + std::to_string(my_mpi_rank) +
+// #endif
+//                                     ".asc");
+//    }
+//}
+
 std::optional<ParameterLib::CoordinateSystem> parseLocalCoordinateSystem(
     std::optional<BaseLib::ConfigTree> const& config,
     std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
@@ -336,8 +385,14 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
                          std::string const& output_directory,
                          std::string const& mesh_directory,
                          [[maybe_unused]] std::string const& script_directory)
-    : _mesh_vec(readMeshes(project_config, mesh_directory))
+    : _mesh_vec(readMeshes(project_config, mesh_directory)),
+      _named_rasters(readRasters(project_config, project_directory,
+                                 GeoLib::AABB(_mesh_vec[0]->getNodes().begin(),
+                                              _mesh_vec[0]->getNodes().end())
+                                     .getMinMaxPoints()))
 {
+    // for debugging raster reading implementation
+    // writeRasters(_named_rasters, output_directory);
     if (auto const python_script =
             //! \ogs_file_param{prj__python_script}
         project_config.getConfigParameterOptional<std::string>("python_script"))
@@ -475,8 +530,8 @@ std::vector<std::string> ProjectData::parseParameters(
          //! \ogs_file_param{prj__parameters__parameter}
          parameters_config.getConfigSubtreeList("parameter"))
     {
-        auto p =
-            ParameterLib::createParameter(parameter_config, _mesh_vec, _curves);
+        auto p = ParameterLib::createParameter(parameter_config, _mesh_vec,
+                                               _named_rasters, _curves);
         if (!names.insert(p->name).second)
         {
             OGS_FATAL("A parameter with name `{:s}' already exists.", p->name);

Applications/ApplicationsLib/ProjectData.h

@@ -32,6 +32,11 @@
 #include <pybind11/embed.h>
 #endif
 
+namespace GeoLib
+{
+struct NamedRaster;
+}
+
 namespace MeshLib
 {
 class Mesh;
@@ -138,6 +143,7 @@ class ProjectData final
         const std::string& output_directory);
 
     std::vector<std::unique_ptr<MeshLib::Mesh>> _mesh_vec;
+    std::vector<GeoLib::NamedRaster> _named_rasters;
     std::vector<std::unique_ptr<ProcessLib::Process>> _processes;
     std::vector<ProcessLib::ProcessVariable> _process_variables;
 

Applications/Utils/GeoTools/addDataToRaster.cpp

@@ -51,6 +51,24 @@ double computeSinXSinY(GeoLib::Point const& point, GeoLib::AABB const& aabb)
                     boost::math::double_constants::pi);
 }
 
+double computeStepFunction(GeoLib::Point const& point, GeoLib::AABB const& aabb)
+{
+    auto const aabb_size = aabb.getMaxPoint() - aabb.getMinPoint();
+    auto const min = aabb.getMinPoint();
+
+    if ((point[0] - min[0]) / aabb_size[0] < 0.5 &&
+        (point[1] - min[1]) / aabb_size[1] < 0.5)
+    {
+        return 1;
+    }
+    if ((point[0] - min[0]) / aabb_size[0] >= 0.5 &&
+        (point[1] - min[1]) / aabb_size[1] >= 0.5)
+    {
+        return 1;
+    }
+    return 0;
+}
+
 int main(int argc, char* argv[])
 {
     TCLAP::CmdLine cmd(
@@ -122,7 +140,8 @@ int main(int argc, char* argv[])
                                      "double value");
 
     cmd.add(ur_y_arg);
-    std::vector<std::string> allowed_functions_vector{"sinxsiny", "exp"};
+    std::vector<std::string> allowed_functions_vector{"sinxsiny", "exp",
+                                                      "step"};
     TCLAP::ValuesConstraint<std::string> allowed_functions(
         allowed_functions_vector);
     TCLAP::ValueArg<std::string> function_arg(
@@ -154,17 +173,34 @@ int main(int argc, char* argv[])
     auto const& header = raster->getHeader();
     auto const& origin = header.origin;
 
+    auto function_selector = [](std::string const& function_string)
+        -> std::function<double(GeoLib::Point const& p,
+                                GeoLib::AABB const& aabb)>
+    {
+        if (function_string == "sinxsiny")
+        {
+            return computeSinXSinY;
+        }
+        if (function_string == "exp")
+        {
+            return compute2DGaussBellCurveValues;
+        }
+        if (function_string == "step")
+        {
+            return computeStepFunction;
+        }
+        OGS_FATAL("Function '{}' isn't implemented.", function_string);
+    };
     std::function<double(GeoLib::Point const& p, GeoLib::AABB const& aabb)>
-        computeFunctionValue = function_arg.getValue() == "sinxsiny"
-                                   ? computeSinXSinY
-                                   : compute2DGaussBellCurveValues;
+        computeFunctionValue = function_selector(function_arg.getValue());
 
     for (std::size_t r = 0; r < header.n_rows; r++)
     {
         for (std::size_t c = 0; c < header.n_cols; c++)
         {
-            GeoLib::Point const p{{origin[0] + header.cell_size * c,
-                                   origin[1] + header.cell_size * r, 0.0}};
+            GeoLib::Point const p{{origin[0] + header.cell_size * (c + 0.5),
+                                   origin[1] + header.cell_size * (r + 0.5),
+                                   0.0}};
             if (!aabb.containsPoint(p, std::numeric_limits<double>::epsilon()))
             {
                 continue;

Applications/Utils/Tests.cmake

@@ -1145,6 +1145,17 @@ AddTest(
     10x10_sinxsiny.asc 10x10_sinxsiny.asc
 )
 
+AddTest(
+    NAME addDataToRaster_10x10_step_function
+    PATH GeoTools/addDataToRaster/
+    WORKING_DIRECTORY ${Data_SOURCE_DIR}/GeoTools/addDataToRaster
+    EXECUTABLE addDataToRaster
+    EXECUTABLE_ARGS --ll_x 1000 --ll_y 100 --ur_x 1100 --ur_y 200 --function step --scaling_value 1 --offset_value 0 -i ${Data_SOURCE_DIR}/GeoTools/createRaster/10x10.asc -o ${Data_BINARY_DIR}/GeoTools/addDataToRaster/10x10_step.asc
+    TESTER diff
+    DIFF_DATA
+    10x10_step.asc 10x10_step.asc
+)
+
 AddTest(
     NAME GMSH2OGS_linearElements
     PATH Utils/GMSH2OGS

Documentation/.vale/Vocab/ogs/accept.txt

@@ -90,6 +90,7 @@ MODFLOW
 MTest
 nbconvert
 nbdime
+NetCDF
 Netbeans
 NetBeans
 Neumann

Documentation/ProjectFile/prj/rasters/i_rasters.md

@@ -0,0 +1 @@
+In the 'rasters' tag, a list of rasters to load is specified.

Documentation/ProjectFile/prj/rasters/raster/i_raster.md

@@ -0,0 +1,7 @@
+In the tag 'raster' the raster to be loaded is specified, i.e. the file name,
+the variable name and the dimension. The parameter name is generated following
+this scheme
+
+`FilenameWithoutExtension_VariableName_Dimension`
+
+This name can then be used to define a parameter of type 'Raster'.

Documentation/ProjectFile/prj/rasters/raster/t_dimension.md

@@ -0,0 +1,2 @@
+The 'dimension' tag is only relevant for NetCDF interface. However, the value of
+the tag is used to form the parameter name.

Documentation/ProjectFile/prj/rasters/raster/t_file.md

@@ -0,0 +1 @@
+The 'file' tag specifies the file name.

Documentation/ProjectFile/prj/rasters/raster/t_variable.md

@@ -0,0 +1,2 @@
+The 'variable' tag is only relevant for the NetCDF interface. However, the value
+of the tag is used to form the parameter name.

GeoLib/CMakeLists.txt

@@ -20,6 +20,14 @@ target_link_libraries(
     PRIVATE range-v3 tet
 )
 
+if(OGS_USE_NETCDF)
+    target_link_libraries(
+        GeoLib
+        PRIVATE ${NETCDF_LIBRARIES_CXX} ${NETCDF_LIBRARIES_C}
+    )
+    target_compile_definitions(GeoLib PRIVATE OGS_USE_NETCDF)
+endif()
+
 foreach(xsd OpenGeoSysGLI.xsd OpenGeoSysSTN.xsd)
     # cmake-lint: disable=E1126
     file(COPY_FILE IO/XmlIO/${xsd} ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${xsd}