OmniSciDB  6686921089
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GeoRaster.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 #ifdef HAVE_SYSTEM_TFS
18 #ifndef __CUDACC__
19 
20 #include <cmath>
21 #include <vector>
22 
23 #include "Shared/Utilities.h"
24 
25 // Allow input types to GeoRaster that are different than class types/output Z type
26 // So we can move everything to the type of T and Z (which can each be either float
27 // or double)
28 template <typename T, typename Z>
29 template <typename T2, typename Z2>
30 GeoRaster<T, Z>::GeoRaster(const Column<T2>& input_x,
31  const Column<T2>& input_y,
32  const Column<Z2>& input_z,
33  const double bin_dim_meters,
34  const bool geographic_coords,
35  const bool align_bins_to_zero_based_grid)
36  : bin_dim_meters_(bin_dim_meters)
37  , geographic_coords_(geographic_coords)
38  , null_sentinel_(std::numeric_limits<Z>::lowest()) {
39  const int64_t input_size{input_z.size()};
40  if (input_size <= 0) {
41  num_bins_ = 0;
42  num_x_bins_ = 0;
43  num_y_bins_ = 0;
44  return;
45  }
46  const auto min_max_x = get_column_min_max(input_x);
47  const auto min_max_y = get_column_min_max(input_y);
48  x_min_ = min_max_x.first;
49  x_max_ = min_max_x.second;
50  y_min_ = min_max_y.first;
51  y_max_ = min_max_y.second;
52 
53  if (align_bins_to_zero_based_grid && !geographic_coords_) {
54  // For implicit, data-defined bounds, we treat the max of the x and y ranges as
55  // inclusive (closed interval), since if the max of the data in either x/y dimensions
56  // is at the first value of the next bin, values at that max will be discarded if we
57  // don't include the final bin. For exmaple, if the input data (perhaps already binned
58  // with a group by query) goes from 0.0 to 40.0 in both x and y directions, we should
59  // have the last x/y bins cover the range [40.0, 50.0), not [30.0, 40.0)
60  align_bins_max_inclusive();
61  }
62  // std::cout << "X range: [" << x_min_ << ", " << x_max_ << "]" << std::endl;
63  // std::cout << "Y range: [" << y_min_ << ", " << y_max_ << "]" << std::endl;
64 
65  calculate_bins_and_scales();
66  compute(input_x, input_y, input_z);
67 }
68 
69 // Allow input types to GeoRaster that are different than class types/output Z type
70 // So we can move everything to the type of T and Z (which can each be either float
71 // or double)
72 template <typename T, typename Z>
73 template <typename T2, typename Z2>
74 GeoRaster<T, Z>::GeoRaster(const Column<T2>& input_x,
75  const Column<T2>& input_y,
76  const Column<Z2>& input_z,
77  const double bin_dim_meters,
78  const bool geographic_coords,
79  const bool align_bins_to_zero_based_grid,
80  const T x_min,
81  const T x_max,
82  const T y_min,
83  const T y_max)
84  : bin_dim_meters_(bin_dim_meters)
85  , geographic_coords_(geographic_coords)
86  , null_sentinel_(std::numeric_limits<Z>::lowest())
87  , x_min_(x_min)
88  , x_max_(x_max)
89  , y_min_(y_min)
90  , y_max_(y_max) {
91  if (align_bins_to_zero_based_grid && !geographic_coords_) {
92  // For explicit, user-defined bounds, we treat the max of the x and y ranges as
93  // exclusive (open interval), since if the user specifies the max x/y as the end of
94  // the bin, they do not intend to add the next full bin For example, if a user
95  // specifies a bin_dim_meters of 10.0 and an x and y range from 0 to 40.0, they almost
96  // assuredly intend for there to be 4 bins in each of the x and y dimensions, with the
97  // last bin of range [30.0, 40.0), not 5 with the final bin's range from [40.0, 50.0)
98  align_bins_max_exclusive();
99  }
100  calculate_bins_and_scales();
101  compute(input_x, input_y, input_z);
102 }
103 
104 template <typename T, typename Z>
105 inline Z GeoRaster<T, Z>::offset_source_z_from_raster_z(const int64_t source_x_bin,
106  const int64_t source_y_bin,
107  const Z source_z_offset) const {
108  if (is_bin_out_of_bounds(source_x_bin, source_y_bin)) {
109  return null_sentinel_;
110  }
111  const Z terrain_z = z_[x_y_bin_to_bin_index(source_x_bin, source_y_bin, num_x_bins_)];
112  if (terrain_z == null_sentinel_) {
113  return terrain_z;
114  }
115  return terrain_z + source_z_offset;
116 }
117 
118 template <typename T, typename Z>
119 inline Z GeoRaster<T, Z>::fill_bin_from_avg_neighbors(const int64_t x_centroid_bin,
120  const int64_t y_centroid_bin,
121  const int64_t bins_radius) const {
122  T val = 0.0;
123  int32_t count = 0;
124  for (int64_t y_bin = y_centroid_bin - bins_radius;
125  y_bin <= y_centroid_bin + bins_radius;
126  y_bin++) {
127  for (int64_t x_bin = x_centroid_bin - bins_radius;
128  x_bin <= x_centroid_bin + bins_radius;
129  x_bin++) {
130  if (x_bin >= 0 && x_bin < num_x_bins_ && y_bin >= 0 && y_bin < num_y_bins_) {
131  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
132  const Z bin_val = z_[bin_idx];
133  if (bin_val != null_sentinel_) {
134  count++;
135  val += bin_val;
136  }
137  }
138  }
139  }
140  return (count == 0) ? null_sentinel_ : val / count;
141 }
142 
143 template <typename T, typename Z>
144 void GeoRaster<T, Z>::align_bins_max_inclusive() {
145  x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
146  x_max_ = std::floor(x_max_ / bin_dim_meters_) * bin_dim_meters_ +
147  bin_dim_meters_; // Snap to end of bin
148  y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
149  y_max_ = std::floor(y_max_ / bin_dim_meters_) * bin_dim_meters_ +
150  bin_dim_meters_; // Snap to end of bin
151 }
152 
153 template <typename T, typename Z>
154 void GeoRaster<T, Z>::align_bins_max_exclusive() {
155  x_min_ = std::floor(x_min_ / bin_dim_meters_) * bin_dim_meters_;
156  x_max_ = std::ceil(x_max_ / bin_dim_meters_) * bin_dim_meters_;
157  y_min_ = std::floor(y_min_ / bin_dim_meters_) * bin_dim_meters_;
158  y_max_ = std::ceil(y_max_ / bin_dim_meters_) * bin_dim_meters_;
159 }
160 
161 template <typename T, typename Z>
162 void GeoRaster<T, Z>::calculate_bins_and_scales() {
163  x_range_ = x_max_ - x_min_;
164  y_range_ = y_max_ - y_min_;
165  if (geographic_coords_) {
166  const T x_centroid = (x_min_ + x_max_) * 0.5;
167  const T y_centroid = (y_min_ + y_max_) * 0.5;
168  x_meters_per_degree_ =
169  distance_in_meters(x_min_, y_centroid, x_max_, y_centroid) / x_range_;
170 
171  y_meters_per_degree_ =
172  distance_in_meters(x_centroid, y_min_, x_centroid, y_max_) / y_range_;
173 
174  num_x_bins_ = x_range_ * x_meters_per_degree_ / bin_dim_meters_;
175  num_y_bins_ = y_range_ * y_meters_per_degree_ / bin_dim_meters_;
176 
177  x_scale_input_to_bin_ = x_meters_per_degree_ / bin_dim_meters_;
178  y_scale_input_to_bin_ = y_meters_per_degree_ / bin_dim_meters_;
179  x_scale_bin_to_input_ = bin_dim_meters_ / x_meters_per_degree_;
180  y_scale_bin_to_input_ = bin_dim_meters_ / y_meters_per_degree_;
181 
182  } else {
183  num_x_bins_ = x_range_ / bin_dim_meters_;
184  num_y_bins_ = y_range_ / bin_dim_meters_;
185 
186  x_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
187  y_scale_input_to_bin_ = 1.0 / bin_dim_meters_;
188  x_scale_bin_to_input_ = bin_dim_meters_;
189  y_scale_bin_to_input_ = bin_dim_meters_;
190  }
191  num_bins_ = num_x_bins_ * num_y_bins_;
192 }
193 
194 template <typename T, typename Z>
195 template <typename T2, typename Z2>
196 void GeoRaster<T, Z>::compute(const Column<T2>& input_x,
197  const Column<T2>& input_y,
198  const Column<Z2>& input_z) {
199  const int64_t input_size{input_z.size()};
200  z_.resize(num_bins_, null_sentinel_);
201  for (int64_t sparse_idx = 0; sparse_idx != input_size; ++sparse_idx) {
202  const int64_t x_bin = get_x_bin(input_x[sparse_idx]);
203  const int64_t y_bin = get_y_bin(input_y[sparse_idx]);
204  if (x_bin < 0 || x_bin >= num_x_bins_ || y_bin < 0 || y_bin >= num_y_bins_) {
205  continue;
206  }
207  // Take the max height for this version, but may want to allow different metrics
208  // like average as well
209  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
210  if (!(input_z.isNull(sparse_idx)) && input_z[sparse_idx] > z_[bin_idx]) {
211  z_[bin_idx] = input_z[sparse_idx];
212  }
213  }
214 }
215 
216 template <typename T, typename Z>
217 int64_t GeoRaster<T, Z>::outputDenseColumns(
219  Column<T>& output_x,
220  Column<T>& output_y,
221  Column<Z>& output_z,
222  const int64_t neighborhood_null_fill_radius) const {
223  mgr.set_output_row_size(num_bins_);
224  for (int64_t y_bin = 0; y_bin < num_y_bins_; ++y_bin) {
225  for (int64_t x_bin = 0; x_bin < num_x_bins_; ++x_bin) {
226  const int64_t bin_idx = x_y_bin_to_bin_index(x_bin, y_bin, num_x_bins_);
227  output_x[bin_idx] = x_min_ + (x_bin + 0.5) * x_scale_bin_to_input_;
228  output_y[bin_idx] = y_min_ + (y_bin + 0.5) * y_scale_bin_to_input_;
229  const Z z_val = z_[bin_idx];
230  if (z_val == null_sentinel_) {
231  output_z.setNull(bin_idx);
232  if (neighborhood_null_fill_radius) {
233  const Z avg_neighbor_value =
234  fill_bin_from_avg_neighbors(x_bin, y_bin, neighborhood_null_fill_radius);
235  if (avg_neighbor_value != null_sentinel_) {
236  output_z[bin_idx] = avg_neighbor_value;
237  }
238  }
239  } else {
240  output_z[bin_idx] = z_[bin_idx];
241  }
242  }
243  }
244  return num_bins_;
245 }
246 
247 #endif // __CUDACC__
248 #endif // HAVE_SYSTEM_TFS
void set_output_row_size(int64_t num_rows)
Definition: OmniSciTypes.h:304
DEVICE int64_t size() const
Definition: OmniSciTypes.h:218
int count
DEVICE bool isNull(int64_t index) const
Definition: OmniSciTypes.h:220
DEVICE void setNull(int64_t index)
Definition: OmniSciTypes.h:221
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.