OGS
createIntermediateRasters.cpp
Go to the documentation of this file.
1
10#include <tclap/CmdLine.h>
11
12#ifdef USE_PETSC
13#include <mpi.h>
14#endif
15
16#include <memory>
17#include <string>
18
19#include "BaseLib/FileTools.h"
21#include "GeoLib/Raster.h"
22#include "InfoLib/GitInfo.h"
23
24int main(int argc, char* argv[])
25{
26 TCLAP::CmdLine cmd(
27 "Takes two DEMs located at the exact same spatial position (but at "
28 "different elevation) and calculates n raster DEMs located at "
29 "equidistant intervals between them (i.e. for n=1, one new raster "
30 "located precisely in the middle will be created).\n\n"
31 "OpenGeoSys-6 software, version " +
33 ".\n"
34 "Copyright (c) 2012-2024, OpenGeoSys Community "
35 "(http://www.opengeosys.org)",
37 TCLAP::ValueArg<std::size_t> number_arg(
38 "n", "number", "number of rasters to be calculated", false, 1, "int");
39 cmd.add(number_arg);
40 TCLAP::ValueArg<std::string> output_arg("o", "output-file",
41 "Raster output file (*.asc)", true,
42 "", "output.asc");
43 cmd.add(output_arg);
44 TCLAP::ValueArg<std::string> input2_arg(
45 "", "file2", "Second DEM-raster file", true, "", "file2.asc");
46 cmd.add(input2_arg);
47 TCLAP::ValueArg<std::string> input1_arg(
48 "", "file1", "First DEM-raster file", true, "", "file1.asc");
49 cmd.add(input1_arg);
50
51 cmd.parse(argc, argv);
52
53#ifdef USE_PETSC
54 MPI_Init(&argc, &argv);
55#endif
56
57 std::unique_ptr<GeoLib::Raster> dem1(
58 FileIO::AsciiRasterInterface::readRaster(input1_arg.getValue()));
59 std::unique_ptr<GeoLib::Raster> dem2(
60 FileIO::AsciiRasterInterface::readRaster(input2_arg.getValue()));
61
62 if (dem1 == nullptr || dem2 == nullptr)
63 {
64#ifdef USE_PETSC
65 MPI_Finalize();
66#endif
67 return 1;
68 }
69
70 GeoLib::RasterHeader const& h1 = dem1->getHeader();
71 GeoLib::RasterHeader const& h2 = dem2->getHeader();
72
73 bool errors_found(false);
74 if (std::abs(h1.origin[0] - h2.origin[0]) > (h1.cell_size / 100.0))
75 {
76 ERR("Origin x-coordinate is not the same in both raster files.\n");
77 errors_found = true;
78 }
79 if (std::abs(h1.origin[1] - h2.origin[1]) > (h1.cell_size / 100.0))
80 {
81 ERR("Origin y-coordinate is not the same in both raster files.\n");
82 errors_found = true;
83 }
84 if (h1.cell_size != h2.cell_size)
85 {
86 ERR("Cellsize is not the same in both raster files.\n");
87 errors_found = true;
88 }
89 if (h1.n_cols != h2.n_cols)
90 {
91 ERR("Raster width is not the same in both raster files.\n");
92 errors_found = true;
93 }
94 if (h1.n_rows != h2.n_rows)
95 {
96 ERR("Raster height is not the same in both raster files.\n");
97 errors_found = true;
98 }
99
100 if (errors_found)
101 {
102#ifdef USE_PETSC
103 MPI_Finalize();
104#endif
105 return 2;
106 }
107
108 std::size_t const n = number_arg.getValue();
109 std::vector<std::vector<double>> raster;
110 for (std::size_t i = 0; i < n; ++i)
111 {
112 std::vector<double> r;
113 r.reserve(h1.n_cols * h1.n_rows);
114 raster.push_back(r);
115 }
116
117 auto it2 = dem2->begin();
118 for (auto it1 = dem1->begin(); it1 != dem1->end(); ++it1)
119 {
120 if (it2 == dem2->end())
121 {
122 ERR("Error: File 2 is shorter than File 1.");
123#ifdef USE_PETSC
124 MPI_Finalize();
125#endif
126 return 3;
127 }
128 if (*it1 == h1.no_data || *it2 == h2.no_data)
129 {
130 for (std::size_t i = 0; i < n; ++i)
131 {
132 raster[i].push_back(h1.no_data);
133 }
134 }
135 else
136 {
137 double const min = std::min(*it1, *it2);
138 double const max = std::max(*it1, *it2);
139 double const step = (max - min) / static_cast<double>(n + 1);
140 for (std::size_t i = 0; i < n; ++i)
141 {
142 raster[i].push_back(max - ((i + 1) * step));
143 }
144 }
145 it2++;
146 }
147 if (it2 != dem2->end())
148 {
149 ERR("Error: File 1 is shorter than File 2.");
150#ifdef USE_PETSC
151 MPI_Finalize();
152#endif
153 return 3;
154 }
155
156 std::string const filename = output_arg.getValue();
157 for (std::size_t i = 0; i < n; ++i)
158 {
159 std::string const basename = BaseLib::dropFileExtension(filename);
160 std::string const ext = BaseLib::getFileExtension(filename);
161
162 GeoLib::Raster r(h1, raster[i].begin(), raster[i].end());
164 r, basename + std::to_string(i) + ext);
165 INFO("Layer {:d} written.", i + 1);
166 }
167#ifdef USE_PETSC
168 MPI_Finalize();
169#endif
170 return EXIT_SUCCESS;
171}
Definition of the AsciiRasterInterface class.
Filename manipulation routines.
Git information.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Definition of the GeoLib::Raster class.
static void writeRasterAsASC(GeoLib::Raster const &raster, std::string const &file_name)
Writes an Esri asc-file.
static GeoLib::Raster * readRaster(std::string const &fname)
Class Raster is used for managing raster data.
Definition Raster.h:49
int main(int argc, char *argv[])
std::string getFileExtension(const std::string &path)
std::string dropFileExtension(std::string const &filename)
GITINFOLIB_EXPORT const std::string ogs_version
Contains the relevant information when storing a geoscientific raster data.
Definition Raster.h:28
MathLib::Point3d origin
Definition Raster.h:32
std::size_t n_cols
Definition Raster.h:29
std::size_t n_rows
Definition Raster.h:30