OGS
AsciiRasterInterface.cpp
Go to the documentation of this file.
1 
14 #include "AsciiRasterInterface.h"
15 
16 #include <tuple>
17 
18 #include "BaseLib/FileTools.h"
19 #include "BaseLib/Logging.h"
20 #include "BaseLib/StringTools.h"
21 
22 namespace FileIO
23 {
25 {
26  std::string ext(BaseLib::getFileExtension(fname));
27  std::transform(ext.begin(), ext.end(), ext.begin(), tolower);
28  if (ext == ".asc")
29  {
30  return getRasterFromASCFile(fname);
31  }
32  if (ext == ".grd")
33  {
34  return getRasterFromSurferFile(fname);
35  }
36  return nullptr;
37 }
38 
40 static double readDoubleFromStream(std::istream& in)
41 {
42  std::string value;
43  in >> value;
44  return std::strtod(BaseLib::replaceString(",", ".", value).c_str(),
45  nullptr);
46 }
47 
50 static std::optional<GeoLib::RasterHeader> readASCHeader(std::ifstream& in)
51 {
52  GeoLib::RasterHeader header;
53 
54  std::string tag;
55  std::string value;
56 
57  in >> tag;
58  if (tag == "ncols")
59  {
60  in >> value;
61  header.n_cols = atoi(value.c_str());
62  }
63  else
64  {
65  return {};
66  }
67 
68  in >> tag;
69  if (tag == "nrows")
70  {
71  in >> value;
72  header.n_rows = atoi(value.c_str());
73  }
74  else
75  {
76  return {};
77  }
78 
79  header.n_depth = 1;
80 
81  in >> tag;
82  if (tag == "xllcorner" || tag == "xllcenter")
83  {
84  header.origin[0] = readDoubleFromStream(in);
85  }
86  else
87  {
88  return {};
89  }
90 
91  in >> tag;
92  if (tag == "yllcorner" || tag == "yllcenter")
93  {
94  header.origin[1] = readDoubleFromStream(in);
95  }
96  else
97  {
98  return {};
99  }
100  header.origin[2] = 0;
101 
102  in >> tag;
103  if (tag == "cellsize")
104  {
105  header.cell_size = readDoubleFromStream(in);
106  }
107  else
108  {
109  return {};
110  }
111 
112  in >> tag;
113  if (tag == "NODATA_value" || tag == "nodata_value")
114  {
115  header.no_data = readDoubleFromStream(in);
116  }
117  else
118  {
119  return {};
120  }
121 
122  return header;
123 }
124 
126  std::string const& fname)
127 {
128  std::ifstream in(fname.c_str());
129 
130  if (!in.is_open())
131  {
132  WARN("Raster::getRasterFromASCFile(): Could not open file {:s}.",
133  fname);
134  return nullptr;
135  }
136 
137  auto const header = readASCHeader(in);
138  if (!header)
139  {
140  WARN(
141  "Raster::getRasterFromASCFile(): Could not read header of file "
142  "{:s}",
143  fname);
144  return nullptr;
145  }
146 
147  std::vector<double> values(header->n_cols * header->n_rows);
148  // read the data into the double-array
149  for (std::size_t j(0); j < header->n_rows; ++j)
150  {
151  const std::size_t idx((header->n_rows - j - 1) * header->n_cols);
152  for (std::size_t i(0); i < header->n_cols; ++i)
153  {
154  values[idx + i] = readDoubleFromStream(in);
155  }
156  }
157 
158  return new GeoLib::Raster(*header, values.begin(), values.end());
159 }
160 
163 static std::optional<std::tuple<GeoLib::RasterHeader, double, double>>
164 readSurferHeader(std::ifstream& in)
165 {
166  std::string tag;
167 
168  in >> tag;
169 
170  if (tag != "DSAA")
171  {
172  ERR("Error in readSurferHeader() - No Surfer file.");
173  return {};
174  }
175 
176  GeoLib::RasterHeader header;
177  in >> header.n_cols >> header.n_rows;
178  double min, max;
179  in >> min >> max;
180  header.origin[0] = min;
181  header.cell_size = (max - min) / static_cast<double>(header.n_cols);
182 
183  in >> min >> max;
184  header.origin[1] = min;
185  header.origin[2] = 0;
186 
187  if (ceil((max - min) / static_cast<double>(header.n_rows)) ==
188  ceil(header.cell_size))
189  {
190  header.cell_size = ceil(header.cell_size);
191  }
192  else
193  {
194  ERR("Error in readSurferHeader() - Anisotropic cellsize detected.");
195  return {};
196  }
197  header.n_depth = 1;
198  header.no_data = min - 1;
199  in >> min >> max;
200 
201  return {{header, min, max}};
202 }
203 
205  std::string const& fname)
206 {
207  std::ifstream in(fname.c_str());
208 
209  if (!in.is_open())
210  {
211  ERR("Raster::getRasterFromSurferFile() - Could not open file {:s}",
212  fname);
213  return nullptr;
214  }
215 
216  auto const optional_header = readSurferHeader(in);
217  if (!optional_header)
218  {
219  ERR("Raster::getRasterFromASCFile() - could not read header of file "
220  "{:s}",
221  fname);
222  return nullptr;
223  }
224 
225  auto const [header, min, max] = *optional_header;
226  const double no_data_val(min - 1);
227  std::vector<double> values(header.n_cols * header.n_rows);
228  // read the data into the double-array
229  for (std::size_t j(0); j < header.n_rows; ++j)
230  {
231  const std::size_t idx(j * header.n_cols);
232  for (std::size_t i(0); i < header.n_cols; ++i)
233  {
234  const double val = readDoubleFromStream(in);
235  values[idx + i] = (val > max || val < min) ? no_data_val : val;
236  }
237  }
238 
239  return new GeoLib::Raster(header, values.begin(), values.end());
240 }
241 
243  std::string const& file_name)
244 {
245  GeoLib::RasterHeader header(raster.getHeader());
246  MathLib::Point3d const& origin(header.origin);
247  unsigned const nCols(header.n_cols);
248  unsigned const nRows(header.n_rows);
249 
250  // write header
251  std::ofstream out(file_name);
252  out << "ncols " << nCols << "\n";
253  out << "nrows " << nRows << "\n";
254  auto const default_precision = out.precision();
255  out.precision(std::numeric_limits<double>::digits10);
256  out << "xllcorner " << origin[0] << "\n";
257  out << "yllcorner " << origin[1] << "\n";
258  out << "cellsize " << header.cell_size << "\n";
259  out.precision(default_precision);
260  out << "NODATA_value " << header.no_data << "\n";
261 
262  // write data
263  double const* const elevation(raster.begin());
264  for (unsigned row(0); row < nRows; ++row)
265  {
266  for (unsigned col(0); col < nCols - 1; ++col)
267  {
268  out << elevation[(nRows - row - 1) * nCols + col] << " ";
269  }
270  out << elevation[(nRows - row) * nCols - 1] << "\n";
271  }
272  out.close();
273 }
274 
276 static bool allRastersExist(std::vector<std::string> const& raster_paths)
277 {
278  for (const auto& raster_path : raster_paths)
279  {
280  std::ifstream file_stream(raster_path, std::ifstream::in);
281  if (!file_stream.good())
282  {
283  ERR("Opening raster file {} failed.", raster_path);
284  return false;
285  }
286  file_stream.close();
287  }
288  return true;
289 }
290 
291 std::optional<std::vector<GeoLib::Raster const*>> readRasters(
292  std::vector<std::string> const& raster_paths)
293 {
294  if (!allRastersExist(raster_paths))
295  {
296  return std::nullopt;
297  }
298 
299  std::vector<GeoLib::Raster const*> rasters;
300  rasters.reserve(raster_paths.size());
301  std::transform(raster_paths.begin(), raster_paths.end(),
302  std::back_inserter(rasters),
303  [](auto const& path)
304  { return FileIO::AsciiRasterInterface::readRaster(path); });
305  return std::make_optional(rasters);
306 }
307 } // end namespace FileIO
Definition of the AsciiRasterInterface class.
Filename manipulation routines.
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Definition of string helper functions.
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
const_iterator begin() const
Definition: Raster.h:94
RasterHeader const & getHeader() const
Returns the complete header information.
Definition: Raster.h:83
std::string getFileExtension(const std::string &path)
Definition: FileTools.cpp:186
std::string replaceString(const std::string &searchString, const std::string &replaceString, std::string stringToReplace)
Definition: StringTools.cpp:50
static std::optional< std::tuple< GeoLib::RasterHeader, double, double > > readSurferHeader(std::ifstream &in)
static double readDoubleFromStream(std::istream &in)
Reads a double replacing comma by point.
std::optional< std::vector< GeoLib::Raster const * > > readRasters(std::vector< std::string > const &raster_paths)
static std::optional< GeoLib::RasterHeader > readASCHeader(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