OGS
addDataToRaster.cpp
Go to the documentation of this file.
1
11#include <tclap/CmdLine.h>
12
13#ifdef USE_PETSC
14#include <mpi.h>
15#endif
16
17#include <algorithm>
18#include <boost/math/constants/constants.hpp>
19#include <cmath>
20#include <memory>
21#include <numeric>
22
24#include "GeoLib/AABB.h"
25#include "GeoLib/Point.h"
26#include "GeoLib/Raster.h"
27
29 GeoLib::AABB const& aabb)
30{
31 auto const sigma_x = (aabb.getMaxPoint() - aabb.getMinPoint())[0] / 3;
32 auto const sigma_y = (aabb.getMaxPoint() - aabb.getMinPoint())[1] / 3;
33
34 auto const mid_point = (aabb.getMaxPoint() + aabb.getMinPoint()) / 2;
35
36 return std::exp(
37 -0.5 * std::pow((point[0] - mid_point[0]), 2) / std::pow(sigma_x, 2) -
38 0.5 * std::pow((point[1] - mid_point[1]), 2) / std::pow(sigma_y, 2));
39}
40
41double computeSinXSinY(GeoLib::Point const& point, GeoLib::AABB const& aabb)
42{
43 auto const aabb_size = aabb.getMaxPoint() - aabb.getMinPoint();
44 auto const offset = aabb.getMinPoint();
45
46 return std::sin((point[0] - offset[0]) / aabb_size[0] *
47 boost::math::double_constants::pi) *
48 std::sin((point[1] - offset[1]) / aabb_size[1] *
49 boost::math::double_constants::pi);
50}
51
52int main(int argc, char* argv[])
53{
54 TCLAP::CmdLine cmd("Add values to raster.", ' ', "0.1");
55
56 TCLAP::ValueArg<std::string> out_raster_arg(
57 "o",
58 "output_raster",
59 "the output raster is stored to a file of this name",
60 true,
61 "",
62 "filename for raster output");
63 cmd.add(out_raster_arg);
64
65 TCLAP::ValueArg<double> scaling_arg(
66 "",
67 "scaling_value",
68 "value the function sin(x pi) sin(y pi) will be scaled with",
69 false,
70 1,
71 "double value");
72 cmd.add(scaling_arg);
73
74 TCLAP::ValueArg<double> offset_arg(
75 "",
76 "offset_value",
77 "constant added to the function 'scaling * sin(x pi) * sin(y pi)'",
78 false,
79 0,
80 "double value");
81 cmd.add(offset_arg);
82
83 TCLAP::ValueArg<double> ll_x_arg(
84 "",
85 "ll_x",
86 "x coordinate of lower left point of axis aligned rectangular region",
87 false,
88 0,
89 "double value");
90 cmd.add(ll_x_arg);
91 TCLAP::ValueArg<double> ll_y_arg(
92 "",
93 "ll_y",
94 "y coordinate of lower left point of axis aligned rectangular region",
95 false,
96 0,
97 "double value");
98 cmd.add(ll_y_arg);
99 TCLAP::ValueArg<double> ur_x_arg("",
100 "ur_x",
101 "x coordinate of the upper right point of "
102 "axis aligned rectangular region",
103 false,
104 0,
105 "double value");
106 cmd.add(ur_x_arg);
107 TCLAP::ValueArg<double> ur_y_arg("",
108 "ur_y",
109 "y coordinate of the upper right point of "
110 "axis aligned rectangular region",
111 false,
112 0,
113 "double value");
114
115 cmd.add(ur_y_arg);
116 std::vector<std::string> allowed_functions_vector{"sinxsiny", "exp"};
117 TCLAP::ValuesConstraint<std::string> allowed_functions(
118 allowed_functions_vector);
119 TCLAP::ValueArg<std::string> function_arg(
120 "f", "function", "Name of the function used to modify the raster", true,
121 "", &allowed_functions);
122 cmd.add(function_arg);
123 TCLAP::ValueArg<std::string> input_arg("i", "input",
124 "Name of the input raster (*.asc)",
125 true, "", "input file name");
126 cmd.add(input_arg);
127
128 cmd.parse(argc, argv);
129
130#ifdef USE_PETSC
131 MPI_Init(&argc, &argv);
132#endif
133
134 std::array input_points = {
135 GeoLib::Point{{ll_x_arg.getValue(), ll_y_arg.getValue(), 0}},
136 GeoLib::Point{{ur_x_arg.getValue(), ur_y_arg.getValue(), 0}}};
137 GeoLib::AABB const aabb{std::begin(input_points), std::end(input_points)};
138
139 auto const s = scaling_arg.getValue();
140 auto const offset = offset_arg.getValue();
141
142 std::unique_ptr<GeoLib::Raster> const raster(
144 input_arg.getValue()));
145 auto const& header = raster->getHeader();
146 auto const& origin = header.origin;
147
148 std::function<double(GeoLib::Point const& p, GeoLib::AABB const& aabb)>
149 computeFunctionValue = function_arg.getValue() == "sinxsiny"
152
153 for (std::size_t r = 0; r < header.n_rows; r++)
154 {
155 for (std::size_t c = 0; c < header.n_cols; c++)
156 {
157 GeoLib::Point const p{{origin[0] + header.cell_size * c,
158 origin[1] + header.cell_size * r, 0.0}};
159 if (!aabb.containsPoint(p, std::numeric_limits<double>::epsilon()))
160 {
161 continue;
162 }
163
164 (*raster)(r, c) += offset + s * computeFunctionValue(p, aabb);
165 }
166 }
167
169 out_raster_arg.getValue());
170#ifdef USE_PETSC
171 MPI_Finalize();
172#endif
173 return EXIT_SUCCESS;
174}
Definition of the AABB class.
Definition of the AsciiRasterInterface class.
Definition of the Point class.
Definition of the GeoLib::Raster class.
int main(int argc, char *argv[])
double computeSinXSinY(GeoLib::Point const &point, GeoLib::AABB const &aabb)
double compute2DGaussBellCurveValues(GeoLib::Point const &point, GeoLib::AABB const &aabb)
static void writeRasterAsASC(GeoLib::Raster const &raster, std::string const &file_name)
Writes an Esri asc-file.
static GeoLib::Raster * getRasterFromASCFile(std::string const &fname)
Reads an ArcGis ASC raster file.
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition: AABB.h:49
Eigen::Vector3d const & getMaxPoint() const
Definition: AABB.h:179
bool containsPoint(T const &pnt, double eps) const
Definition: AABB.h:136
Eigen::Vector3d const & getMinPoint() const
Definition: AABB.h:172
static const double p
static const double r
static const double s