22 #include <tbb/parallel_for.h>
23 #include <tbb/task_arena.h>
33 template <
typename T,
typename Z>
34 template <
typename T2,
typename Z2>
39 const double bin_dim_meters,
40 const bool geographic_coords,
41 const bool align_bins_to_zero_based_grid)
42 : raster_agg_type_(raster_agg_type)
43 , bin_dim_meters_(bin_dim_meters)
44 , geographic_coords_(geographic_coords)
47 const int64_t input_size{input_z.
size()};
48 if (input_size <= 0) {
78 template <
typename T,
typename Z>
79 template <
typename T2,
typename Z2>
84 const double bin_dim_meters,
85 const bool geographic_coords,
86 const bool align_bins_to_zero_based_grid,
91 : raster_agg_type_(raster_agg_type)
92 , bin_dim_meters_(bin_dim_meters)
93 , geographic_coords_(geographic_coords)
113 template <
typename T,
typename Z>
115 const int64_t source_y_bin,
116 const Z source_z_offset)
const {
117 if (is_bin_out_of_bounds(source_x_bin, source_y_bin)) {
118 return null_sentinel_;
121 if (terrain_z == null_sentinel_) {
124 return terrain_z + source_z_offset;
127 template <
typename T,
typename Z>
129 x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
130 x_max_ = std::floor(x_max_ / bin_dim_meters_) * bin_dim_meters_ +
132 y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
133 y_max_ = std::floor(y_max_ / bin_dim_meters_) * bin_dim_meters_ +
137 template <
typename T,
typename Z>
139 x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
140 x_max_ = std::ceil(x_max_ / bin_dim_meters_) * bin_dim_meters_;
141 y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
142 y_max_ = std::ceil(y_max_ / bin_dim_meters_) * bin_dim_meters_;
145 template <
typename T,
typename Z>
147 x_range_ = x_max_ - x_min_;
148 y_range_ = y_max_ - y_min_;
149 if (geographic_coords_) {
150 const T x_centroid = (x_min_ + x_max_) * 0.5;
151 const T y_centroid = (y_min_ + y_max_) * 0.5;
152 x_meters_per_degree_ =
155 y_meters_per_degree_ =
158 num_x_bins_ = x_range_ * x_meters_per_degree_ / bin_dim_meters_;
159 num_y_bins_ = y_range_ * y_meters_per_degree_ / bin_dim_meters_;
161 x_scale_input_to_bin_ = x_meters_per_degree_ / bin_dim_meters_;
162 y_scale_input_to_bin_ = y_meters_per_degree_ / bin_dim_meters_;
163 x_scale_bin_to_input_ = bin_dim_meters_ / x_meters_per_degree_;
164 y_scale_bin_to_input_ = bin_dim_meters_ / y_meters_per_degree_;
167 num_x_bins_ = x_range_ / bin_dim_meters_;
168 num_y_bins_ = y_range_ / bin_dim_meters_;
170 x_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
171 y_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
172 x_scale_bin_to_input_ = bin_dim_meters_;
173 y_scale_bin_to_input_ = bin_dim_meters_;
175 num_bins_ = num_x_bins_ * num_y_bins_;
178 template <
typename T,
typename Z>
179 template <RasterAggType AggType,
typename T2,
typename Z2>
183 std::vector<Z>& output_z) {
185 const int64_t input_size{input_z.
size()};
187 ? std::numeric_limits<Z>::max()
188 : std::numeric_limits<Z>::lowest();
189 output_z.resize(num_bins_, agg_sentinel);
191 for (int64_t sparse_idx = 0; sparse_idx != input_size; ++sparse_idx) {
192 const int64_t x_bin = get_x_bin(input_x[sparse_idx]);
193 const int64_t y_bin = get_y_bin(input_y[sparse_idx]);
194 if (x_bin < 0 || x_bin >= num_x_bins_ || y_bin < 0 || y_bin >= num_y_bins_) {
198 compute_agg(input_z[sparse_idx], output_z[bin_idx], inline_null_value<Z2>());
200 for (int64_t bin_idx = 0; bin_idx != num_bins_; ++bin_idx) {
201 if (output_z[bin_idx] == agg_sentinel) {
202 output_z[bin_idx] = null_sentinel_;
207 template <
typename T,
typename Z>
208 template <RasterAggType AggType>
210 const std::vector<std::vector<Z>>& z_inputs,
211 std::vector<Z>& output_z,
212 const Z agg_sentinel) {
213 const size_t num_inputs = z_inputs.size();
214 output_z.resize(num_bins_, agg_sentinel);
218 tbb::blocked_range<size_t>(0, num_bins_), [&](
const tbb::blocked_range<size_t>& r) {
219 const size_t start_idx = r.begin();
220 const size_t end_idx = r.end();
221 for (
size_t bin_idx = start_idx; bin_idx != end_idx; ++bin_idx) {
222 for (
size_t input_idx = 0; input_idx < num_inputs; ++input_idx) {
223 reduction_agg(z_inputs[input_idx][bin_idx], output_z[bin_idx], agg_sentinel);
229 [&](
const tbb::blocked_range<size_t>& r) {
230 const size_t start_idx = r.begin();
231 const size_t end_idx = r.end();
232 for (
size_t bin_idx = start_idx; bin_idx != end_idx; ++bin_idx) {
233 if (output_z[bin_idx] == agg_sentinel) {
234 output_z[bin_idx] = null_sentinel_;
240 template <
typename T,
typename Z>
241 template <RasterAggType AggType,
typename T2,
typename Z2>
245 std::vector<Z>& output_z,
247 const size_t input_size = input_z.
size();
248 const size_t max_thread_count = std::thread::hardware_concurrency();
249 const size_t num_threads_by_input_elements =
250 std::min(max_thread_count,
251 ((input_size + max_inputs_per_thread - 1) / max_inputs_per_thread));
252 const size_t num_threads_by_output_size =
254 const size_t num_threads =
255 std::min(num_threads_by_input_elements, num_threads_by_output_size);
256 if (num_threads <= 1) {
257 computeSerialImpl<AggType, T2, Z2>(input_x, input_y, input_z, output_z);
262 std::vector<std::vector<Z>> per_thread_z_outputs(num_threads);
265 ? std::numeric_limits<Z>::max()
266 : std::numeric_limits<Z>::lowest();
269 [&](
const tbb::blocked_range<size_t>& r) {
270 for (
size_t t = r.begin(); t != r.end(); ++t) {
271 per_thread_z_outputs[t].resize(num_bins_, agg_sentinel);
276 tbb::task_arena limited_arena(num_threads);
277 limited_arena.execute([&] {
279 tbb::blocked_range<int64_t>(0, input_size),
280 [&](
const tbb::blocked_range<int64_t>& r) {
281 const int64_t start_idx = r.begin();
282 const int64_t end_idx = r.end();
283 size_t thread_idx = tbb::this_task_arena::current_thread_index();
284 std::vector<Z>& this_thread_z_output = per_thread_z_outputs[thread_idx];
286 for (int64_t sparse_idx = start_idx; sparse_idx != end_idx; ++sparse_idx) {
287 const int64_t x_bin = get_x_bin(input_x[sparse_idx]);
288 const int64_t y_bin = get_y_bin(input_y[sparse_idx]);
289 if (x_bin < 0 || x_bin >= num_x_bins_ || y_bin < 0 || y_bin >= num_y_bins_) {
293 compute_agg(input_z[sparse_idx],
294 this_thread_z_output[bin_idx],
295 inline_null_value<Z2>());
303 computeParallelReductionAggImpl<RasterAggType::SUM>(
304 per_thread_z_outputs, output_z, agg_sentinel);
306 computeParallelReductionAggImpl<AggType>(
307 per_thread_z_outputs, output_z, agg_sentinel);
311 template <
typename T,
typename Z>
312 template <
typename T2,
typename Z2>
317 switch (raster_agg_type_) {
319 computeParallelImpl<RasterAggType::COUNT, T2, Z2>(
324 computeParallelImpl<RasterAggType::MIN, T2, Z2>(
329 computeParallelImpl<RasterAggType::MAX, T2, Z2>(
334 computeParallelImpl<RasterAggType::SUM, T2, Z2>(
339 computeParallelImpl<RasterAggType::SUM, T2, Z2>(
341 computeParallelImpl<RasterAggType::COUNT, T2, Z2>(
344 [&](
const tbb::blocked_range<size_t>& r) {
345 const size_t start_idx = r.begin();
346 const size_t end_idx = r.end();
347 for (
size_t bin_idx = start_idx; bin_idx != end_idx;
351 if (counts_[bin_idx] > 1) {
352 z_[bin_idx] /= counts_[bin_idx];
364 template <
typename T,
typename Z>
366 const int64_t x_centroid_bin,
367 const int64_t y_centroid_bin,
368 const int64_t bins_radius)
const {
371 for (int64_t y_bin = y_centroid_bin - bins_radius;
372 y_bin <= y_centroid_bin + bins_radius;
374 for (int64_t x_bin = x_centroid_bin - bins_radius;
375 x_bin <= x_centroid_bin + bins_radius;
377 if (x_bin >= 0 && x_bin < num_x_bins_ && y_bin >= 0 && y_bin < num_y_bins_) {
379 const Z bin_val = z_[bin_idx];
380 if (bin_val != null_sentinel_) {
387 return (count == 0) ? null_sentinel_ : val / count;
391 const int64_t kernel_size = fill_radius * 2 + 1;
392 std::vector<double> gaussian_kernel(kernel_size);
393 const double expr = 1.0 / (sigma * sqrt(2.0 *
M_PI));
394 for (int64_t kernel_idx = -fill_radius; kernel_idx <= fill_radius; ++kernel_idx) {
395 gaussian_kernel[kernel_idx + fill_radius] =
396 expr * exp((kernel_idx * kernel_idx) / (-2.0 * (sigma * sigma)));
398 return gaussian_kernel;
401 template <
typename T,
typename Z>
403 const int64_t neighborhood_fill_radius,
404 const bool fill_only_nulls) {
405 const double sigma = neighborhood_fill_radius * 2.0 / 6.0;
406 const auto gaussian_kernel =
408 std::vector<Z> only_nulls_source_z;
409 if (fill_only_nulls) {
411 only_nulls_source_z.resize(z_.size());
413 [&](
const tbb::blocked_range<int64_t>& r) {
414 const int64_t end_bin_idx = r.end();
415 for (int64_t bin_idx = r.begin(); bin_idx < end_bin_idx;
417 only_nulls_source_z[bin_idx] = z_[bin_idx];
421 auto& source_z = fill_only_nulls ? only_nulls_source_z : z_;
423 std::vector<Z> new_z(num_bins_);
425 tbb::blocked_range<int64_t>(0, num_y_bins_),
426 [&](
const tbb::blocked_range<int64_t>& r) {
427 const int64_t end_y_bin = r.end();
428 for (int64_t y_start_bin = r.begin(); y_start_bin != end_y_bin; ++y_start_bin) {
429 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
432 for (int64_t y_offset = -neighborhood_fill_radius;
433 y_offset <= neighborhood_fill_radius;
435 const auto y_bin = y_start_bin + y_offset;
436 if (y_bin >= 0 && y_bin < num_y_bins_) {
438 const Z bin_val = source_z[bin_idx];
439 if (bin_val != null_sentinel_) {
440 const auto gaussian_weight =
441 gaussian_kernel[y_offset + neighborhood_fill_radius];
442 val += bin_val * gaussian_weight;
443 sum_weights += gaussian_weight;
448 new_z[bin_idx] = sum_weights > 0.0 ? (val / sum_weights) : null_sentinel_;
452 source_z.swap(new_z);
454 tbb::blocked_range<int64_t>(0, num_x_bins_),
455 [&](
const tbb::blocked_range<int64_t>& r) {
456 const int64_t end_x_bin = r.end();
457 for (int64_t x_start_bin = r.begin(); x_start_bin != end_x_bin; ++x_start_bin) {
458 for (int64_t y_bin = 0; y_bin < num_y_bins_; ++y_bin) {
461 for (int64_t x_offset = -neighborhood_fill_radius;
462 x_offset <= neighborhood_fill_radius;
464 const auto x_bin = x_start_bin + x_offset;
465 if (x_bin >= 0 && x_bin < num_x_bins_) {
467 const Z bin_val = source_z[bin_idx];
468 if (bin_val != null_sentinel_) {
469 const auto gaussian_weight =
470 gaussian_kernel[x_offset + neighborhood_fill_radius];
471 val += bin_val * gaussian_weight;
472 sum_weights += gaussian_weight;
477 new_z[bin_idx] = sum_weights > 0.0 ? (val / sum_weights) : null_sentinel_;
481 if (fill_only_nulls) {
485 [&](
const tbb::blocked_range<int64_t>& r) {
486 const int64_t end_bin_idx = r.end();
487 for (int64_t bin_idx = r.begin(); bin_idx < end_bin_idx;
489 if (z_[bin_idx] == null_sentinel_) {
490 z_[bin_idx] = new_z[bin_idx];
492 only_nulls_source_z[bin_idx] = z_[bin_idx];
496 source_z.swap(new_z);
500 template <
typename T,
typename Z>
501 template <RasterAggType AggType>
503 const int64_t x_centroid_bin,
504 const int64_t y_centroid_bin,
505 const int64_t bins_radius,
508 : std::numeric_limits<Z>::lowest();
510 for (int64_t y_bin = y_centroid_bin - bins_radius;
511 y_bin <= y_centroid_bin + bins_radius;
513 for (int64_t x_bin = x_centroid_bin - bins_radius;
514 x_bin <= x_centroid_bin + bins_radius;
516 if (x_bin >= 0 && x_bin < num_x_bins_ && y_bin >= 0 && y_bin < num_y_bins_) {
518 compute_agg(z_[bin_idx], val, null_sentinel_);
522 return val == agg_sentinel ? null_sentinel_ : val;
525 template <
typename T,
typename Z>
526 template <RasterAggType AggType>
528 const int64_t neighborhood_fill_radius,
529 const bool fill_only_nulls) {
533 std::vector<Z> new_z(num_bins_);
536 [&](
const tbb::blocked_range<int64_t>& r) {
537 const int64_t end_y_bin = r.end();
538 for (int64_t y_bin = r.begin(); y_bin != end_y_bin; ++y_bin) {
539 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
540 const int64_t bin_idx =
542 const Z z_val = z_[bin_idx];
543 if (!fill_only_nulls || z_val == null_sentinel_) {
545 new_z[bin_idx] = fill_bin_from_avg_box_neighborhood(
546 x_bin, y_bin, neighborhood_fill_radius);
548 new_z[bin_idx] = fill_bin_from_box_neighborhood(
549 x_bin, y_bin, neighborhood_fill_radius, compute_agg);
552 new_z[bin_idx] = z_val;
560 template <
typename T,
typename Z>
562 const bool fill_only_nulls,
571 switch (raster_fill_agg_type) {
573 fill_bins_from_box_neighborhood<RasterAggType::COUNT>(neighborhood_fill_radius,
578 fill_bins_from_box_neighborhood<RasterAggType::MIN>(neighborhood_fill_radius,
583 fill_bins_from_box_neighborhood<RasterAggType::MAX>(neighborhood_fill_radius,
588 fill_bins_from_box_neighborhood<RasterAggType::SUM>(neighborhood_fill_radius,
593 fill_bins_from_box_neighborhood<RasterAggType::BOX_AVG>(neighborhood_fill_radius,
598 fill_bins_from_gaussian_neighborhood(neighborhood_fill_radius, fill_only_nulls);
607 template <
typename T,
typename Z>
611 const int64_t num_bins_radius,
612 std::vector<Z>& neighboring_bins)
const {
613 const size_t num_bins_per_dim = num_bins_radius * 2 + 1;
614 CHECK_EQ(neighboring_bins.size(), num_bins_per_dim * num_bins_per_dim);
615 const int64_t end_y_bin_idx = y_bin + num_bins_radius;
616 const int64_t end_x_bin_idx = x_bin + num_bins_radius;
617 size_t output_bin = 0;
618 for (int64_t y_bin_idx = y_bin - num_bins_radius; y_bin_idx <= end_y_bin_idx;
620 for (int64_t x_bin_idx = x_bin - num_bins_radius; x_bin_idx <= end_x_bin_idx;
622 if (x_bin_idx < 0 || x_bin_idx >= num_x_bins_ || y_bin_idx < 0 ||
623 y_bin_idx >= num_y_bins_) {
627 neighboring_bins[output_bin++] = z_[bin_idx];
628 if (z_[bin_idx] == null_sentinel_) {
636 template <
typename T,
typename Z>
638 const std::vector<Z>& neighboring_cells,
639 const bool compute_slope_in_degrees)
const {
641 ((neighboring_cells[8] + 2 * neighboring_cells[5] + neighboring_cells[2]) -
642 (neighboring_cells[6] + 2 * neighboring_cells[3] + neighboring_cells[0])) /
643 (8 * bin_dim_meters_);
645 ((neighboring_cells[6] + 2 * neighboring_cells[7] + neighboring_cells[8]) -
646 (neighboring_cells[0] + 2 * neighboring_cells[1] + neighboring_cells[2])) /
647 (8 * bin_dim_meters_);
648 const Z slope = sqrt(dz_dx * dz_dx + dz_dy * dz_dy);
649 std::pair<Z, Z> slope_and_aspect;
650 slope_and_aspect.first =
652 if (slope < 0.0001) {
653 slope_and_aspect.second = null_sentinel_;
655 const Z aspect_degrees =
657 slope_and_aspect.second = aspect_degrees + 180.0;
660 return slope_and_aspect;
663 template <
typename T,
typename Z>
667 const bool compute_slope_in_degrees)
const {
671 tbb::blocked_range<int64_t>(0, num_y_bins_),
672 [&](
const tbb::blocked_range<int64_t>& r) {
673 std::vector<Z> neighboring_z_vals(9);
674 for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
675 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
676 const bool not_null =
677 get_nxn_neighbors_if_not_null(x_bin, y_bin, 1, neighboring_z_vals);
683 const auto slope_and_aspect = calculate_slope_and_aspect_of_cell(
684 neighboring_z_vals, compute_slope_in_degrees);
685 slope[bin_idx] = slope_and_aspect.first;
686 if (slope_and_aspect.second == null_sentinel_) {
689 aspect[bin_idx] = slope_and_aspect.second;
697 template <
typename T,
typename Z>
703 const int64_t neighborhood_null_fill_radius)
const {
707 tbb::blocked_range<int64_t>(0, num_y_bins_),
708 [&](
const tbb::blocked_range<int64_t>& r) {
709 for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
710 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
712 output_x[bin_idx] = x_min_ + (x_bin + 0.5) * x_scale_bin_to_input_;
713 output_y[bin_idx] = y_min_ + (y_bin + 0.5) * y_scale_bin_to_input_;
714 const Z z_val = z_[bin_idx];
715 if (z_val == null_sentinel_) {
717 if (neighborhood_null_fill_radius) {
718 const Z avg_neighbor_value = fill_bin_from_avg_box_neighborhood(
719 x_bin, y_bin, neighborhood_null_fill_radius);
720 if (avg_neighbor_value != null_sentinel_) {
721 output_z[bin_idx] = avg_neighbor_value;
725 output_z[bin_idx] = z_[bin_idx];
733 template <
typename T,
typename Z>
741 tbb::blocked_range<int64_t>(0, num_y_bins_),
742 [&](
const tbb::blocked_range<int64_t>& r) {
743 for (int64_t y_bin = r.begin(); y_bin != r.end(); ++y_bin) {
744 for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
746 output_x[bin_idx] = x_min_ + (x_bin + 0.5) * x_scale_bin_to_input_;
747 output_y[bin_idx] = y_min_ + (y_bin + 0.5) * y_scale_bin_to_input_;
748 const Z z_val = z_[bin_idx];
749 if (z_val == null_sentinel_) {
752 output_z[bin_idx] = z_[bin_idx];
760 template <
typename T,
typename Z>
764 mgr.
set_metadata(
"geo_raster_x_min", static_cast<double>(x_min_));
765 mgr.
set_metadata(
"geo_raster_y_min", static_cast<double>(y_min_));
767 static_cast<double>(x_scale_input_to_bin_));
769 static_cast<double>(y_scale_input_to_bin_));
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
void setMetadata(TableFunctionManager &mgr) 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)
GeoRaster(const Column< T2 > &input_x, const Column< T2 > &input_y, const Column< Z2 > &input_z, const RasterAggType raster_agg_type, const double bin_dim_meters, const bool geographic_coords, const bool align_bins_to_zero_based_grid)
int64_t outputDenseColumns(TableFunctionManager &mgr, Column< T > &output_x, Column< T > &output_y, Column< Z > &output_z) const
DEVICE int64_t size() const
void fill_bins_from_neighbors(const int64_t neighborhood_fill_radius, const bool fill_only_nulls, const RasterAggType raster_agg_type=RasterAggType::GAUSS_AVG)
void computeParallelReductionAggImpl(const std::vector< std::vector< Z >> &z_inputs, std::vector< Z > &output_z, const Z agg_sentinel)
const size_t max_inputs_per_thread
DEVICE TextEncodingDict inline_null_value()
const size_t max_temp_output_entries
void align_bins_max_exclusive()
std::vector< double > generate_1d_gaussian_kernel(const int64_t fill_radius, double sigma)
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 fill_bins_from_gaussian_neighborhood(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
void computeParallelImpl(const Column< T2 > &input_x, const Column< T2 > &input_y, const Column< Z2 > &input_z, std::vector< Z > &output_z, const size_t max_inputs_per_thread)
Z fill_bin_from_box_neighborhood(const int64_t x_centroid_bin, const int64_t y_centroid_bin, const int64_t bins_radius, const ComputeAgg< AggType > &compute_agg) const
void fill_bins_from_box_neighborhood(const int64_t neighborhood_fill_radius, const bool fill_only_nulls)
DEVICE void setNull(int64_t index)
void compute(const Column< T2 > &input_x, const Column< T2 > &input_y, const Column< Z2 > &input_z, const size_t max_inputs_per_thread)
Z fill_bin_from_avg_box_neighborhood(const int64_t x_centroid_bin, const int64_t y_centroid_bin, const int64_t bins_radius) const
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.
void set_metadata(const std::string &key, const T &value)
#define DEBUG_TIMER(name)
void align_bins_max_inclusive()
void calculate_bins_and_scales()
void computeSerialImpl(const Column< T2 > &input_x, const Column< T2 > &input_y, const Column< Z2 > &input_z, std::vector< Z > &output_z)