OmniSciDB  91042dcc5b
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GeoRasterTableFunctions.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2021 MapD Technologies, Inc.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 #ifndef __CUDACC__
18 
19 #include <cmath>
20 #include <vector>
21 
22 #include <tbb/parallel_for.h>
23 #include <tbb/task_arena.h>
24 #include <tbb/task_group.h>
26 #include "Shared/math_consts.h"
27 
28 const size_t max_inputs_per_thread = 1000000L;
29 const size_t max_temp_output_entries = 200000000L;
30 
31 // Allow input types to GeoRaster that are different than class types/output Z type
32 // So we can move everything to the type of T and Z (which can each be either float
33 // or double)
34 template <typename T, typename Z>
35 template <typename T2, typename Z2>
37  const Column<T2>& input_y,
38  const Column<Z2>& input_z,
39  const double bin_dim_meters,
40  const bool geographic_coords,
41  const bool align_bins_to_zero_based_grid)
42  : bin_dim_meters_(bin_dim_meters)
43  , geographic_coords_(geographic_coords)
44  , null_sentinel_(std::numeric_limits<Z>::lowest()) {
45  auto timer = DEBUG_TIMER(__func__);
46  const int64_t input_size{input_z.size()};
47  if (input_size <= 0) {
48  num_bins_ = 0;
49  num_x_bins_ = 0;
50  num_y_bins_ = 0;
51  return;
52  }
53  const auto min_max_x = get_column_min_max(input_x);
54  const auto min_max_y = get_column_min_max(input_y);
55  x_min_ = min_max_x.first;
56  x_max_ = min_max_x.second;
57  y_min_ = min_max_y.first;
58  y_max_ = min_max_y.second;
59 
60  if (align_bins_to_zero_based_grid && !geographic_coords_) {
61  // For implicit, data-defined bounds, we treat the max of the x and y ranges as
62  // inclusive (closed interval), since if the max of the data in either x/y dimensions
63  // is at the first value of the next bin, values at that max will be discarded if we
64  // don't include the final bin. For exmaple, if the input data (perhaps already binned
65  // with a group by query) goes from 0.0 to 40.0 in both x and y directions, we should
66  // have the last x/y bins cover the range [40.0, 50.0), not [30.0, 40.0)
68  }
69 
71  computeParallel(input_x, input_y, input_z, max_inputs_per_thread);
72 }
73 
74 // Allow input types to GeoRaster that are different than class types/output Z type
75 // So we can move everything to the type of T and Z (which can each be either float
76 // or double)
77 template <typename T, typename Z>
78 template <typename T2, typename Z2>
80  const Column<T2>& input_y,
81  const Column<Z2>& input_z,
82  const double bin_dim_meters,
83  const bool geographic_coords,
84  const bool align_bins_to_zero_based_grid,
85  const T x_min,
86  const T x_max,
87  const T y_min,
88  const T y_max)
89  : bin_dim_meters_(bin_dim_meters)
90  , geographic_coords_(geographic_coords)
91  , null_sentinel_(std::numeric_limits<Z>::lowest())
92  , x_min_(x_min)
93  , x_max_(x_max)
94  , y_min_(y_min)
95  , y_max_(y_max) {
96  auto timer = DEBUG_TIMER(__func__);
97  if (align_bins_to_zero_based_grid && !geographic_coords_) {
98  // For explicit, user-defined bounds, we treat the max of the x and y ranges as
99  // exclusive (open interval), since if the user specifies the max x/y as the end of
100  // the bin, they do not intend to add the next full bin For example, if a user
101  // specifies a bin_dim_meters of 10.0 and an x and y range from 0 to 40.0, they almost
102  // assuredly intend for there to be 4 bins in each of the x and y dimensions, with the
103  // last bin of range [30.0, 40.0), not 5 with the final bin's range from [40.0, 50.0)
105  }
107  computeParallel(input_x, input_y, input_z, max_inputs_per_thread);
108 }
109 
110 template <typename T, typename Z>
111 inline Z GeoRaster<T, Z>::offset_source_z_from_raster_z(const int64_t source_x_bin,
112  const int64_t source_y_bin,
113  const Z source_z_offset) const {
114  if (is_bin_out_of_bounds(source_x_bin, source_y_bin)) {
115  return null_sentinel_;
116  }
117  const Z terrain_z = z_[x_y_bin_to_bin_index(source_x_bin, source_y_bin, num_x_bins_)];
118  if (terrain_z == null_sentinel_) {
119  return terrain_z;
120  }
121  return terrain_z + source_z_offset;
122 }
123 
124 template <typename T, typename Z>
125 inline Z GeoRaster<T, Z>::fill_bin_from_avg_neighbors(const int64_t x_centroid_bin,
126  const int64_t y_centroid_bin,
127  const int64_t bins_radius) const {
128  T val = 0.0;
129  int32_t count = 0;
130  for (int64_t y_bin = y_centroid_bin - bins_radius;
131  y_bin <= y_centroid_bin + bins_radius;
132  y_bin++) {
133  for (int64_t x_bin = x_centroid_bin - bins_radius;
134  x_bin <= x_centroid_bin + bins_radius;
135  x_bin++) {
136  if (x_bin >= 0 && x_bin < num_x_bins_ && y_bin >= 0 && y_bin < num_y_bins_) {
137  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
138  const Z bin_val = z_[bin_idx];
139  if (bin_val != null_sentinel_) {
140  count++;
141  val += bin_val;
142  }
143  }
144  }
145  }
146  return (count == 0) ? null_sentinel_ : val / count;
147 }
148 
149 template <typename T, typename Z>
151  x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
152  x_max_ = std::floor(x_max_ / bin_dim_meters_) * bin_dim_meters_ +
153  bin_dim_meters_; // Snap to end of bin
154  y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
155  y_max_ = std::floor(y_max_ / bin_dim_meters_) * bin_dim_meters_ +
156  bin_dim_meters_; // Snap to end of bin
157 }
158 
159 template <typename T, typename Z>
161  x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
162  x_max_ = std::ceil(x_max_ / bin_dim_meters_) * bin_dim_meters_;
163  y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
164  y_max_ = std::ceil(y_max_ / bin_dim_meters_) * bin_dim_meters_;
165 }
166 
167 template <typename T, typename Z>
169  x_range_ = x_max_ - x_min_;
170  y_range_ = y_max_ - y_min_;
171  if (geographic_coords_) {
172  const T x_centroid = (x_min_ + x_max_) * 0.5;
173  const T y_centroid = (y_min_ + y_max_) * 0.5;
174  x_meters_per_degree_ =
175  distance_in_meters(x_min_, y_centroid, x_max_, y_centroid) / x_range_;
176 
177  y_meters_per_degree_ =
178  distance_in_meters(x_centroid, y_min_, x_centroid, y_max_) / y_range_;
179 
180  num_x_bins_ = x_range_ * x_meters_per_degree_ / bin_dim_meters_;
181  num_y_bins_ = y_range_ * y_meters_per_degree_ / bin_dim_meters_;
182 
183  x_scale_input_to_bin_ = x_meters_per_degree_ / bin_dim_meters_;
184  y_scale_input_to_bin_ = y_meters_per_degree_ / bin_dim_meters_;
185  x_scale_bin_to_input_ = bin_dim_meters_ / x_meters_per_degree_;
186  y_scale_bin_to_input_ = bin_dim_meters_ / y_meters_per_degree_;
187 
188  } else {
189  num_x_bins_ = x_range_ / bin_dim_meters_;
190  num_y_bins_ = y_range_ / bin_dim_meters_;
191 
192  x_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
193  y_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
194  x_scale_bin_to_input_ = bin_dim_meters_;
195  y_scale_bin_to_input_ = bin_dim_meters_;
196  }
197  num_bins_ = num_x_bins_ * num_y_bins_;
198 }
199 
200 template <typename T, typename Z>
201 template <typename T2, typename Z2>
203  const Column<T2>& input_y,
204  const Column<Z2>& input_z) {
205  auto timer = DEBUG_TIMER(__func__);
206  const int64_t input_size{input_z.size()};
207  z_.resize(num_bins_, null_sentinel_);
208  for (int64_t sparse_idx = 0; sparse_idx != input_size; ++sparse_idx) {
209  const int64_t x_bin = get_x_bin(input_x[sparse_idx]);
210  const int64_t y_bin = get_y_bin(input_y[sparse_idx]);
211  if (x_bin < 0 || x_bin >= num_x_bins_ || y_bin < 0 || y_bin >= num_y_bins_) {
212  continue;
213  }
214  // Take the max height for this version, but may want to allow different metrics
215  // like average as well
216  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
217  if (!(input_z.isNull(sparse_idx)) && input_z[sparse_idx] > z_[bin_idx]) {
218  z_[bin_idx] = input_z[sparse_idx];
219  }
220  }
221 }
222 
223 template <typename T, typename Z>
224 template <typename T2, typename Z2>
226  const Column<T2>& input_y,
227  const Column<Z2>& input_z,
228  const size_t max_inputs_per_thread) {
229  const size_t input_size = input_z.size();
230  const size_t max_thread_count = std::thread::hardware_concurrency();
231  const size_t num_threads_by_input_elements =
232  std::min(max_thread_count,
233  ((input_size + max_inputs_per_thread - 1) / max_inputs_per_thread));
234  const size_t num_threads_by_output_size =
235  std::min(max_thread_count, ((max_temp_output_entries + num_bins_ - 1) / num_bins_));
236  const size_t num_threads =
237  std::min(num_threads_by_input_elements, num_threads_by_output_size);
238  if (num_threads <= 1) {
239  compute(input_x, input_y, input_z);
240  return;
241  }
242  auto timer = DEBUG_TIMER(__func__);
243 
244  std::vector<std::vector<Z>> per_thread_z_outputs(num_threads);
245 
246  tbb::parallel_for(tbb::blocked_range<size_t>(0, num_threads),
247  [&](const tbb::blocked_range<size_t>& r) {
248  for (size_t t = r.begin(); t != r.end(); ++t) {
249  per_thread_z_outputs[t].resize(num_bins_, null_sentinel_);
250  }
251  });
252 
253  tbb::task_arena limited_arena(num_threads);
254  tbb::task_group tg;
255  limited_arena.execute([&] {
256  tg.run([&] {
258  tbb::blocked_range<int64_t>(0, input_size),
259  [&](const tbb::blocked_range<int64_t>& r) {
260  const int64_t start_idx = r.begin();
261  const int64_t end_idx = r.end();
262  size_t thread_idx = tbb::this_task_arena::current_thread_index();
263  std::vector<Z>& this_thread_z_output = per_thread_z_outputs[thread_idx];
264 
265  for (int64_t sparse_idx = start_idx; sparse_idx != end_idx; ++sparse_idx) {
266  const int64_t x_bin = get_x_bin(input_x[sparse_idx]);
267  const int64_t y_bin = get_y_bin(input_y[sparse_idx]);
268  if (x_bin < 0 || x_bin >= num_x_bins_ || y_bin < 0 ||
269  y_bin >= num_y_bins_) {
270  continue;
271  }
272  // Take the max height for this version, but may want to allow different
273  // metrics like average as well
274  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
275  if (!(input_z.isNull(sparse_idx)) &&
276  input_z[sparse_idx] > this_thread_z_output[bin_idx]) {
277  this_thread_z_output[bin_idx] = input_z[sparse_idx];
278  }
279  }
280  });
281  });
282  });
283 
284  limited_arena.execute([&] { tg.wait(); });
285 
286  z_.resize(num_bins_, null_sentinel_);
287 
289  tbb::blocked_range<size_t>(0, num_bins_), [&](const tbb::blocked_range<size_t>& r) {
290  const size_t start_idx = r.begin();
291  const size_t end_idx = r.end();
292  for (size_t bin_idx = start_idx; bin_idx != end_idx; ++bin_idx) {
293  for (size_t thread_idx = 0; thread_idx < num_threads; ++thread_idx) {
294  const T thread_bin_z_output = per_thread_z_outputs[thread_idx][bin_idx];
295  if (thread_bin_z_output != null_sentinel_ &&
296  thread_bin_z_output > z_[bin_idx]) {
297  z_[bin_idx] = thread_bin_z_output;
298  }
299  }
300  }
301  });
302 }
303 
304 template <typename T, typename Z>
305 void GeoRaster<T, Z>::fill_bins_from_neighbors(const int64_t neighborhood_fill_radius,
306  const bool fill_only_nulls) {
307  auto timer = DEBUG_TIMER(__func__);
308  std::vector<Z> new_z(num_bins_);
309  tbb::parallel_for(tbb::blocked_range<int64_t>(0, num_y_bins_),
310  [&](const tbb::blocked_range<int64_t>& r) {
311  for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
312  for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
313  const int64_t bin_idx =
314  x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
315  const Z z_val = z_[bin_idx];
316  if (!fill_only_nulls || z_val == null_sentinel_) {
317  new_z[bin_idx] = fill_bin_from_avg_neighbors(
318  x_bin, y_bin, neighborhood_fill_radius);
319  } else {
320  new_z[bin_idx] = z_val;
321  }
322  }
323  }
324  });
325  z_.swap(new_z);
326 }
327 
328 template <typename T, typename Z>
330  const int64_t x_bin,
331  const int64_t y_bin,
332  const int64_t num_bins_radius,
333  std::vector<Z>& neighboring_bins) const {
334  const size_t num_bins_per_dim = num_bins_radius * 2 + 1;
335  CHECK_EQ(neighboring_bins.size(), num_bins_per_dim * num_bins_per_dim);
336  const int64_t end_y_bin_idx = y_bin + num_bins_radius;
337  const int64_t end_x_bin_idx = x_bin + num_bins_radius;
338  size_t output_bin = 0;
339  for (int64_t y_bin_idx = y_bin - num_bins_radius; y_bin_idx <= end_y_bin_idx;
340  ++y_bin_idx) {
341  for (int64_t x_bin_idx = x_bin - num_bins_radius; x_bin_idx <= end_x_bin_idx;
342  ++x_bin_idx) {
343  if (x_bin_idx < 0 || x_bin_idx >= num_x_bins_ || y_bin_idx < 0 ||
344  y_bin_idx >= num_y_bins_) {
345  return false;
346  }
347  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin_idx, y_bin_idx, num_x_bins_);
348  neighboring_bins[output_bin++] = z_[bin_idx];
349  if (z_[bin_idx] == null_sentinel_) {
350  return false;
351  }
352  }
353  }
354  return true; // not_null
355 }
356 
357 template <typename T, typename Z>
359  const std::vector<Z>& neighboring_cells,
360  const bool compute_slope_in_degrees) const {
361  const Z dz_dx =
362  ((neighboring_cells[8] + 2 * neighboring_cells[5] + neighboring_cells[2]) -
363  (neighboring_cells[6] + 2 * neighboring_cells[3] + neighboring_cells[0])) /
364  (8 * bin_dim_meters_);
365  const Z dz_dy =
366  ((neighboring_cells[6] + 2 * neighboring_cells[7] + neighboring_cells[8]) -
367  (neighboring_cells[0] + 2 * neighboring_cells[1] + neighboring_cells[2])) /
368  (8 * bin_dim_meters_);
369  const Z slope = sqrt(dz_dx * dz_dx + dz_dy * dz_dy);
370  std::pair<Z, Z> slope_and_aspect;
371  slope_and_aspect.first =
372  compute_slope_in_degrees ? atan(slope) * math_consts::radians_to_degrees : slope;
373  if (slope < 0.0001) {
374  slope_and_aspect.second = null_sentinel_;
375  } else {
376  const Z aspect_degrees =
377  math_consts::radians_to_degrees * atan2(dz_dx, dz_dy); // -180.0 to 180.0
378  slope_and_aspect.second = aspect_degrees + 180.0;
379  // aspect_degrees < 0.0 ? 180.0 + aspect_degrees : aspect_degrees;
380  }
381  return slope_and_aspect;
382 }
383 
384 template <typename T, typename Z>
386  Column<Z>& slope,
387  Column<Z>& aspect,
388  const bool compute_slope_in_degrees) const {
389  auto timer = DEBUG_TIMER(__func__);
390  CHECK_EQ(slope.size(), num_bins_);
392  tbb::blocked_range<int64_t>(0, num_y_bins_),
393  [&](const tbb::blocked_range<int64_t>& r) {
394  std::vector<Z> neighboring_z_vals(9); // 3X3 calc
395  for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
396  for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
397  const bool not_null =
398  get_nxn_neighbors_if_not_null(x_bin, y_bin, 1, neighboring_z_vals);
399  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
400  if (!not_null) {
401  slope.setNull(bin_idx);
402  aspect.setNull(bin_idx);
403  } else {
404  const auto slope_and_aspect = calculate_slope_and_aspect_of_cell(
405  neighboring_z_vals, compute_slope_in_degrees);
406  slope[bin_idx] = slope_and_aspect.first;
407  if (slope_and_aspect.second == null_sentinel_) {
408  aspect.setNull(bin_idx);
409  } else {
410  aspect[bin_idx] = slope_and_aspect.second;
411  }
412  }
413  }
414  }
415  });
416 }
417 
418 template <typename T, typename Z>
421  Column<T>& output_x,
422  Column<T>& output_y,
423  Column<Z>& output_z,
424  const int64_t neighborhood_null_fill_radius) const {
425  auto timer = DEBUG_TIMER(__func__);
426  mgr.set_output_row_size(num_bins_);
428  tbb::blocked_range<int64_t>(0, num_y_bins_),
429  [&](const tbb::blocked_range<int64_t>& r) {
430  for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
431  for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
432  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
433  output_x[bin_idx] = x_min_ + (x_bin + 0.5) * x_scale_bin_to_input_;
434  output_y[bin_idx] = y_min_ + (y_bin + 0.5) * y_scale_bin_to_input_;
435  const Z z_val = z_[bin_idx];
436  if (z_val == null_sentinel_) {
437  output_z.setNull(bin_idx);
438  if (neighborhood_null_fill_radius) {
439  const Z avg_neighbor_value = fill_bin_from_avg_neighbors(
440  x_bin, y_bin, neighborhood_null_fill_radius);
441  if (avg_neighbor_value != null_sentinel_) {
442  output_z[bin_idx] = avg_neighbor_value;
443  }
444  }
445  } else {
446  output_z[bin_idx] = z_[bin_idx];
447  }
448  }
449  }
450  });
451  return num_bins_;
452 }
453 
454 template <typename T, typename Z>
456  Column<T>& output_x,
457  Column<T>& output_y,
458  Column<Z>& output_z) const {
459  auto timer = DEBUG_TIMER(__func__);
460  mgr.set_output_row_size(num_bins_);
462  tbb::blocked_range<int64_t>(0, num_y_bins_),
463  [&](const tbb::blocked_range<int64_t>& r) {
464  for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
465  for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
466  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
467  output_x[bin_idx] = x_min_ + (x_bin + 0.5) * x_scale_bin_to_input_;
468  output_y[bin_idx] = y_min_ + (y_bin + 0.5) * y_scale_bin_to_input_;
469  const Z z_val = z_[bin_idx];
470  if (z_val == null_sentinel_) {
471  output_z.setNull(bin_idx);
472  } else {
473  output_z[bin_idx] = z_[bin_idx];
474  }
475  }
476  }
477  });
478  return num_bins_;
479 }
480 
481 #endif // __CUDACC__
void set_output_row_size(int64_t num_rows)
Definition: OmniSciTypes.h:334
#define CHECK_EQ(x, y)
Definition: Logger.h:219
constexpr double radians_to_degrees
Definition: math_consts.h:22
Z offset_source_z_from_raster_z(const int64_t source_x_bin, const int64_t source_y_bin, const Z source_z_offset) const
bool get_nxn_neighbors_if_not_null(const int64_t x_bin, const int64_t y_bin, const int64_t num_bins_radius, std::vector< Z > &neighboring_bins) const
std::pair< Z, Z > calculate_slope_and_aspect_of_cell(const std::vector< Z > &neighboring_cells, const bool compute_slope_in_degrees) const
int64_t outputDenseColumns(TableFunctionManager &mgr, Column< T > &output_x, Column< T > &output_y, Column< Z > &output_z) const
DEVICE int64_t size() const
Definition: OmniSciTypes.h:248
TEMPLATE_NOINLINE std::pair< T, T > get_column_min_max(const Column< T > &col)
void compute(const Column< T2 > &input_x, const Column< T2 > &input_y, const Column< Z2 > &input_z)
Z fill_bin_from_avg_neighbors(const int64_t x_centroid_bin, const int64_t y_centroid_bin, const int64_t bins_radius) const
const size_t max_inputs_per_thread
void computeParallel(const Column< T2 > &input_x, const Column< T2 > &input_y, const Column< Z2 > &input_z, const size_t max_inputs_per_thread)
GeoRaster(const Column< T2 > &input_x, const Column< T2 > &input_y, const Column< Z2 > &input_z, const double bin_dim_meters, const bool geographic_coords, const bool align_bins_to_zero_based_grid)
int count
const size_t max_temp_output_entries
void fill_bins_from_neighbors(const int64_t neighborhood_fill_radius, const bool fill_only_nulls)
const bool geographic_coords_
void calculate_slope_and_aspect(Column< Z > &slope, Column< Z > &aspect, const bool compute_slope_in_degrees) const
DEVICE bool isNull(int64_t index) const
Definition: OmniSciTypes.h:250
DEVICE void setNull(int64_t index)
Definition: OmniSciTypes.h:251
void parallel_for(const blocked_range< Int > &range, const Body &body, const Partitioner &p=Partitioner())
EXTENSION_NOINLINE double distance_in_meters(const double fromlon, const double fromlat, const double tolon, const double tolat)
Computes the distance, in meters, between two WGS-84 positions.
#define DEBUG_TIMER(name)
Definition: Logger.h:358
char * t
int64_t x_y_bin_to_bin_index(const int64_t x_bin, const int64_t y_bin, const int64_t num_x_bins)