OGS
AsciiRasterInterface.cpp
Go to the documentation of this file.
1
15
16#include <fstream>
17#include <tuple>
18
19#include "BaseLib/FileTools.h"
20#include "BaseLib/Logging.h"
21#include "BaseLib/StringTools.h"
22#include "GeoLib/Point.h"
23
24namespace FileIO
25{
27{
28 std::string ext(BaseLib::getFileExtension(fname));
29 std::transform(ext.begin(), ext.end(), ext.begin(), tolower);
30 if (ext == ".asc")
31 {
32 return getRasterFromASCFile(fname);
33 }
34 if (ext == ".grd")
35 {
36 return getRasterFromSurferFile(fname);
37 }
38 if (ext == ".xyz")
39 {
40 return getRasterFromXyzFile(fname);
41 }
42 return nullptr;
43}
44
46static double readDoubleFromStream(std::istream& in)
47{
48 std::string value;
49 in >> value;
50 return std::strtod(BaseLib::replaceString(",", ".", value).c_str(),
51 nullptr);
52}
53
56static std::optional<GeoLib::RasterHeader> readASCHeader(std::ifstream& in)
57{
59
60 std::string tag;
61 std::string value;
62
63 in >> tag;
64 if (tag == "ncols")
65 {
66 in >> value;
67 header.n_cols = atoi(value.c_str());
68 }
69 else
70 {
71 return {};
72 }
73
74 in >> tag;
75 if (tag == "nrows")
76 {
77 in >> value;
78 header.n_rows = atoi(value.c_str());
79 }
80 else
81 {
82 return {};
83 }
84
85 header.n_depth = 1;
86
87 in >> tag;
88 if (tag == "xllcorner" || tag == "xllcenter")
89 {
90 header.origin[0] = readDoubleFromStream(in);
91 }
92 else
93 {
94 return {};
95 }
96
97 in >> tag;
98 if (tag == "yllcorner" || tag == "yllcenter")
99 {
100 header.origin[1] = readDoubleFromStream(in);
101 }
102 else
103 {
104 return {};
105 }
106 header.origin[2] = 0;
107
108 in >> tag;
109 if (tag == "cellsize")
110 {
111 header.cell_size = readDoubleFromStream(in);
112 }
113 else
114 {
115 return {};
116 }
117
118 in >> tag;
119 if (tag == "NODATA_value" || tag == "nodata_value")
120 {
121 header.no_data = readDoubleFromStream(in);
122 }
123 else
124 {
125 return {};
126 }
127
128 return header;
129}
130
132 std::string const& fname)
133{
134 std::ifstream in(fname.c_str());
135
136 if (!in.is_open())
137 {
138 WARN("Raster::getRasterFromASCFile(): Could not open file {:s}.",
139 fname);
140 return nullptr;
141 }
142
143 auto const header = readASCHeader(in);
144 if (!header)
145 {
146 WARN(
147 "Raster::getRasterFromASCFile(): Could not read header of file "
148 "{:s}",
149 fname);
150 return nullptr;
151 }
152
153 std::vector<double> values(header->n_cols * header->n_rows);
154 // read the data into the double-array
155 for (std::size_t j(0); j < header->n_rows; ++j)
156 {
157 const std::size_t idx((header->n_rows - j - 1) * header->n_cols);
158 for (std::size_t i(0); i < header->n_cols; ++i)
159 {
160 values[idx + i] = readDoubleFromStream(in);
161 }
162 }
163
164 return new GeoLib::Raster(*header, values.begin(), values.end());
165}
166
169static std::optional<std::tuple<GeoLib::RasterHeader, double, double>>
170readSurferHeader(std::ifstream& in)
171{
172 std::string tag;
173
174 in >> tag;
175
176 if (tag != "DSAA")
177 {
178 ERR("Error in readSurferHeader() - No Surfer file.");
179 return {};
180 }
181
183 in >> header.n_cols >> header.n_rows;
184 double min, max;
185 in >> min >> max;
186 header.origin[0] = min;
187 header.cell_size = (max - min) / static_cast<double>(header.n_cols);
188
189 in >> min >> max;
190 header.origin[1] = min;
191 header.origin[2] = 0;
192
193 if (ceil((max - min) / static_cast<double>(header.n_rows)) ==
194 ceil(header.cell_size))
195 {
196 header.cell_size = ceil(header.cell_size);
197 }
198 else
199 {
200 ERR("Error in readSurferHeader() - Anisotropic cellsize detected.");
201 return {};
202 }
203 header.n_depth = 1;
204 header.no_data = -9999;
205 in >> min >> max;
206
207 return {{header, min, max}};
208}
209
211 std::string const& fname)
212{
213 std::ifstream in(fname.c_str());
214
215 if (!in.is_open())
216 {
217 ERR("Raster::getRasterFromSurferFile() - Could not open file {:s}",
218 fname);
219 return nullptr;
220 }
221
222 auto const optional_header = readSurferHeader(in);
223 if (!optional_header)
224 {
225 ERR("Raster::getRasterFromASCFile() - could not read header of file "
226 "{:s}",
227 fname);
228 return nullptr;
229 }
230
231 auto const [header, min, max] = *optional_header;
232 std::vector<double> values(header.n_cols * header.n_rows);
233 // read the data into the double-array
234 for (std::size_t j(0); j < header.n_rows; ++j)
235 {
236 const std::size_t idx(j * header.n_cols);
237 for (std::size_t i(0); i < header.n_cols; ++i)
238 {
239 const double val = readDoubleFromStream(in);
240 values[idx + i] = (val > max || val < min) ? header.no_data : val;
241 }
242 }
243
244 return new GeoLib::Raster(header, values.begin(), values.end());
245}
246
247std::optional<std::array<double, 3>> readCoordinates(std::istream& in)
248{
249 std::string line("");
250 if (std::getline(in, line))
251 {
252 std::stringstream str_stream(line);
253 std::array<double, 3> coords;
254 str_stream >> coords[0] >> coords[1] >> coords[2];
255 return std::make_optional(coords);
256 }
257 return std::nullopt;
258}
259
261 std::string const& fname)
262{
263 std::ifstream in(fname.c_str());
264 if (!in.is_open())
265 {
266 ERR("Raster::getRasterFromXyzFile() - Could not open file {:s}", fname);
267 return nullptr;
268 }
269
270 auto coords = readCoordinates(in);
271 if (coords == std::nullopt)
272 {
273 return nullptr;
274 }
275
276 std::vector<double> values;
277 values.push_back((*coords)[2]);
278
279 auto coords2 = readCoordinates(in);
280 if (coords2 == std::nullopt)
281 {
282 return nullptr;
283 }
284 values.push_back((*coords2)[2]);
286 0, 0, 1, GeoLib::Point(*coords), (*coords2)[0] - (*coords)[0], -9999};
287
288 std::size_t n_cols = 2, n_rows = 1;
289 while ((coords = readCoordinates(in)))
290 {
291 values.push_back((*coords)[2]);
292 if ((*coords)[0] > (*coords2)[0])
293 {
294 if ((*coords)[0] - (*coords2)[0] != header.cell_size)
295 {
296 ERR("Varying cell sizes or unordered pixel values found. "
297 "Aborting...");
298 return nullptr;
299 }
300 n_cols++;
301 }
302 else // new line
303 {
304 if ((*coords)[1] - (*coords2)[1] != header.cell_size)
305 {
306 ERR("Varying cell sizes or unordered pixel values found. "
307 "Aborting...");
308 return nullptr;
309 }
310 n_rows++;
311 // define #columns
312 if (header.n_cols == 0)
313 {
314 header.n_cols = n_cols;
315 }
316 // just check if #columns is consistent
317 else
318 {
319 if (n_cols != header.n_cols)
320 {
321 ERR("Different number of pixels per line. Aborting!");
322 return nullptr;
323 }
324 }
325 n_cols = 1;
326 }
327 coords2 = coords;
328 }
329 header.n_rows = n_rows;
330 if (header.n_cols == 0)
331 {
332 ERR("Could not determine raster size. Note that minimum allowed raster "
333 "size is 2 x 2 pixels.");
334 return nullptr;
335 }
336 return new GeoLib::Raster(header, values.begin(), values.end());
337}
338
340 std::string const& file_name)
341{
342 GeoLib::RasterHeader header(raster.getHeader());
343 MathLib::Point3d const& origin(header.origin);
344 unsigned const nCols(header.n_cols);
345 unsigned const nRows(header.n_rows);
346
347 // write header
348 std::ofstream out(file_name);
349 out << "ncols " << nCols << "\n";
350 out << "nrows " << nRows << "\n";
351 auto const default_precision = out.precision();
352 out.precision(std::numeric_limits<double>::digits10);
353 out << "xllcorner " << origin[0] << "\n";
354 out << "yllcorner " << origin[1] << "\n";
355 out << "cellsize " << header.cell_size << "\n";
356 out.precision(default_precision);
357 out << "NODATA_value " << header.no_data << "\n";
358
359 // write data
360 double const* const elevation(raster.begin());
361 for (unsigned row(0); row < nRows; ++row)
362 {
363 for (unsigned col(0); col < nCols - 1; ++col)
364 {
365 out << elevation[(nRows - row - 1) * nCols + col] << " ";
366 }
367 out << elevation[(nRows - row) * nCols - 1] << "\n";
368 }
369 out.close();
370}
371
373static bool allRastersExist(std::vector<std::string> const& raster_paths)
374{
375 for (const auto& raster_path : raster_paths)
376 {
377 std::ifstream file_stream(raster_path, std::ifstream::in);
378 if (!file_stream.good())
379 {
380 ERR("Opening raster file {} failed.", raster_path);
381 return false;
382 }
383 file_stream.close();
384 }
385 return true;
386}
387
388std::optional<std::vector<GeoLib::Raster const*>> readRasters(
389 std::vector<std::string> const& raster_paths)
390{
391 if (!allRastersExist(raster_paths))
392 {
393 return std::nullopt;
394 }
395
396 std::vector<GeoLib::Raster const*> rasters;
397 rasters.reserve(raster_paths.size());
398 std::transform(raster_paths.begin(), raster_paths.end(),
399 std::back_inserter(rasters),
400 [](auto const& path)
401 { return FileIO::AsciiRasterInterface::readRaster(path); });
402 return std::make_optional(rasters);
403}
404} // end namespace FileIO
Definition of the AsciiRasterInterface class.
Filename manipulation routines.
Definition of the Point class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:40
Definition of string helper functions.
static GeoLib::Raster * getRasterFromXyzFile(std::string const &fname)
Reads a XYZ raster file.
static GeoLib::Raster * getRasterFromSurferFile(std::string const &fname)
Reads a Surfer GRD raster file.
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.
static GeoLib::Raster * getRasterFromASCFile(std::string const &fname)
Reads an ArcGis ASC raster file.
Class Raster is used for managing raster data.
Definition: Raster.h:45
RasterHeader const & getHeader() const
Returns the complete header information.
Definition: Raster.h:83
const_iterator begin() const
Definition: Raster.h:94
std::string getFileExtension(const std::string &path)
Definition: FileTools.cpp:190
std::string replaceString(const std::string &searchString, const std::string &replaceString, std::string stringToReplace)
Definition: StringTools.cpp:49
static std::optional< GeoLib::RasterHeader > readASCHeader(std::ifstream &in)
static double readDoubleFromStream(std::istream &in)
Reads a double replacing comma by point.
std::optional< std::array< double, 3 > > readCoordinates(std::istream &in)
std::optional< std::vector< GeoLib::Raster const * > > readRasters(std::vector< std::string > const &raster_paths)
static std::optional< std::tuple< GeoLib::RasterHeader, double, double > > readSurferHeader(std::ifstream &in)
static bool allRastersExist(std::vector< std::string > const &raster_paths)
Checks if all raster files actually exist.
Contains the relevant information when storing a geoscientific raster data.
Definition: Raster.h:27
std::size_t n_depth
Definition: Raster.h:30
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