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 
20 #include "BaseLib/FileTools.h"
21 #include "GeoLib/Raster.h"
22 #include "InfoLib/GitInfo.h"
23 
24 int 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-2022, 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:34
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:44
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)
Reads raster file by detecting type based on extension and then calling the appropriate method.
Class Raster is used for managing raster data.
Definition: Raster.h:45
int main(int argc, char *argv[])
std::string getFileExtension(const std::string &path)
Definition: FileTools.cpp:187
std::string dropFileExtension(std::string const &filename)
Definition: FileTools.cpp:170
GITINFOLIB_EXPORT const std::string ogs_version
static const double r
Contains the relevant information when storing a geoscientific raster data.
Definition: Raster.h:27
double cell_size
Definition: Raster.h:32
MathLib::Point3d origin
Definition: Raster.h:31
std::size_t n_cols
Definition: Raster.h:28
std::size_t n_rows
Definition: Raster.h:29