OmniSciDB  2e3a973ef4
GeosRuntime.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2019 OmniSci 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 ENABLE_GEOS
18 
19 #ifndef __CUDACC__
20 
21 #include <cstdarg>
22 #include <mutex>
23 
24 #include "Geospatial/Compression.h"
25 #include "Geospatial/Types.h"
27 #include "Shared/checked_alloc.h"
28 #include "Shared/funcannotations.h"
29 
30 #define GEOS_USE_ONLY_R_API
31 #include <geos_c.h>
32 
33 using namespace Geospatial;
34 
35 using WKB = std::vector<uint8_t>;
36 
37 #define MAX_GEOS_MESSAGE_LEN 200
38 
39 static std::mutex geos_log_info_mutex;
40 static std::mutex geos_log_error_mutex;
41 
42 // called by GEOS on notice
43 static void geos_notice_handler(const char* fmt, ...) {
44  char buffer[MAX_GEOS_MESSAGE_LEN];
45  va_list args;
46  va_start(args, fmt);
47  vsnprintf(buffer, MAX_GEOS_MESSAGE_LEN, fmt, args);
48  va_end(args);
49  {
50  std::lock_guard<std::mutex> guard(geos_log_info_mutex);
51  LOG(INFO) << "GEOS Notice: " << std::string(buffer);
52  }
53 }
54 
55 // called by GEOS on error
56 static void geos_error_handler(const char* fmt, ...) {
57  va_list args;
58  char buffer[MAX_GEOS_MESSAGE_LEN];
59  va_start(args, fmt);
60  vsnprintf(buffer, MAX_GEOS_MESSAGE_LEN, fmt, args);
61  va_end(args);
62  {
63  std::lock_guard<std::mutex> guard(geos_log_error_mutex);
64  LOG(ERROR) << "GEOS Error: " << std::string(buffer);
65  }
66 }
67 
68 GEOSContextHandle_t create_context() {
69 #if GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR < 5
70  GEOSContextHandle_t context = initGEOS_r(geos_notice_handler, geos_error_handler);
71  CHECK(context);
72  return context;
73 #else
74  GEOSContextHandle_t context = GEOS_init_r();
75  CHECK(context);
76  GEOSContext_setNoticeHandler_r(context, geos_notice_handler);
77  GEOSContext_setErrorHandler_r(context, geos_error_handler);
78  return context;
79 #endif
80 }
81 
82 void destroy_context(GEOSContextHandle_t context) {
83  CHECK(context);
84 #if GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR < 5
85  finishGEOS_r(context);
86 #else
87  GEOS_finish_r(context);
88 #endif
89 }
90 
91 bool toWkb(WKB& wkb,
92  int type, // internal geometry type
93  int8_t* coords,
94  int64_t coords_size,
95  int32_t* meta1, // e.g. ring_sizes
96  int64_t meta1_size, // e.g. num_rings
97  int32_t* meta2, // e.g. rings (number of rings in each poly)
98  int64_t meta2_size, // e.g. num_polys
99  int32_t ic // input compression
100 ) {
101  // decompressed double coords
102  auto cv = Geospatial::decompress_coords<double, int32_t>(ic, coords, coords_size);
103  if (static_cast<SQLTypes>(type) == kPOINT) {
104  GeoPoint point(*cv);
105  return point.getWkb(wkb);
106  }
107  if (static_cast<SQLTypes>(type) == kLINESTRING) {
108  GeoLineString linestring(*cv);
109  return linestring.getWkb(wkb);
110  }
111  std::vector<int32_t> meta1v(meta1, meta1 + meta1_size);
112  if (static_cast<SQLTypes>(type) == kPOLYGON) {
113  GeoPolygon poly(*cv, meta1v);
114  return poly.getWkb(wkb);
115  }
116  std::vector<int32_t> meta2v(meta2, meta2 + meta2_size);
117  if (static_cast<SQLTypes>(type) == kMULTIPOLYGON) {
118  // Recognize GEOMETRYCOLLECTION EMPTY encoding
119  // MULTIPOLYGON (((0 0,0.00000012345 0.0,0.0 0.00000012345,0 0)))
120  // Used to pass along EMPTY from ST_Intersection to ST_IsEmpty for example
121  if (meta1_size == 1 && meta2_size == 1) {
122  const std::vector<double> ecv = {0.0, 0.0, 0.00000012345, 0.0, 0.0, 0.00000012345};
123  if (*cv == ecv) {
124  GeoGeometryCollection empty("GEOMETRYCOLLECTION EMPTY");
125  return empty.getWkb(wkb);
126  }
127  }
128  GeoMultiPolygon mpoly(*cv, meta1v, meta2v);
129  return mpoly.getWkb(wkb);
130  }
131  return false;
132 }
133 
134 // Conversion form wkb to internal vector representation.
135 // Each vector components is malloced, caller is reponsible for freeing.
136 bool fromWkb(WKB& wkb,
137  int* result_type,
138  int8_t** result_coords,
139  int64_t* result_coords_size,
140  int32_t** result_meta1,
141  int64_t* result_meta1_size,
142  int32_t** result_meta2,
143  int64_t* result_meta2_size) {
145 
146  // Get the column values
147  std::vector<double> coords{};
148  std::vector<int32_t> ring_sizes{};
149  std::vector<int32_t> poly_rings{};
150  std::vector<double> bounds{};
151 
152  // Forcing MULTIPOLYGON result until we can handle any geo.
153  if (result->isEmpty()) {
154  // Generate a tiny polygon around POINT(0 0), make it a multipolygon
155  // MULTIPOLYGON (((0 0,0.00000012345 0.0,0.0 0.00000012345,0 0)))
156  // to simulate an empty result
157  coords = {0.0, 0.0, 0.00000012345, 0.0, 0.0, 0.00000012345};
158  ring_sizes.push_back(3);
159  poly_rings.push_back(1);
160  } else if (auto result_point = dynamic_cast<GeoPoint*>(result.get())) {
161  result_point->getColumns(coords);
162  // Generate a tiny polygon around the point, make it a multipolygon
163  coords.push_back(coords[0] + 0.0000001);
164  coords.push_back(coords[1]);
165  coords.push_back(coords[0]);
166  coords.push_back(coords[1] + 0.0000001);
167  ring_sizes.push_back(3);
168  poly_rings.push_back(ring_sizes.size());
169  } else if (auto result_poly = dynamic_cast<GeoPolygon*>(result.get())) {
170  result_poly->getColumns(coords, ring_sizes, bounds);
171  // Convert to a 1-polygon multipolygon
172  poly_rings.push_back(ring_sizes.size());
173  } else if (auto result_mpoly = dynamic_cast<GeoMultiPolygon*>(result.get())) {
174  result_mpoly->getColumns(coords, ring_sizes, poly_rings, bounds);
175  } else {
176  return false;
177  }
178 
179  // TODO: consider using a single buffer to hold all components,
180  // instead of allocating and registering each component buffer separately
181 
182  *result_type = static_cast<int>(kMULTIPOLYGON);
183 
184  *result_coords = nullptr;
185  int64_t size = coords.size() * sizeof(double);
186  if (size > 0) {
187  auto buf = checked_malloc(size);
188  std::memcpy(buf, coords.data(), size);
189  *result_coords = reinterpret_cast<int8_t*>(buf);
190  }
191  *result_coords_size = size;
192 
193  *result_meta1 = nullptr;
194  size = ring_sizes.size() * sizeof(int32_t);
195  if (size > 0) {
196  auto buf = checked_malloc(size);
197  std::memcpy(buf, ring_sizes.data(), size);
198  *result_meta1 = reinterpret_cast<int32_t*>(buf);
199  }
200  *result_meta1_size = ring_sizes.size();
201 
202  *result_meta2 = nullptr;
203  size = poly_rings.size() * sizeof(int32_t);
204  if (size > 0) {
205  auto buf = checked_malloc(size);
206  std::memcpy(buf, poly_rings.data(), size);
207  *result_meta2 = reinterpret_cast<int32_t*>(buf);
208  }
209  *result_meta2_size = poly_rings.size();
210 
211  return true;
212 }
213 
214 GEOSGeometry* postprocess(GEOSContextHandle_t context, GEOSGeometry* g) {
215  if (g && GEOSisEmpty_r(context, g) == 0) {
216  auto type = GEOSGeomTypeId_r(context, g);
217  if (type != -1) {
218  if (type != GEOS_POINT && type != GEOS_POLYGON && type != GEOS_MULTIPOLYGON) {
219  int quadsegs = 1; // coarse
220  double tiny_distance = 0.000000001;
221  auto ng = GEOSBuffer_r(context, g, tiny_distance, quadsegs);
222  GEOSGeom_destroy_r(context, g);
223  return ng;
224  }
225  }
226  }
227  return g;
228 }
229 #endif
230 
231 extern "C" bool Geos_Wkb_Wkb(int op,
232  int arg1_type,
233  int8_t* arg1_coords,
234  int64_t arg1_coords_size,
235  int32_t* arg1_meta1,
236  int64_t arg1_meta1_size,
237  int32_t* arg1_meta2,
238  int64_t arg1_meta2_size,
239  // TODO: add meta3 args to support generic geometries
240  int32_t arg1_ic,
241  int32_t arg1_srid,
242  int arg2_type,
243  int8_t* arg2_coords,
244  int64_t arg2_coords_size,
245  int32_t* arg2_meta1,
246  int64_t arg2_meta1_size,
247  int32_t* arg2_meta2,
248  int64_t arg2_meta2_size,
249  // TODO: add meta3 args to support generic geometries
250  int32_t arg2_ic,
251  int32_t arg2_srid,
252  // TODO: add transform args
253  int* result_type,
254  int8_t** result_coords,
255  int64_t* result_coords_size,
256  int32_t** result_meta1,
257  int64_t* result_meta1_size,
258  int32_t** result_meta2,
259  int64_t* result_meta2_size) {
260 #ifndef __CUDACC__
261  // Get the result geo
262  // What if intersection is not a POLYGON? POINT? LINESTRING, MULTIPOLYGON?
263  // What if intersection is empty? Return null buffer pointers? Return false?
264  // What if geos fails?
265 
266  WKB wkb1{};
267  if (!toWkb(wkb1,
268  arg1_type,
269  arg1_coords,
270  arg1_coords_size,
271  arg1_meta1,
272  arg1_meta1_size,
273  arg1_meta2,
274  arg1_meta2_size,
275  arg1_ic)) {
276  return false;
277  }
278  WKB wkb2{};
279  if (!toWkb(wkb2,
280  arg2_type,
281  arg2_coords,
282  arg2_coords_size,
283  arg2_meta1,
284  arg2_meta1_size,
285  arg2_meta2,
286  arg2_meta2_size,
287  arg2_ic)) {
288  return false;
289  }
290  auto status = false;
291  auto context = create_context();
292  if (!context) {
293  return status;
294  }
295  auto* g1 = GEOSGeomFromWKB_buf_r(context, wkb1.data(), wkb1.size());
296  if (g1) {
297  GEOSSetSRID_r(context, g1, arg1_srid);
298  auto* g2 = GEOSGeomFromWKB_buf_r(context, wkb2.data(), wkb2.size());
299  if (g2) {
300  GEOSSetSRID_r(context, g2, arg2_srid);
301  GEOSGeometry* g = nullptr;
302  if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kINTERSECTION) {
303  g = GEOSIntersection_r(context, g1, g2);
304  } else if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kDIFFERENCE) {
305  g = GEOSDifference_r(context, g1, g2);
306  } else if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kUNION) {
307  g = GEOSUnion_r(context, g1, g2);
308  }
309  g = postprocess(context, g);
310  if (g) {
311  size_t wkb_size = 0ULL;
312  auto wkb_buf = GEOSGeomToWKB_buf_r(context, g, &wkb_size);
313  if (wkb_buf && wkb_size > 0ULL) {
314  WKB wkb(wkb_buf, wkb_buf + wkb_size);
315  free(wkb_buf);
316  status = fromWkb(wkb,
317  result_type,
318  result_coords,
319  result_coords_size,
320  result_meta1,
321  result_meta1_size,
322  result_meta2,
323  result_meta2_size);
324  }
325  GEOSGeom_destroy_r(context, g);
326  }
327  GEOSGeom_destroy_r(context, g2);
328  }
329  GEOSGeom_destroy_r(context, g1);
330  }
331  destroy_context(context);
332  return status;
333 #else
334  return false;
335 #endif
336 }
337 
338 extern "C" bool Geos_Wkb_double(int op,
339  int arg1_type,
340  int8_t* arg1_coords,
341  int64_t arg1_coords_size,
342  int32_t* arg1_meta1,
343  int64_t arg1_meta1_size,
344  int32_t* arg1_meta2,
345  int64_t arg1_meta2_size,
346  // TODO: add meta3 args to support generic geometries
347  int32_t arg1_ic,
348  int32_t arg1_srid,
349  double arg2,
350  // TODO: add transform args
351  int* result_type,
352  int8_t** result_coords,
353  int64_t* result_coords_size,
354  int32_t** result_meta1,
355  int64_t* result_meta1_size,
356  int32_t** result_meta2,
357  int64_t* result_meta2_size) {
358 #ifndef __CUDACC__
359  WKB wkb1{};
360  if (!toWkb(wkb1,
361  arg1_type,
362  arg1_coords,
363  arg1_coords_size,
364  arg1_meta1,
365  arg1_meta1_size,
366  arg1_meta2,
367  arg1_meta2_size,
368  arg1_ic)) {
369  return false;
370  }
371 
372  auto status = false;
373  auto context = create_context();
374  if (!context) {
375  return status;
376  }
377  auto* g1 = GEOSGeomFromWKB_buf_r(context, wkb1.data(), wkb1.size());
378  if (g1) {
379  GEOSSetSRID_r(context, g1, arg1_srid);
380  GEOSGeometry* g = nullptr;
381  if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kBUFFER) {
382  int quadsegs = 8; // default
383  g = GEOSBuffer_r(context, g1, arg2, quadsegs);
384  }
385  g = postprocess(context, g);
386  if (g) {
387  size_t wkb_size = 0ULL;
388  auto wkb_buf = GEOSGeomToWKB_buf_r(context, g, &wkb_size);
389  if (wkb_buf && wkb_size > 0ULL) {
390  WKB wkb(wkb_buf, wkb_buf + wkb_size);
391  free(wkb_buf);
392  status = fromWkb(wkb,
393  result_type,
394  result_coords,
395  result_coords_size,
396  result_meta1,
397  result_meta1_size,
398  result_meta2,
399  result_meta2_size);
400  }
401  GEOSGeom_destroy_r(context, g);
402  }
403  GEOSGeom_destroy_r(context, g1);
404  }
405  destroy_context(context);
406  return status;
407 #else
408  return false;
409 #endif
410 }
411 
412 extern "C" bool Geos_Wkb(int op,
413  int arg_type,
414  int8_t* arg_coords,
415  int64_t arg_coords_size,
416  int32_t* arg_meta1,
417  int64_t arg_meta1_size,
418  int32_t* arg_meta2,
419  int64_t arg_meta2_size,
420  // TODO: add meta3 args to support generic geometries
421  int32_t arg_ic,
422  int32_t arg_srid,
423  bool* result) {
424 #ifndef __CUDACC__
425  WKB wkb1{};
426  if (!result || !toWkb(wkb1,
427  arg_type,
428  arg_coords,
429  arg_coords_size,
430  arg_meta1,
431  arg_meta1_size,
432  arg_meta2,
433  arg_meta2_size,
434  arg_ic)) {
435  return false;
436  }
437 
438  auto status = false;
439  auto context = create_context();
440  if (!context) {
441  return status;
442  }
443  auto* g1 = GEOSGeomFromWKB_buf_r(context, wkb1.data(), wkb1.size());
444  if (g1) {
445  GEOSSetSRID_r(context, g1, arg_srid);
446  if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kISEMPTY) {
447  *result = GEOSisEmpty_r(context, g1);
448  status = true;
449  } else if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kISVALID) {
450  *result = GEOSisValid_r(context, g1);
451  status = true;
452  }
453  GEOSGeom_destroy_r(context, g1);
454  }
455  destroy_context(context);
456  return status;
457 #else
458  return false;
459 #endif
460 }
461 
462 #endif
static std::unique_ptr< GeoBase > createGeoType(const std::string &wkt_or_wkb_hex)
Definition: Types.cpp:683
#define LOG(tag)
Definition: Logger.h:188
bool Geos_Wkb(int op, int arg_type, int8_t *arg_coords, int64_t arg_coords_size, int32_t *arg_meta1, int64_t arg_meta1_size, int32_t *arg_meta2, int64_t arg_meta2_size, int32_t arg_ic, int32_t arg_srid, bool *result)
void * checked_malloc(const size_t size)
Definition: checked_alloc.h:44
bool Geos_Wkb_Wkb(int op, int arg1_type, int8_t *arg1_coords, int64_t arg1_coords_size, int32_t *arg1_meta1, int64_t arg1_meta1_size, int32_t *arg1_meta2, int64_t arg1_meta2_size, int32_t arg1_ic, int32_t arg1_srid, int arg2_type, int8_t *arg2_coords, int64_t arg2_coords_size, int32_t *arg2_meta1, int64_t arg2_meta1_size, int32_t *arg2_meta2, int64_t arg2_meta2_size, int32_t arg2_ic, int32_t arg2_srid, int *result_type, int8_t **result_coords, int64_t *result_coords_size, int32_t **result_meta1, int64_t *result_meta1_size, int32_t **result_meta2, int64_t *result_meta2_size)
bool Geos_Wkb_double(int op, int arg1_type, int8_t *arg1_coords, int64_t arg1_coords_size, int32_t *arg1_meta1, int64_t arg1_meta1_size, int32_t *arg1_meta2, int64_t arg1_meta2_size, int32_t arg1_ic, int32_t arg1_srid, double arg2, int *result_type, int8_t **result_coords, int64_t *result_coords_size, int32_t **result_meta1, int64_t *result_meta1_size, int32_t **result_meta2, int64_t *result_meta2_size)
std::shared_ptr< std::vector< double > > decompress_coords< double, int32_t >(const int32_t &ic, const int8_t *coords, const size_t coords_sz)
#define CHECK(condition)
Definition: Logger.h:197