22 #include <tbb/parallel_for.h>
23 #include <tbb/task_arena.h>
33 template <
typename T,
typename Z>
34 template <
typename T2,
typename Z2>
38 const double bin_dim_meters,
39 const bool geographic_coords,
40 const bool align_bins_to_zero_based_grid)
41 : bin_dim_meters_(bin_dim_meters)
42 , geographic_coords_(geographic_coords)
43 , null_sentinel_(std::numeric_limits<Z>::lowest()) {
45 const int64_t input_size{input_z.
size()};
46 if (input_size <= 0) {
76 template <
typename T,
typename Z>
77 template <
typename T2,
typename Z2>
81 const double bin_dim_meters,
82 const bool geographic_coords,
83 const bool align_bins_to_zero_based_grid,
88 : bin_dim_meters_(bin_dim_meters)
89 , geographic_coords_(geographic_coords)
90 , null_sentinel_(std::numeric_limits<Z>::lowest())
109 template <
typename T,
typename Z>
111 const int64_t source_y_bin,
112 const Z source_z_offset)
const {
113 if (is_bin_out_of_bounds(source_x_bin, source_y_bin)) {
114 return null_sentinel_;
117 if (terrain_z == null_sentinel_) {
120 return terrain_z + source_z_offset;
123 template <
typename T,
typename Z>
125 const int64_t y_centroid_bin,
126 const int64_t bins_radius)
const {
129 for (int64_t y_bin = y_centroid_bin - bins_radius;
130 y_bin <= y_centroid_bin + bins_radius;
132 for (int64_t x_bin = x_centroid_bin - bins_radius;
133 x_bin <= x_centroid_bin + bins_radius;
135 if (x_bin >= 0 && x_bin < num_x_bins_ && y_bin >= 0 && y_bin < num_y_bins_) {
137 const Z bin_val = z_[bin_idx];
138 if (bin_val != null_sentinel_) {
145 return (count == 0) ? null_sentinel_ : val / count;
148 template <
typename T,
typename Z>
150 x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
151 x_max_ = std::floor(x_max_ / bin_dim_meters_) * bin_dim_meters_ +
153 y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
154 y_max_ = std::floor(y_max_ / bin_dim_meters_) * bin_dim_meters_ +
158 template <
typename T,
typename Z>
160 x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
161 x_max_ = std::ceil(x_max_ / bin_dim_meters_) * bin_dim_meters_;
162 y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
163 y_max_ = std::ceil(y_max_ / bin_dim_meters_) * bin_dim_meters_;
166 template <
typename T,
typename Z>
168 x_range_ = x_max_ - x_min_;
169 y_range_ = y_max_ - y_min_;
170 if (geographic_coords_) {
171 const T x_centroid = (x_min_ + x_max_) * 0.5;
172 const T y_centroid = (y_min_ + y_max_) * 0.5;
173 x_meters_per_degree_ =
176 y_meters_per_degree_ =
179 num_x_bins_ = x_range_ * x_meters_per_degree_ / bin_dim_meters_;
180 num_y_bins_ = y_range_ * y_meters_per_degree_ / bin_dim_meters_;
182 x_scale_input_to_bin_ = x_meters_per_degree_ / bin_dim_meters_;
183 y_scale_input_to_bin_ = y_meters_per_degree_ / bin_dim_meters_;
184 x_scale_bin_to_input_ = bin_dim_meters_ / x_meters_per_degree_;
185 y_scale_bin_to_input_ = bin_dim_meters_ / y_meters_per_degree_;
188 num_x_bins_ = x_range_ / bin_dim_meters_;
189 num_y_bins_ = y_range_ / bin_dim_meters_;
191 x_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
192 y_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
193 x_scale_bin_to_input_ = bin_dim_meters_;
194 y_scale_bin_to_input_ = bin_dim_meters_;
196 num_bins_ = num_x_bins_ * num_y_bins_;
199 template <
typename T,
typename Z>
200 template <
typename T2,
typename Z2>
205 const int64_t input_size{input_z.
size()};
206 z_.resize(num_bins_, null_sentinel_);
207 for (int64_t sparse_idx = 0; sparse_idx != input_size; ++sparse_idx) {
208 const int64_t x_bin = get_x_bin(input_x[sparse_idx]);
209 const int64_t y_bin = get_y_bin(input_y[sparse_idx]);
210 if (x_bin < 0 || x_bin >= num_x_bins_ || y_bin < 0 || y_bin >= num_y_bins_) {
216 if (!(input_z.
isNull(sparse_idx)) && input_z[sparse_idx] > z_[bin_idx]) {
217 z_[bin_idx] = input_z[sparse_idx];
222 template <
typename T,
typename Z>
223 template <
typename T2,
typename Z2>
228 const size_t input_size = input_z.
size();
229 const size_t max_thread_count = std::thread::hardware_concurrency();
230 const size_t num_threads_by_input_elements =
231 std::min(max_thread_count,
232 ((input_size + max_inputs_per_thread - 1) / max_inputs_per_thread));
233 const size_t num_threads_by_output_size =
235 const size_t num_threads =
236 std::min(num_threads_by_input_elements, num_threads_by_output_size);
237 if (num_threads <= 1) {
238 compute(input_x, input_y, input_z);
243 std::vector<std::vector<Z>> per_thread_z_outputs(num_threads);
246 [&](
const tbb::blocked_range<size_t>& r) {
247 for (
size_t t = r.begin(); t != r.end(); ++t) {
248 per_thread_z_outputs[t].resize(num_bins_, null_sentinel_);
252 tbb::task_arena limited_arena(num_threads);
253 limited_arena.execute([&] {
255 tbb::blocked_range<int64_t>(0, input_size),
256 [&](
const tbb::blocked_range<int64_t>& r) {
257 const int64_t start_idx = r.begin();
258 const int64_t end_idx = r.end();
259 size_t thread_idx = tbb::this_task_arena::current_thread_index();
260 std::vector<Z>& this_thread_z_output = per_thread_z_outputs[thread_idx];
262 for (int64_t sparse_idx = start_idx; sparse_idx != end_idx; ++sparse_idx) {
263 const int64_t x_bin = get_x_bin(input_x[sparse_idx]);
264 const int64_t y_bin = get_y_bin(input_y[sparse_idx]);
265 if (x_bin < 0 || x_bin >= num_x_bins_ || y_bin < 0 || y_bin >= num_y_bins_) {
271 if (!(input_z.
isNull(sparse_idx)) &&
272 input_z[sparse_idx] > this_thread_z_output[bin_idx]) {
273 this_thread_z_output[bin_idx] = input_z[sparse_idx];
279 z_.resize(num_bins_, null_sentinel_);
282 tbb::blocked_range<size_t>(0, num_bins_), [&](
const tbb::blocked_range<size_t>& r) {
283 const size_t start_idx = r.begin();
284 const size_t end_idx = r.end();
285 for (
size_t bin_idx = start_idx; bin_idx != end_idx; ++bin_idx) {
286 for (
size_t thread_idx = 0; thread_idx < num_threads; ++thread_idx) {
287 const T thread_bin_z_output = per_thread_z_outputs[thread_idx][bin_idx];
288 if (thread_bin_z_output != null_sentinel_ &&
289 thread_bin_z_output > z_[bin_idx]) {
290 z_[bin_idx] = thread_bin_z_output;
297 template <
typename T,
typename Z>
299 const bool fill_only_nulls) {
301 std::vector<Z> new_z(num_bins_);
303 [&](
const tbb::blocked_range<int64_t>& r) {
304 for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
305 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
306 const int64_t bin_idx =
308 const Z z_val = z_[bin_idx];
309 if (!fill_only_nulls || z_val == null_sentinel_) {
310 new_z[bin_idx] = fill_bin_from_avg_neighbors(
311 x_bin, y_bin, neighborhood_fill_radius);
313 new_z[bin_idx] = z_val;
321 template <
typename T,
typename Z>
325 const int64_t num_bins_radius,
326 std::vector<Z>& neighboring_bins)
const {
327 const size_t num_bins_per_dim = num_bins_radius * 2 + 1;
328 CHECK_EQ(neighboring_bins.size(), num_bins_per_dim * num_bins_per_dim);
329 const int64_t end_y_bin_idx = y_bin + num_bins_radius;
330 const int64_t end_x_bin_idx = x_bin + num_bins_radius;
331 size_t output_bin = 0;
332 for (int64_t y_bin_idx = y_bin - num_bins_radius; y_bin_idx <= end_y_bin_idx;
334 for (int64_t x_bin_idx = x_bin - num_bins_radius; x_bin_idx <= end_x_bin_idx;
336 if (x_bin_idx < 0 || x_bin_idx >= num_x_bins_ || y_bin_idx < 0 ||
337 y_bin_idx >= num_y_bins_) {
341 neighboring_bins[output_bin++] = z_[bin_idx];
342 if (z_[bin_idx] == null_sentinel_) {
350 template <
typename T,
typename Z>
352 const std::vector<Z>& neighboring_cells,
353 const bool compute_slope_in_degrees)
const {
355 ((neighboring_cells[8] + 2 * neighboring_cells[5] + neighboring_cells[2]) -
356 (neighboring_cells[6] + 2 * neighboring_cells[3] + neighboring_cells[0])) /
357 (8 * bin_dim_meters_);
359 ((neighboring_cells[6] + 2 * neighboring_cells[7] + neighboring_cells[8]) -
360 (neighboring_cells[0] + 2 * neighboring_cells[1] + neighboring_cells[2])) /
361 (8 * bin_dim_meters_);
362 const Z slope = sqrt(dz_dx * dz_dx + dz_dy * dz_dy);
363 std::pair<Z, Z> slope_and_aspect;
364 slope_and_aspect.first =
366 if (slope < 0.0001) {
367 slope_and_aspect.second = null_sentinel_;
369 const Z aspect_degrees =
371 slope_and_aspect.second = aspect_degrees + 180.0;
374 return slope_and_aspect;
377 template <
typename T,
typename Z>
381 const bool compute_slope_in_degrees)
const {
385 tbb::blocked_range<int64_t>(0, num_y_bins_),
386 [&](
const tbb::blocked_range<int64_t>& r) {
387 std::vector<Z> neighboring_z_vals(9);
388 for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
389 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
390 const bool not_null =
391 get_nxn_neighbors_if_not_null(x_bin, y_bin, 1, neighboring_z_vals);
397 const auto slope_and_aspect = calculate_slope_and_aspect_of_cell(
398 neighboring_z_vals, compute_slope_in_degrees);
399 slope[bin_idx] = slope_and_aspect.first;
400 if (slope_and_aspect.second == null_sentinel_) {
403 aspect[bin_idx] = slope_and_aspect.second;
411 template <
typename T,
typename Z>
417 const int64_t neighborhood_null_fill_radius)
const {
421 tbb::blocked_range<int64_t>(0, num_y_bins_),
422 [&](
const tbb::blocked_range<int64_t>& r) {
423 for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
424 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
426 output_x[bin_idx] = x_min_ + (x_bin + 0.5) * x_scale_bin_to_input_;
427 output_y[bin_idx] = y_min_ + (y_bin + 0.5) * y_scale_bin_to_input_;
428 const Z z_val = z_[bin_idx];
429 if (z_val == null_sentinel_) {
431 if (neighborhood_null_fill_radius) {
432 const Z avg_neighbor_value = fill_bin_from_avg_neighbors(
433 x_bin, y_bin, neighborhood_null_fill_radius);
434 if (avg_neighbor_value != null_sentinel_) {
435 output_z[bin_idx] = avg_neighbor_value;
439 output_z[bin_idx] = z_[bin_idx];
447 template <
typename T,
typename Z>
455 tbb::blocked_range<int64_t>(0, num_y_bins_),
456 [&](
const tbb::blocked_range<int64_t>& r) {
457 for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
458 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
460 output_x[bin_idx] = x_min_ + (x_bin + 0.5) * x_scale_bin_to_input_;
461 output_y[bin_idx] = y_min_ + (y_bin + 0.5) * y_scale_bin_to_input_;
462 const Z z_val = z_[bin_idx];
463 if (z_val == null_sentinel_) {
466 output_z[bin_idx] = z_[bin_idx];
void set_output_row_size(int64_t num_rows)
constexpr double radians_to_degrees
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
NEVER_INLINE HOST std::pair< T, T > get_column_min_max(const Column< T > &col)
int64_t outputDenseColumns(TableFunctionManager &mgr, Column< T > &output_x, Column< T > &output_y, Column< Z > &output_z) const
DEVICE int64_t size() const
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)
const size_t max_temp_output_entries
void fill_bins_from_neighbors(const int64_t neighborhood_fill_radius, const bool fill_only_nulls)
void align_bins_max_exclusive()
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
DEVICE void setNull(int64_t index)
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)
void align_bins_max_inclusive()
int64_t x_y_bin_to_bin_index(const int64_t x_bin, const int64_t y_bin, const int64_t num_x_bins)
void calculate_bins_and_scales()