OmniSciDB  8fa3bf436f
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  int32_t* best_planar_srid_ptr) {
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  if (best_planar_srid_ptr) {
106  // A non-NULL pointer signifies a request to find the best planar srid
107  // to transform this WGS84 geometry to.
108  *best_planar_srid_ptr = point.getBestPlanarSRID();
109  if (!point.transform(4326, *best_planar_srid_ptr)) {
110  return false;
111  }
112  }
113  return point.getWkb(wkb);
114  }
115  if (static_cast<SQLTypes>(type) == kLINESTRING) {
116  GeoLineString linestring(*cv);
117  if (best_planar_srid_ptr) {
118  // A non-NULL pointer signifies a request to find the best planar srid
119  // to transform this WGS84 geometry to, based on geometry's centroid.
120  *best_planar_srid_ptr = linestring.getBestPlanarSRID();
121  if (!linestring.transform(4326, *best_planar_srid_ptr)) {
122  return false;
123  }
124  }
125  return linestring.getWkb(wkb);
126  }
127  std::vector<int32_t> meta1v(meta1, meta1 + meta1_size);
128  if (static_cast<SQLTypes>(type) == kPOLYGON) {
129  GeoPolygon poly(*cv, meta1v);
130  if (best_planar_srid_ptr) {
131  // A non-NULL pointer signifies a request to find the best planar srid
132  // to transform this WGS84 geometry to, based on geometry's centroid.
133  *best_planar_srid_ptr = poly.getBestPlanarSRID();
134  if (!poly.transform(4326, *best_planar_srid_ptr)) {
135  return false;
136  };
137  }
138  return poly.getWkb(wkb);
139  }
140  std::vector<int32_t> meta2v(meta2, meta2 + meta2_size);
141  if (static_cast<SQLTypes>(type) == kMULTIPOLYGON) {
142  // Recognize GEOMETRYCOLLECTION EMPTY encoding
143  // MULTIPOLYGON (((0 0,0.00000012345 0.0,0.0 0.00000012345,0 0)))
144  // Used to pass along EMPTY from ST_Intersection to ST_IsEmpty for example
145  if (meta1_size == 1 && meta2_size == 1) {
146  const std::vector<double> ecv = {0.0, 0.0, 0.00000012345, 0.0, 0.0, 0.00000012345};
147  if (*cv == ecv) {
148  GeoGeometryCollection empty("GEOMETRYCOLLECTION EMPTY");
149  return empty.getWkb(wkb);
150  }
151  }
152  GeoMultiPolygon mpoly(*cv, meta1v, meta2v);
153  if (best_planar_srid_ptr) {
154  // A non-NULL pointer signifies a request to find the best planar srid
155  // to transform this WGS84 geometry to, based on geometry's centroid.
156  *best_planar_srid_ptr = mpoly.getBestPlanarSRID();
157  if (!mpoly.transform(4326, *best_planar_srid_ptr)) {
158  return false;
159  };
160  }
161  return mpoly.getWkb(wkb);
162  }
163  return false;
164 }
165 
166 // Conversion form wkb to internal vector representation.
167 // Each vector components is malloced, caller is reponsible for freeing.
168 bool fromWkb(WKB& wkb,
169  int* result_type,
170  int8_t** result_coords,
171  int64_t* result_coords_size,
172  int32_t** result_meta1,
173  int64_t* result_meta1_size,
174  int32_t** result_meta2,
175  int64_t* result_meta2_size,
176  int32_t* best_planar_srid_ptr) {
178  if (best_planar_srid_ptr && !result->isEmpty()) {
179  // If original geometry has previously been projected to planar srid,
180  // need to transform back to WGS84
181  if (!result->transform(*best_planar_srid_ptr, 4326)) {
182  return false;
183  }
184  }
185 
186  // Get the column values
187  std::vector<double> coords{};
188  std::vector<int32_t> ring_sizes{};
189  std::vector<int32_t> poly_rings{};
190  std::vector<double> bounds{};
191 
192  // Forcing MULTIPOLYGON result until we can handle any geo.
193  if (result->isEmpty()) {
194  // Generate a tiny polygon around POINT(0 0), make it a multipolygon
195  // MULTIPOLYGON (((0 0,0.00000012345 0.0,0.0 0.00000012345,0 0)))
196  // to simulate an empty result
197  coords = {0.0, 0.0, 0.00000012345, 0.0, 0.0, 0.00000012345};
198  ring_sizes.push_back(3);
199  poly_rings.push_back(1);
200  } else if (auto result_point = dynamic_cast<GeoPoint*>(result.get())) {
201  result_point->getColumns(coords);
202  // Generate a tiny polygon around the point, make it a multipolygon
203  coords.push_back(coords[0] + 0.0000001);
204  coords.push_back(coords[1]);
205  coords.push_back(coords[0]);
206  coords.push_back(coords[1] + 0.0000001);
207  ring_sizes.push_back(3);
208  poly_rings.push_back(ring_sizes.size());
209  } else if (auto result_poly = dynamic_cast<GeoPolygon*>(result.get())) {
210  result_poly->getColumns(coords, ring_sizes, bounds);
211  // Convert to a 1-polygon multipolygon
212  poly_rings.push_back(ring_sizes.size());
213  } else if (auto result_mpoly = dynamic_cast<GeoMultiPolygon*>(result.get())) {
214  result_mpoly->getColumns(coords, ring_sizes, poly_rings, bounds);
215  } else {
216  return false;
217  }
218 
219  // TODO: consider using a single buffer to hold all components,
220  // instead of allocating and registering each component buffer separately
221 
222  *result_type = static_cast<int>(kMULTIPOLYGON);
223 
224  *result_coords = nullptr;
225  int64_t size = coords.size() * sizeof(double);
226  if (size > 0) {
227  auto buf = checked_malloc(size);
228  std::memcpy(buf, coords.data(), size);
229  *result_coords = reinterpret_cast<int8_t*>(buf);
230  }
231  *result_coords_size = size;
232 
233  *result_meta1 = nullptr;
234  size = ring_sizes.size() * sizeof(int32_t);
235  if (size > 0) {
236  auto buf = checked_malloc(size);
237  std::memcpy(buf, ring_sizes.data(), size);
238  *result_meta1 = reinterpret_cast<int32_t*>(buf);
239  }
240  *result_meta1_size = ring_sizes.size();
241 
242  *result_meta2 = nullptr;
243  size = poly_rings.size() * sizeof(int32_t);
244  if (size > 0) {
245  auto buf = checked_malloc(size);
246  std::memcpy(buf, poly_rings.data(), size);
247  *result_meta2 = reinterpret_cast<int32_t*>(buf);
248  }
249  *result_meta2_size = poly_rings.size();
250 
251  return true;
252 }
253 
254 GEOSGeometry* postprocess(GEOSContextHandle_t context, GEOSGeometry* g) {
255  if (g && GEOSisEmpty_r(context, g) == 0) {
256  auto type = GEOSGeomTypeId_r(context, g);
257  if (type != -1) {
258  if (type != GEOS_POINT && type != GEOS_POLYGON && type != GEOS_MULTIPOLYGON) {
259  int quadsegs = 1; // coarse
260  double tiny_distance = 0.000000001;
261  auto ng = GEOSBuffer_r(context, g, tiny_distance, quadsegs);
262  GEOSGeom_destroy_r(context, g);
263  return ng;
264  }
265  }
266  }
267  return g;
268 }
269 #endif
270 
271 extern "C" RUNTIME_EXPORT bool Geos_Wkb_Wkb(
272  int op,
273  int arg1_type,
274  int8_t* arg1_coords,
275  int64_t arg1_coords_size,
276  int32_t* arg1_meta1,
277  int64_t arg1_meta1_size,
278  int32_t* arg1_meta2,
279  int64_t arg1_meta2_size,
280  // TODO: add meta3 args to support generic geometries
281  int32_t arg1_ic,
282  int32_t arg1_srid,
283  int arg2_type,
284  int8_t* arg2_coords,
285  int64_t arg2_coords_size,
286  int32_t* arg2_meta1,
287  int64_t arg2_meta1_size,
288  int32_t* arg2_meta2,
289  int64_t arg2_meta2_size,
290  // TODO: add meta3 args to support generic geometries
291  int32_t arg2_ic,
292  int32_t arg2_srid,
293  // TODO: add transform args
294  int* result_type,
295  int8_t** result_coords,
296  int64_t* result_coords_size,
297  int32_t** result_meta1,
298  int64_t* result_meta1_size,
299  int32_t** result_meta2,
300  int64_t* result_meta2_size) {
301 #ifndef __CUDACC__
302  // Get the result geo
303  // What if intersection is not a POLYGON? POINT? LINESTRING, MULTIPOLYGON?
304  // What if intersection is empty? Return null buffer pointers? Return false?
305  // What if geos fails?
306 
307  WKB wkb1{};
308  if (!toWkb(wkb1,
309  arg1_type,
310  arg1_coords,
311  arg1_coords_size,
312  arg1_meta1,
313  arg1_meta1_size,
314  arg1_meta2,
315  arg1_meta2_size,
316  arg1_ic,
317  nullptr)) {
318  return false;
319  }
320  WKB wkb2{};
321  if (!toWkb(wkb2,
322  arg2_type,
323  arg2_coords,
324  arg2_coords_size,
325  arg2_meta1,
326  arg2_meta1_size,
327  arg2_meta2,
328  arg2_meta2_size,
329  arg2_ic,
330  nullptr)) {
331  return false;
332  }
333  auto status = false;
334  auto context = create_context();
335  if (!context) {
336  return status;
337  }
338  auto* g1 = GEOSGeomFromWKB_buf_r(context, wkb1.data(), wkb1.size());
339  if (g1) {
340  auto* g2 = GEOSGeomFromWKB_buf_r(context, wkb2.data(), wkb2.size());
341  if (g2) {
342  GEOSGeometry* g = nullptr;
343  if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kINTERSECTION) {
344  g = GEOSIntersection_r(context, g1, g2);
345  } else if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kDIFFERENCE) {
346  g = GEOSDifference_r(context, g1, g2);
347  } else if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kUNION) {
348  g = GEOSUnion_r(context, g1, g2);
349  }
350  g = postprocess(context, g);
351  if (g) {
352  size_t wkb_size = 0ULL;
353  auto wkb_buf = GEOSGeomToWKB_buf_r(context, g, &wkb_size);
354  if (wkb_buf && wkb_size > 0ULL) {
355  WKB wkb(wkb_buf, wkb_buf + wkb_size);
356  free(wkb_buf);
357  status = fromWkb(wkb,
358  result_type,
359  result_coords,
360  result_coords_size,
361  result_meta1,
362  result_meta1_size,
363  result_meta2,
364  result_meta2_size,
365  nullptr);
366  }
367  GEOSGeom_destroy_r(context, g);
368  }
369  GEOSGeom_destroy_r(context, g2);
370  }
371  GEOSGeom_destroy_r(context, g1);
372  }
373  destroy_context(context);
374  return status;
375 #else
376  return false;
377 #endif
378 }
379 
380 extern "C" RUNTIME_EXPORT bool Geos_Wkb_double(
381  int op,
382  int arg1_type,
383  int8_t* arg1_coords,
384  int64_t arg1_coords_size,
385  int32_t* arg1_meta1,
386  int64_t arg1_meta1_size,
387  int32_t* arg1_meta2,
388  int64_t arg1_meta2_size,
389  // TODO: add meta3 args to support generic geometries
390  int32_t arg1_ic,
391  int32_t arg1_srid,
392  double arg2,
393  // TODO: add transform args
394  int* result_type,
395  int8_t** result_coords,
396  int64_t* result_coords_size,
397  int32_t** result_meta1,
398  int64_t* result_meta1_size,
399  int32_t** result_meta2,
400  int64_t* result_meta2_size) {
401 #ifndef __CUDACC__
402  int32_t best_planar_srid;
403  int32_t* best_planar_srid_ptr = nullptr;
404  if (arg1_srid == 4326 && static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kBUFFER) {
405  // Use the best (location-based) planar transform to project geometry before
406  // running geos operation, back-project the result to 4326
407  best_planar_srid_ptr = &best_planar_srid;
408  }
409  WKB wkb1{};
410  // Project to best planar srid before running certain geos ops
411  if (!toWkb(wkb1,
412  arg1_type,
413  arg1_coords,
414  arg1_coords_size,
415  arg1_meta1,
416  arg1_meta1_size,
417  arg1_meta2,
418  arg1_meta2_size,
419  arg1_ic,
420  best_planar_srid_ptr)) {
421  return false;
422  }
423 
424  auto status = false;
425  auto context = create_context();
426  if (!context) {
427  return status;
428  }
429  auto* g1 = GEOSGeomFromWKB_buf_r(context, wkb1.data(), wkb1.size());
430  if (g1) {
431  GEOSGeometry* g = nullptr;
432  if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kBUFFER) {
433  int quadsegs = 8; // default
434  g = GEOSBuffer_r(context, g1, arg2, quadsegs);
435  }
436  g = postprocess(context, g);
437  if (g) {
438  size_t wkb_size = 0ULL;
439  auto wkb_buf = GEOSGeomToWKB_buf_r(context, g, &wkb_size);
440  if (wkb_buf && wkb_size > 0ULL) {
441  WKB wkb(wkb_buf, wkb_buf + wkb_size);
442  free(wkb_buf);
443  // Back-project the result from planar to 4326 if necessary
444  status = fromWkb(wkb,
445  result_type,
446  result_coords,
447  result_coords_size,
448  result_meta1,
449  result_meta1_size,
450  result_meta2,
451  result_meta2_size,
452  best_planar_srid_ptr);
453  }
454  GEOSGeom_destroy_r(context, g);
455  }
456  GEOSGeom_destroy_r(context, g1);
457  }
458  destroy_context(context);
459  return status;
460 #else
461  return false;
462 #endif
463 }
464 
465 extern "C" RUNTIME_EXPORT bool Geos_Wkb(
466  int op,
467  int arg_type,
468  int8_t* arg_coords,
469  int64_t arg_coords_size,
470  int32_t* arg_meta1,
471  int64_t arg_meta1_size,
472  int32_t* arg_meta2,
473  int64_t arg_meta2_size,
474  // TODO: add meta3 args to support generic geometries
475  int32_t arg_ic,
476  int32_t arg_srid,
477  bool* result) {
478 #ifndef __CUDACC__
479  WKB wkb1{};
480  if (!result || !toWkb(wkb1,
481  arg_type,
482  arg_coords,
483  arg_coords_size,
484  arg_meta1,
485  arg_meta1_size,
486  arg_meta2,
487  arg_meta2_size,
488  arg_ic,
489  nullptr)) {
490  return false;
491  }
492 
493  auto status = false;
494  auto context = create_context();
495  if (!context) {
496  return status;
497  }
498  auto* g1 = GEOSGeomFromWKB_buf_r(context, wkb1.data(), wkb1.size());
499  if (g1) {
500  if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kISEMPTY) {
501  *result = GEOSisEmpty_r(context, g1);
502  status = true;
503  } else if (static_cast<GeoBase::GeoOp>(op) == GeoBase::GeoOp::kISVALID) {
504  *result = GEOSisValid_r(context, g1);
505  status = true;
506  }
507  GEOSGeom_destroy_r(context, g1);
508  }
509  destroy_context(context);
510  return status;
511 #else
512  return false;
513 #endif
514 }
515 
516 #endif
static std::unique_ptr< GeoBase > createGeoType(const std::string &wkt_or_wkb_hex)
Definition: Types.cpp:889
RUNTIME_EXPORT 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)
#define LOG(tag)
Definition: Logger.h:194
RUNTIME_EXPORT 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)
RUNTIME_EXPORT 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)
void * checked_malloc(const size_t size)
Definition: checked_alloc.h:45
#define RUNTIME_EXPORT
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:203