8#include <compile_time_options.h>
14#include <deal.II/base/parameter_acceptor.h>
59 : ParameterAcceptor(subsection)
61 filename_ =
"ryujin.tif";
62 this->add_parameter(
"filename", filename_,
"GeoTIFF: image file to load");
64 transformation_ = {0., 0.01, 0., 0., 0., 0.01};
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 "
75 transformation_allow_out_of_bounds_queries_ =
false;
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.");
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.");
91 transformation_use_geotiff_origin_ =
false;
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.");
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");
106 height_scaling_ = 1.0;
107 this->add_parameter(
"height scaling",
109 "GeoTIFF: choose base point for height normalization "
110 "that is set to 0.: none, minimum, average, maximum");
112 const auto set_up = [&] {
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};
127 this->parse_parameters_call_back.connect(set_up);
141 DEAL_II_ALWAYS_INLINE
inline double
150 if constexpr (dim >= 1)
154 if constexpr (dim >= 2)
157 const auto &[di, dj] = apply_inverse_transformation(x, y);
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;
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 <<
")"
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"));
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]);
196 std::cout <<
"index bounding box: (" << i_left <<
"," << j_left
197 <<
") and (" << i_right <<
"," << j_right <<
")" << std::endl;
201 const double i_ratio = std::fmod(di, 1.);
202 const double j_ratio = std::fmod(dj, 1.);
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]];
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]];
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;
213 return height_scaling_ * (v_jl * (1. - j_ratio) + v_jr * j_ratio);
227 void read_in_raster()
const
230 auto dataset_handle = GDALOpen(filename_.c_str(), GA_ReadOnly);
231 AssertThrow(dataset_handle,
232 dealii::ExcMessage(
"GDAL error: file not found"));
234 auto dataset = GDALDataset::FromHandle(dataset_handle);
235 Assert(dataset, dealii::ExcInternalError());
237 const auto driver = dataset->GetDriver();
239 driver_name_ = driver->GetMetadataItem(GDAL_DMD_LONGNAME);
240 if (dataset->GetProjectionRef() !=
nullptr)
241 driver_projection_ = dataset->GetProjectionRef();
246 dataset->GetRasterCount() == 1,
248 "GDAL driver error: currently we only support one raster"));
250 const auto raster_band = dataset->GetRasterBand(1);
252 AssertThrow(dataset->GetRasterXSize() == raster_band->GetXSize() &&
253 dataset->GetRasterYSize() == raster_band->GetYSize(),
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."));
267 raster_offset_ = {0, 0};
268 raster_size_ = {dataset->GetRasterXSize(), dataset->GetRasterYSize()};
270 raster_.resize(raster_size_[0] * raster_size_[1]);
271 const auto error_code = raster_band->RasterIO(
284 AssertThrow(error_code == 0,
286 "GDAL driver error: error reading in geotiff file"));
298 if (transformation_use_geotiff_) {
300 dataset->GetGeoTransform(affine_transformation_.data()) == CE_None;
302 dealii::ExcMessage(
"GDAL driver error: no geo transform "
303 "present in geotiff file"));
305 affine_transformation_ = transformation_;
307 affine_transformation_[2] *= -1.;
308 affine_transformation_[5] *= -1.;
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];
334 inverse_affine_transformation_[0] = affine_transformation_[0];
335 inverse_affine_transformation_[3] = affine_transformation_[3];
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];
346 GDALClose(dataset_handle);
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;
371 shift = *std::min_element(std::begin(raster_), std::end(raster_));
373 shift = *std::max_element(std::begin(raster_), std::end(raster_));
376 dealii::ExcInternalError());
377 const auto sum = std::reduce(std::begin(raster_), std::end(raster_));
378 shift = sum / raster_.size();
381 std::for_each(std::begin(raster_),
383 [&](
auto &element) { element -= shift; });
387 static constexpr auto message =
388 "ryujin has to be configured with GDAL support in order to read in "
390 AssertThrow(
false, dealii::ExcMessage(message));
396 DEAL_II_ALWAYS_INLINE
inline std::array<double, 2>
397 apply_transformation(
const double i,
const double j)
const
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;
406 DEAL_II_ALWAYS_INLINE
inline std::array<double, 2>
407 apply_inverse_transformation(
const double x,
const double y)
const
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]);
417 std::string filename_;
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_;
424 double height_scaling_;
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_;
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)