ryujin 2.1.1 revision dbf0e3ba7acdb60b6d558e4257815df4a8f8daf9
geotiff_reader.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2024 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include "convenience_macros.h"
11#include "lazy.h"
12#include "patterns_conversion.h"
13
14#include <deal.II/base/parameter_acceptor.h>
15
16#include <numeric>
17
18#ifdef WITH_GDAL
19#include <cpl_conv.h>
20#include <gdal.h>
21#include <gdal_priv.h>
22#endif
23
24namespace ryujin
25{
28 none,
30 minimum,
32 average,
35 };
36} // namespace ryujin
37
38#ifndef DOXYGEN
44#endif
45
46namespace ryujin
47{
55 class GeoTIFFReader : dealii::ParameterAcceptor
56 {
57 public:
58 GeoTIFFReader(const std::string subsection)
59 : ParameterAcceptor(subsection)
60 {
61 filename_ = "ryujin.tif";
62 this->add_parameter("filename", filename_, "GeoTIFF: image file to load");
63
64 transformation_ = {0., 0.01, 0., 0., 0., 0.01};
65 this->add_parameter(
66 "transformation",
67 transformation_,
68 "Array \"t[]\" describing an affine transformation between image "
69 "space (indices i and j from bottom left) and real coordinates (x "
70 "and y): x = t[0] + t[1] * i + t[2] * j, and y = t[3] + t[4] * i + "
71 "t[5] * j. (This transformation sets the origin of the image space "
72 "into the bottom left corner with index i to the right and index j "
73 "up)");
74
75 transformation_allow_out_of_bounds_queries_ = false;
76 this->add_parameter(
77 "transformation allow out of bounds queries",
78 transformation_allow_out_of_bounds_queries_,
79 "GeoTIFF: allow out-of-bounds queries. When set to true, the reader "
80 "returns constant extended values for coordinates that are outside "
81 "of the image range.");
82
83 transformation_use_geotiff_ = true;
84 this->add_parameter("transformation use geotiff",
85 transformation_use_geotiff_,
86 "GeoTIFF: read in transformation from GeoTIFF for "
87 "constructing the affine transformation. If set to "
88 "false the manually specified transformation "
89 "parameters will be used instead.");
90
91 transformation_use_geotiff_origin_ = false;
92 this->add_parameter(
93 "transformation use geotiff origin",
94 transformation_use_geotiff_origin_,
95 "GeoTIFF: read in affine shift (i.e., position of "
96 "lower left corner) from GeoTIFF for constructing "
97 "the affine transformation. If set to false the origin specified "
98 "in the transformation parameter will be used instead.");
99
100 height_normalization_ = HeightNormalization::minimum;
101 this->add_parameter("height normalization",
102 height_normalization_,
103 "GeoTIFF: choose base point for height normalization "
104 "that is set to 0.: none, minimum, average, maximum");
105
106 height_scaling_ = 1.0;
107 this->add_parameter("height scaling",
108 height_scaling_,
109 "GeoTIFF: choose base point for height normalization "
110 "that is set to 0.: none, minimum, average, maximum");
111
112 const auto set_up = [&] {
113#ifdef WITH_GDAL
114 /* Initial GDAL and reset all data: */
115 GDALAllRegister();
116 driver_name_ = "";
117 driver_projection_ = "";
118 affine_transformation_ = {0, 0, 0, 0, 0, 0};
119 inverse_affine_transformation_ = {0, 0, 0, 0, 0, 0};
120 raster_offset_ = {0, 0};
121 raster_size_ = {0, 0};
122 raster_.clear();
123#endif
124 };
125
126 set_up();
127 this->parse_parameters_call_back.connect(set_up);
128 }
129
130 /*
131 * Query height information for a given position. The method is
132 * templated in dim, so that it can be called with a @p point having an
133 * arbitrary dimension. The method, however, only only uses the first
134 * two coordinates, x and y, from @p point, translates (x, y) into
135 * image coordinates (i, j) and returns an interpolated height value.
136 *
137 * @note If dim < 2, then the y coordinate is set to zero. Similarly,
138 * if dim == 0, then the tiff raster is queried at coordinate (0, 0).
139 */
140 template <int dim>
141 DEAL_II_ALWAYS_INLINE inline double
142 compute_height(const dealii::Point<dim> &point) const
143 {
144 geotiff_guard_.ensure_initialized([&]() {
145 read_in_raster();
146 return true;
147 });
148
149 double x = 0;
150 if constexpr (dim >= 1)
151 x = point[0];
152
153 double y = 0;
154 if constexpr (dim >= 2)
155 y = point[1];
156
157 const auto &[di, dj] = apply_inverse_transformation(x, y);
158
159 /* Check that we are in bounds: */
160
161 const bool in_bounds =
162 di > -1.0 && di < static_cast<double>(raster_size_[0]) + 1.0 &&
163 dj > -1.0 && dj < static_cast<double>(raster_size_[1]) + 1.0;
164
165#ifdef DEBUG_OUTPUT
166 if (!in_bounds) {
167 std::cout << std::setprecision(16);
168 std::cout << "Queried point out of bounds." << std::endl;
169 std::cout << "Point: " << point << std::endl;
170 std::cout << "Transformed coordinates: (" << di << "," << dj << ")"
171 << std::endl;
172 }
173#endif
174
175 AssertThrow(
176 transformation_allow_out_of_bounds_queries_ || in_bounds,
177 dealii::ExcMessage("Raster error: The requested point is outside "
178 "the image boundary of the geotiff file"));
179
180 /*
181 * Use a simple bilinear interpolation and ensure we never go below
182 * the minimum or above the maximum index.
183 */
184
185 const auto i_left = std::min(
186 std::max(static_cast<int>(std::floor(di)), 0), raster_size_[0]);
187 const auto i_right = std::min(
188 std::max(static_cast<int>(std::ceil(di)), 0), raster_size_[0]);
189 const auto j_left = std::min(
190 std::max(static_cast<int>(std::floor(dj)), 0), raster_size_[1]);
191 const auto j_right = std::min(
192 std::max(static_cast<int>(std::ceil(dj)), 0), raster_size_[1]);
193
194#ifdef DEBUG_OUTPUT
195 if (!in_bounds) {
196 std::cout << "index bounding box: (" << i_left << "," << j_left
197 << ") and (" << i_right << "," << j_right << ")" << std::endl;
198 }
199#endif
200
201 const double i_ratio = std::fmod(di, 1.);
202 const double j_ratio = std::fmod(dj, 1.);
203
204 const auto v_iljl = raster_[i_left + j_left * raster_size_[0]];
205 const auto v_irjl = raster_[i_right + j_left * raster_size_[0]];
206
207 const auto v_iljr = raster_[i_left + j_right * raster_size_[0]];
208 const auto v_irjr = raster_[i_right + j_right * raster_size_[0]];
209
210 const auto v_jl = v_iljl * (1. - i_ratio) + v_irjl * i_ratio;
211 const auto v_jr = v_iljr * (1. - i_ratio) + v_irjr * i_ratio;
212
213 return height_scaling_ * (v_jl * (1. - j_ratio) + v_jr * j_ratio);
214 }
215
216 /*
217 * Return the affine transformation information that is stored in the
218 * GeoTIFF image.
219 */
221
225
226 private:
227 void read_in_raster() const
228 {
229#ifdef WITH_GDAL
230 auto dataset_handle = GDALOpen(filename_.c_str(), GA_ReadOnly);
231 AssertThrow(dataset_handle,
232 dealii::ExcMessage("GDAL error: file not found"));
233
234 auto dataset = GDALDataset::FromHandle(dataset_handle);
235 Assert(dataset, dealii::ExcInternalError());
236
237 const auto driver = dataset->GetDriver();
238
239 driver_name_ = driver->GetMetadataItem(GDAL_DMD_LONGNAME);
240 if (dataset->GetProjectionRef() != nullptr)
241 driver_projection_ = dataset->GetProjectionRef();
242
243 /* For now we support one raster in the dataset: */
244
245 AssertThrow(
246 dataset->GetRasterCount() == 1,
247 dealii::ExcMessage(
248 "GDAL driver error: currently we only support one raster"));
249
250 const auto raster_band = dataset->GetRasterBand(1);
251
252 AssertThrow(dataset->GetRasterXSize() == raster_band->GetXSize() &&
253 dataset->GetRasterYSize() == raster_band->GetYSize(),
254 dealii::ExcMessage(
255 "GDAL driver error: the raster band has a different "
256 "dimension than the (global) raster dimension of the "
257 "geotiff image. This is not supported."));
258
259 /*
260 * FIXME: For now, we simply read in the entire geotiff file on
261 * each rank. In order to save memory for very large files it would
262 * be possible to create a bounding box for the all active cells of
263 * the triangulation and then only read in a small region for which
264 * we actually need data.
265 */
266
267 raster_offset_ = {0, 0};
268 raster_size_ = {dataset->GetRasterXSize(), dataset->GetRasterYSize()};
269
270 raster_.resize(raster_size_[0] * raster_size_[1]);
271 const auto error_code = raster_band->RasterIO(
272 GF_Read,
273 raster_offset_[0], /* x-offset of image region */
274 raster_offset_[1], /* y-offset of image region */
275 raster_size_[0], /* x-size of image region */
276 raster_size_[1], /* y-size of image region */
277 raster_.data(),
278 raster_size_[0], /* x-size of target buffer */
279 raster_size_[1], /* y-size of target buffer */
280 GDT_Float32,
281 0,
282 0);
283
284 AssertThrow(error_code == 0,
285 dealii::ExcMessage(
286 "GDAL driver error: error reading in geotiff file"));
287
288 /*
289 * Read in the affine transformation from the geotiff image.
290 *
291 * Note that this transformation differs from the one we use in the
292 * parameter file: GDAL uses typical image orientation: the origin
293 * of the dataset is in the "top left" corner (instead of bottom
294 * left) and the first (column) index goes to the right and the
295 * second (row) index goes down.
296 */
297
298 if (transformation_use_geotiff_) {
299 const auto success =
300 dataset->GetGeoTransform(affine_transformation_.data()) == CE_None;
301 AssertThrow(success,
302 dealii::ExcMessage("GDAL driver error: no geo transform "
303 "present in geotiff file"));
304 } else {
305 affine_transformation_ = transformation_;
306 /* Flip sign for j index (y-coordinate): */
307 affine_transformation_[2] *= -1.;
308 affine_transformation_[5] *= -1.;
309 }
310
311 /*
312 * Ensure that (i=0, j=raster_size[1]-1) corresponds to the user
313 * supplied (transformation_[0], transformation_[3]).
314 */
315 if (transformation_use_geotiff_ == false ||
316 transformation_use_geotiff_origin_ == false) {
317 const auto j_max = raster_size_[1] - 1;
318 affine_transformation_[0] =
319 transformation_[0] - j_max * affine_transformation_[2];
320 affine_transformation_[3] =
321 transformation_[3] - j_max * affine_transformation_[5];
322 }
323
324 /*
325 * Compute inverse transformation of
326 *
327 * x = t[0] + t[1] * i + t[2] * j, y = t[3] + t[4] * i + t[5] * j.
328 *
329 * namely:
330 *
331 * i = it[1] * (x - it[0]) + it[2] * (y - it[3])
332 * j = it[4] * (x - it[0]) + it[5] * (y - it[3])
333 */
334 inverse_affine_transformation_[0] = affine_transformation_[0];
335 inverse_affine_transformation_[3] = affine_transformation_[3];
336
337 const auto determinant =
338 affine_transformation_[1] * affine_transformation_[5] -
339 affine_transformation_[2] * affine_transformation_[4];
340 const auto inv = 1. / determinant;
341 inverse_affine_transformation_[1] = inv * affine_transformation_[5];
342 inverse_affine_transformation_[2] = inv * (-affine_transformation_[2]);
343 inverse_affine_transformation_[4] = inv * (-affine_transformation_[4]);
344 inverse_affine_transformation_[5] = inv * affine_transformation_[1];
345
346 GDALClose(dataset_handle);
347
348#ifdef DEBUG_OUTPUT
349 std::cout << std::setprecision(16);
350 std::cout << "GDAL: driver name = " << driver_name_;
351 std::cout << "\nGDAL: projection = " << driver_projection_;
352 std::cout << "\nGDAL: transformation =";
353 for (const auto &it : affine_transformation_)
354 std::cout << " " << it;
355 std::cout << "\nGDAL: inverse trafo =";
356 for (const auto &it : inverse_affine_transformation_)
357 std::cout << " " << it;
358 std::cout << "\nGDAL: raster offset =";
359 for (const auto &it : raster_offset_)
360 std::cout << " " << it;
361 std::cout << "\nGDAL: raster size =";
362 for (const auto &it : raster_size_)
363 std::cout << " " << it;
364 std::cout << std::endl;
365#endif
366
367 if (height_normalization_ != HeightNormalization::none) {
368 float shift = 0.;
369
370 if (height_normalization_ == HeightNormalization::minimum)
371 shift = *std::min_element(std::begin(raster_), std::end(raster_));
372 else if (height_normalization_ == HeightNormalization::maximum)
373 shift = *std::max_element(std::begin(raster_), std::end(raster_));
374 else {
375 Assert(height_normalization_ == HeightNormalization::average,
376 dealii::ExcInternalError());
377 const auto sum = std::reduce(std::begin(raster_), std::end(raster_));
378 shift = sum / raster_.size();
379 }
380
381 std::for_each(std::begin(raster_),
382 std::end(raster_),
383 [&](auto &element) { element -= shift; });
384 }
385
386#else
387 static constexpr auto message =
388 "ryujin has to be configured with GDAL support in order to read in "
389 "GeoTIFF images";
390 AssertThrow(false, dealii::ExcMessage(message));
391 __builtin_trap();
392#endif
393 }
394
395
396 DEAL_II_ALWAYS_INLINE inline std::array<double, 2>
397 apply_transformation(const double i, const double j) const
398 {
399 const auto &at = affine_transformation_;
400 const double x = at[0] + at[1] * i + at[2] * j;
401 const double y = at[3] + at[4] * i + at[5] * j;
402 return {x, y};
403 }
404
405
406 DEAL_II_ALWAYS_INLINE inline std::array<double, 2>
407 apply_inverse_transformation(const double x, const double y) const
408 {
409 const auto &iat = inverse_affine_transformation_;
410 const double i = iat[1] * (x - iat[0]) + iat[2] * (y - iat[3]);
411 const double j = iat[4] * (x - iat[0]) + iat[5] * (y - iat[3]);
412 return {i, j};
413 }
414
415 /* Runtime parameters: */
416
417 std::string filename_;
418
419 std::array<double, 6> transformation_;
420 bool transformation_allow_out_of_bounds_queries_;
421 bool transformation_use_geotiff_;
422 bool transformation_use_geotiff_origin_;
423 HeightNormalization height_normalization_;
424 double height_scaling_;
425
426 /* GDAL data structures: */
427
428 //
429 // We use a Lazy<t> wrapper for lazy initialization with efficient
430 // Schmidt's double checking. We simply ignore the bool type here.
431 //
432 Lazy<bool> geotiff_guard_;
433 mutable std::string driver_name_;
434 mutable std::string driver_projection_;
435 mutable std::array<double, 6> affine_transformation_;
436 mutable std::array<double, 6> inverse_affine_transformation_;
437 mutable std::array<int, 2> raster_offset_;
438 mutable std::array<int, 2> raster_size_;
439 mutable std::vector<float> raster_;
440 };
441} // namespace ryujin
const auto & raster_offset() const
const auto & raster_size() const
const auto & height_scaling() const
const auto & affine_transformation() const
GeoTIFFReader(const std::string subsection)
DEAL_II_ALWAYS_INLINE double compute_height(const dealii::Point< dim > &point) const
void ensure_initialized(const Callable &creator) const
#define ACCESSOR_READ_ONLY(member)
HeightNormalization