OmniSciDB  2e3a973ef4
Types.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2017 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 #include "Geospatial/Types.h"
18 
19 #include <limits>
20 
21 #include <gdal.h>
22 #include <ogr_geometry.h>
23 #include <ogrsf_frmts.h>
24 
25 #include "Logger/Logger.h"
26 #include "Shared/sqltypes.h"
27 
37 namespace {
38 constexpr auto DOUBLE_MAX = std::numeric_limits<double>::max();
39 constexpr auto DOUBLE_MIN = std::numeric_limits<double>::lowest();
40 
41 struct Coords {
42  double x;
43  double y;
44  Coords(double x, double y) : x(x), y(y) {}
45 };
46 
47 struct BoundingBox {
50  BoundingBox() : min(DOUBLE_MAX, DOUBLE_MAX), max(DOUBLE_MIN, DOUBLE_MIN) {}
51 
52  void update(double x, double y) {
53  if (x < min.x) {
54  min.x = x;
55  }
56  if (y < min.y) {
57  min.y = y;
58  }
59  if (x > max.x) {
60  max.x = x;
61  }
62  if (y > max.y) {
63  max.y = y;
64  }
65  }
66 };
67 
68 int process_poly_ring(OGRLinearRing* ring,
69  std::vector<double>& coords,
70  BoundingBox* bbox) {
71  double last_x = DOUBLE_MAX, last_y = DOUBLE_MAX;
72  size_t first_index = coords.size();
73  int num_points_added = 0;
74  int num_points_in_ring = ring->getNumPoints();
75  if (num_points_in_ring < 3) {
77  "PolyRing",
78  "All poly rings must have more than 3 points. Found ring with " +
79  std::to_string(num_points_in_ring) + " points.");
80  }
81  for (auto i = 0; i < num_points_in_ring; i++) {
82  OGRPoint point;
83  ring->getPoint(i, &point);
84  last_x = point.getX();
85  last_y = point.getY();
86  coords.push_back(last_x);
87  coords.push_back(last_y);
88  if (bbox) {
89  bbox->update(last_x, last_y);
90  }
91  num_points_added++;
92  }
93  // Store all rings as open rings (implicitly assumes all rings are closed)
94  if ((coords[first_index] == last_x) && (coords[first_index + 1] == last_y)) {
95  coords.pop_back();
96  coords.pop_back();
97  num_points_added--;
98  if (num_points_added < 3) {
100  "PolyRing",
101  "All exterior rings must have more than 3 points. Found ring with " +
102  std::to_string(num_points_added) + " points.");
103  }
104  }
105  return num_points_added;
106 }
107 
108 } // namespace
109 
110 namespace Geospatial {
111 
112 std::string GeoTypesError::OGRErrorToStr(const int ogr_err) {
113  switch (ogr_err) {
114  case OGRERR_NOT_ENOUGH_DATA:
115  return std::string("not enough input data");
116  case OGRERR_NOT_ENOUGH_MEMORY:
117  return std::string("not enough memory");
118  case OGRERR_UNSUPPORTED_GEOMETRY_TYPE:
119  return std::string("unsupported geometry type");
120  case OGRERR_UNSUPPORTED_OPERATION:
121  return std::string("unsupported operation");
122  case OGRERR_CORRUPT_DATA:
123  return std::string("corrupt input data");
124  case OGRERR_FAILURE:
125  return std::string("ogr failure");
126  case OGRERR_UNSUPPORTED_SRS:
127  return std::string("unsupported spatial reference system");
128  case OGRERR_INVALID_HANDLE:
129  return std::string("invalid file handle");
130 #if (GDAL_VERSION_MAJOR > 1)
131  case OGRERR_NON_EXISTING_FEATURE:
132  return std::string("feature does not exist in input geometry");
133 #endif
134  default:
135  return std::string("Unknown OGOR error encountered: ") + std::to_string(ogr_err);
136  }
137 }
138 
139 GeoBase::~GeoBase() {
140  // Note: Removing the geometry object that was pulled from an OGRFeature results in a
141  // segfault. If we are wrapping around a pre-existing OGRGeometry object, we let the
142  // caller manage the memory.
143  if (geom_ && owns_geom_obj_) {
144  OGRGeometryFactory::destroyGeometry(geom_);
145  }
146 }
147 
148 OGRErr GeoBase::createFromWktString(const std::string& wkt, OGRGeometry** geom) {
149 #if (GDAL_VERSION_MAJOR > 2) || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 3)
150  OGRErr ogr_status = OGRGeometryFactory::createFromWkt(wkt.c_str(), nullptr, geom);
151 #else
152  auto data = (char*)wkt.c_str();
153  OGRErr ogr_status = OGRGeometryFactory::createFromWkt(&data, NULL, geom);
154 #endif
155  return ogr_status;
156 }
157 
158 std::string GeoBase::getWktString() const {
159  char* wkt = nullptr;
160  geom_->exportToWkt(&wkt);
161  CHECK(wkt);
162  std::string wkt_str(wkt);
163  CPLFree(wkt);
164  return wkt_str;
165 }
166 
167 OGRErr GeoBase::createFromWkb(const std::vector<uint8_t>& wkb, OGRGeometry** geom) {
168 #if (GDAL_VERSION_MAJOR > 2) || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 3)
169  OGRErr ogr_status =
170  OGRGeometryFactory::createFromWkb(wkb.data(), nullptr, geom, wkb.size());
171  return ogr_status;
172 #else
173  CHECK(false);
174 #endif
175 }
176 
177 // GeoBase could also generate GEOS geometries through geom_->exportToGEOS(),
178 // conversion to WKB and subsequent load of WKB into GEOS could be eliminated.
179 
180 bool GeoBase::getWkb(std::vector<uint8_t>& wkb) const {
181  auto size = geom_->WkbSize();
182  if (size > 0) {
183  wkb.resize(size);
184  geom_->exportToWkb(wkbNDR, wkb.data());
185  return true;
186  }
187  return false;
188 }
189 
190 bool GeoBase::isEmpty() const {
191  return geom_ && geom_->IsEmpty();
192 }
193 
194 bool GeoBase::operator==(const GeoBase& other) const {
195  if (!this->geom_ || !other.geom_) {
196  return false;
197  }
198  return this->geom_->Equals(other.geom_);
199 }
200 
201 // Run a specific geo operation on this and other geometries
202 std::unique_ptr<GeoBase> GeoBase::run(GeoBase::GeoOp op, const GeoBase& other) const {
203  OGRGeometry* result = nullptr;
204  // Invalid geometries are derailing geo op performance,
205  // Checking validity before running an operation doesn't have lesser penalty.
206  // DOES NOT HELP: if (geom_->IsValid() && other.geom_->IsValid()) {
207  switch (op) {
208  case GeoBase::GeoOp::kINTERSECTION:
209  result = geom_->Intersection(other.geom_);
210  break;
211  case GeoBase::GeoOp::kDIFFERENCE:
212  result = geom_->Difference(other.geom_);
213  break;
214  case GeoBase::GeoOp::kUNION:
215  result = geom_->Union(other.geom_);
216  break;
217  default:
218  break;
219  }
220  // TODO: Need to handle empty/non-POLYGON result
221  if (!result || result->IsEmpty() ||
222  !(result->getGeometryType() == wkbPolygon ||
223  result->getGeometryType() == wkbMultiPolygon)) {
224  // throw GeoTypesError(std::string(OGRGeometryTypeToName(geom_->getGeometryType())),
225  // "Currenly don't support invalid or empty result");
226  // return GeoTypesFactory::createGeoType("POLYGON EMPTY");
227  // Not supporting EMPTY polygons, return a dot polygon
228  return GeoTypesFactory::createGeoType(
229  "MULTIPOLYGON(((0 0,0.0000001 0,0 0.0000001)))");
230  }
231  return GeoTypesFactory::createGeoType(result);
232 }
233 
234 // Run a specific geo operation on this and other geometries
235 std::unique_ptr<GeoBase> GeoBase::optimized_run(GeoBase::GeoOp op,
236  const GeoBase& other) const {
237  OGRGeometry* result = nullptr;
238  // Loop through polys combinations, check validity, do intersections
239  // where needed, return a union of all intersections
240  auto gc1 = geom_->toGeometryCollection();
241  auto gc2 = other.geom_->toGeometryCollection();
242  if (!gc1 || !gc2 || gc1->IsEmpty() || gc2->IsEmpty()) {
243  return nullptr;
244  }
245  for (int i1 = 0; i1 < gc1->getNumGeometries(); i1++) {
246  auto g1 = gc1->getGeometryRef(i1);
247  // Validity check is very slow
248  if (!g1 || g1->IsEmpty() /*|| !g1->IsValid()*/) {
249  continue;
250  }
251  OGREnvelope ge1;
252  g1->getEnvelope(&ge1);
253  for (int i2 = 0; i2 < gc2->getNumGeometries(); i2++) {
254  auto g2 = gc2->getGeometryRef(i2);
255  // Validity check is very slow
256  if (!g2 || g2->IsEmpty() /*|| !g2->IsValid()*/) {
257  continue;
258  }
259  // Check for bounding box overlap
260  OGREnvelope ge2;
261  g2->getEnvelope(&ge2);
262  if (!ge1.Intersects(ge2)) {
263  continue;
264  }
265  // Do intersection
266  auto intermediate_result = g1->Intersection(g2);
267  if (!intermediate_result || intermediate_result->IsEmpty()) {
268  continue;
269  }
270  if (!result) {
271  result = intermediate_result;
272  } else {
273  result = result->Union(intermediate_result);
274  }
275  }
276  }
277 
278  // TODO: Need to handle empty/non-POLYGON result
279  if (!result || result->IsEmpty() ||
280  !(result->getGeometryType() == wkbPolygon ||
281  result->getGeometryType() == wkbMultiPolygon)) {
282  // throw GeoTypesError(std::string(OGRGeometryTypeToName(geom_->getGeometryType())),
283  // "Currenly don't support invalid or empty result");
284  // return GeoTypesFactory::createGeoType("POLYGON EMPTY");
285  // Not supporting EMPTY polygons, return a dot polygon
286  return GeoTypesFactory::createGeoType(
287  "MULTIPOLYGON(((0 0,0.0000001 0,0 0.0000001)))");
288  }
289  return GeoTypesFactory::createGeoType(result);
290 }
291 
292 // Run a specific geo operation on this geometry and a double param
293 std::unique_ptr<GeoBase> GeoBase::run(GeoBase::GeoOp op, double param) const {
294  OGRGeometry* result = nullptr;
295  switch (op) {
296  case GeoBase::GeoOp::kBUFFER:
297  result = geom_->Buffer(param);
298  break;
299  default:
300  break;
301  }
302  // TODO: Need to handle empty/non-POLYGON result
303  if (!result || result->IsEmpty() ||
304  !(result->getGeometryType() == wkbPolygon ||
305  result->getGeometryType() == wkbMultiPolygon)) {
306  // throw GeoTypesError(std::string(OGRGeometryTypeToName(geom_->getGeometryType())),
307  // "Currenly don't support invalid or empty result");
308  // return GeoTypesFactory::createGeoType("POLYGON EMPTY");
309  // Not supporting EMPTY polygons, return a dot polygon
310  return GeoTypesFactory::createGeoType(
311  "MULTIPOLYGON(((0 0,0.0000001 0,0 0.0000001)))");
312  }
313  return GeoTypesFactory::createGeoType(result);
314 }
315 
316 // Run a specific predicate operation on this geometry
317 bool GeoBase::run(GeoBase::GeoOp op) const {
318  auto result = false;
319  switch (op) {
320  case GeoBase::GeoOp::kISVALID:
321  result = geom_->IsValid();
322  break;
323  case GeoBase::GeoOp::kISEMPTY:
324  result = isEmpty();
325  break;
326  default:
327  break;
328  }
329  return result;
330 }
331 
332 GeoPoint::GeoPoint(const std::vector<double>& coords) {
333  if (coords.size() != 2) {
334  throw GeoTypesError("Point",
335  "Incorrect coord size of " + std::to_string(coords.size()) +
336  " supplied. Expected 2.");
337  }
338  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbPoint);
339  OGRPoint* point = dynamic_cast<OGRPoint*>(geom_);
340  CHECK(point);
341  point->setX(coords[0]);
342  point->setY(coords[1]);
343 }
344 
345 GeoPoint::GeoPoint(const std::string& wkt) {
346  const auto err = GeoBase::createFromWktString(wkt, &geom_);
347  if (err != OGRERR_NONE) {
348  throw GeoTypesError("Point", err);
349  }
350  CHECK(geom_);
351  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbPoint) {
352  throw GeoTypesError("Point",
353  "Unexpected geometry type from WKT string: " +
354  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
355  }
356 }
357 
358 void GeoPoint::getColumns(std::vector<double>& coords) const {
359  const auto point_geom = dynamic_cast<OGRPoint*>(geom_);
360  CHECK(point_geom);
361 
362  if (point_geom->IsEmpty()) {
363  // until the run-time can handle empties
364  throw GeoTypesError("Point", "'EMPTY' not supported");
365  // we cannot yet handle NULL fixed-length array
366  // so we have to store sentinel values instead
367  coords.push_back(NULL_DOUBLE);
368  coords.push_back(NULL_DOUBLE);
369  return;
370  }
371 
372  coords.push_back(point_geom->getX());
373  coords.push_back(point_geom->getY());
374 }
375 
376 GeoLineString::GeoLineString(const std::vector<double>& coords) {
377  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbLineString);
378  OGRLineString* line = dynamic_cast<OGRLineString*>(geom_);
379  CHECK(line);
380  for (size_t i = 0; i < coords.size(); i += 2) {
381  line->addPoint(coords[i], coords[i + 1]);
382  }
383 }
384 
385 GeoLineString::GeoLineString(const std::string& wkt) {
386  const auto err = GeoBase::createFromWktString(wkt, &geom_);
387  if (err != OGRERR_NONE) {
388  throw GeoTypesError("LineString", err);
389  }
390  CHECK(geom_);
391  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbLineString) {
392  throw GeoTypesError("LineString",
393  "Unexpected geometry type from WKT string: " +
394  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
395  }
396 }
397 
398 void GeoLineString::getColumns(std::vector<double>& coords,
399  std::vector<double>& bounds) const {
400  auto linestring_geom = dynamic_cast<OGRLineString*>(geom_);
401  CHECK(linestring_geom);
402 
403  if (linestring_geom->IsEmpty()) {
404  // until the run-time can handle empties
405  throw GeoTypesError("LineString", "'EMPTY' not supported");
406  // return null bounds
407  bounds.push_back(NULL_DOUBLE);
408  bounds.push_back(NULL_DOUBLE);
409  bounds.push_back(NULL_DOUBLE);
410  bounds.push_back(NULL_DOUBLE);
411  return;
412  }
413 
414  BoundingBox bbox;
415  for (auto i = 0; i < linestring_geom->getNumPoints(); i++) {
416  OGRPoint point;
417  linestring_geom->getPoint(i, &point);
418  double x = point.getX();
419  double y = point.getY();
420  coords.push_back(x);
421  coords.push_back(y);
422  bbox.update(x, y);
423  }
424  bounds.push_back(bbox.min.x);
425  bounds.push_back(bbox.min.y);
426  bounds.push_back(bbox.max.x);
427  bounds.push_back(bbox.max.y);
428 }
429 
430 GeoPolygon::GeoPolygon(const std::vector<double>& coords,
431  const std::vector<int32_t>& ring_sizes) {
432  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbPolygon);
433  OGRPolygon* poly = dynamic_cast<OGRPolygon*>(geom_);
434  CHECK(poly);
435 
436  size_t coords_ctr = 0;
437  for (size_t r = 0; r < ring_sizes.size(); r++) {
438  OGRLinearRing ring;
439  const auto ring_sz = ring_sizes[r];
440  for (auto i = 0; i < 2 * ring_sz; i += 2) {
441  ring.addPoint(coords[coords_ctr + i], coords[coords_ctr + i + 1]);
442  }
443  ring.addPoint(coords[coords_ctr], coords[coords_ctr + 1]);
444  coords_ctr += 2 * ring_sz;
445  poly->addRing(&ring);
446  }
447 }
448 
449 GeoPolygon::GeoPolygon(const std::string& wkt) {
450  const auto err = GeoBase::createFromWktString(wkt, &geom_);
451  if (err != OGRERR_NONE) {
452  throw GeoTypesError("Polygon", err);
453  }
454  CHECK(geom_);
455  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbPolygon) {
456  throw GeoTypesError("Polygon",
457  "Unexpected geometry type from WKT string: " +
458  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
459  }
460 }
461 
462 void GeoPolygon::getColumns(std::vector<double>& coords,
463  std::vector<int32_t>& ring_sizes,
464  std::vector<double>& bounds) const {
465  const auto poly_geom = dynamic_cast<OGRPolygon*>(geom_);
466  CHECK(poly_geom);
467 
468  if (poly_geom->IsEmpty()) {
469  // until the run-time can handle empties
470  throw GeoTypesError("Polygon", "'EMPTY' not supported");
471  // return null bounds
472  bounds.push_back(NULL_DOUBLE);
473  bounds.push_back(NULL_DOUBLE);
474  bounds.push_back(NULL_DOUBLE);
475  bounds.push_back(NULL_DOUBLE);
476  return;
477  }
478 
479  BoundingBox bbox;
480  const auto exterior_ring = poly_geom->getExteriorRing();
481  CHECK(exterior_ring);
482  // All exterior rings are imported CCW
483  if (exterior_ring->isClockwise()) {
484  exterior_ring->reverseWindingOrder();
485  }
486  const auto num_points_added = process_poly_ring(exterior_ring, coords, &bbox);
487  ring_sizes.push_back(num_points_added);
488  for (auto r = 0; r < poly_geom->getNumInteriorRings(); r++) {
489  auto interior_ring = poly_geom->getInteriorRing(r);
490  CHECK(interior_ring);
491  // All interior rings are imported CW
492  if (!interior_ring->isClockwise()) {
493  interior_ring->reverseWindingOrder();
494  }
495  const auto num_points_added = process_poly_ring(interior_ring, coords, nullptr);
496  ring_sizes.push_back(num_points_added);
497  }
498  bounds.push_back(bbox.min.x);
499  bounds.push_back(bbox.min.y);
500  bounds.push_back(bbox.max.x);
501  bounds.push_back(bbox.max.y);
502 }
503 
504 int32_t GeoPolygon::getNumInteriorRings() const {
505  const auto poly_geom = dynamic_cast<OGRPolygon*>(geom_);
506  CHECK(poly_geom);
507  return poly_geom->getNumInteriorRings();
508 }
509 
510 GeoMultiPolygon::GeoMultiPolygon(const std::vector<double>& coords,
511  const std::vector<int32_t>& ring_sizes,
512  const std::vector<int32_t>& poly_rings) {
513  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbMultiPolygon);
514  OGRMultiPolygon* multipoly = dynamic_cast<OGRMultiPolygon*>(geom_);
515  CHECK(multipoly);
516 
517  size_t ring_ctr = 0;
518  size_t coords_ctr = 0;
519  for (const auto& rings_in_poly : poly_rings) {
520  OGRPolygon poly;
521  for (auto r = 0; r < rings_in_poly; r++) {
522  OGRLinearRing ring;
523  const auto ring_sz = ring_sizes[ring_ctr];
524  for (auto i = 0; i < 2 * ring_sz; i += 2) {
525  ring.addPoint(coords[coords_ctr + i], coords[coords_ctr + i + 1]);
526  }
527  ring.addPoint(coords[coords_ctr], coords[coords_ctr + 1]);
528  coords_ctr += 2 * ring_sz;
529  poly.addRing(&ring);
530  ring_ctr++;
531  }
532  multipoly->addGeometry(&poly);
533  }
534 }
535 
536 GeoMultiPolygon::GeoMultiPolygon(const std::string& wkt) {
537  const auto err = GeoBase::createFromWktString(wkt, &geom_);
538  if (err != OGRERR_NONE) {
539  throw GeoTypesError("MultiPolygon", err);
540  }
541  CHECK(geom_);
542  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbMultiPolygon) {
543  throw GeoTypesError("MultiPolygon",
544  "Unexpected geometry type from WKT string: " +
545  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
546  }
547 }
548 
549 void GeoMultiPolygon::getColumns(std::vector<double>& coords,
550  std::vector<int32_t>& ring_sizes,
551  std::vector<int32_t>& poly_rings,
552  std::vector<double>& bounds) const {
553  const auto mpoly = dynamic_cast<OGRMultiPolygon*>(geom_);
554  CHECK(mpoly);
555 
556  if (mpoly->IsEmpty()) {
557  // until the run-time can handle empties
558  throw GeoTypesError("MultiPolygon", "'EMPTY' not supported");
559  // return null bounds
560  bounds.push_back(NULL_DOUBLE);
561  bounds.push_back(NULL_DOUBLE);
562  bounds.push_back(NULL_DOUBLE);
563  bounds.push_back(NULL_DOUBLE);
564  return;
565  }
566 
567  BoundingBox bbox;
568  for (auto p = 0; p < mpoly->getNumGeometries(); p++) {
569  const auto mpoly_geom = mpoly->getGeometryRef(p);
570  CHECK(mpoly_geom);
571  const auto poly_geom = dynamic_cast<OGRPolygon*>(mpoly_geom);
572  if (!poly_geom) {
573  throw GeoTypesError("MultiPolygon",
574  "Failed to read polygon geometry from multipolygon");
575  }
576  const auto exterior_ring = poly_geom->getExteriorRing();
577  CHECK(exterior_ring);
578  // All exterior rings are imported CCW
579  if (exterior_ring->isClockwise()) {
580  exterior_ring->reverseWindingOrder();
581  }
582  const auto num_points_added = process_poly_ring(exterior_ring, coords, &bbox);
583  ring_sizes.push_back(num_points_added);
584 
585  for (auto r = 0; r < poly_geom->getNumInteriorRings(); r++) {
586  auto interior_ring = poly_geom->getInteriorRing(r);
587  CHECK(interior_ring);
588  // All interior rings are imported CW
589  if (!interior_ring->isClockwise()) {
590  interior_ring->reverseWindingOrder();
591  }
592  const auto num_points_added = process_poly_ring(interior_ring, coords, nullptr);
593  ring_sizes.push_back(num_points_added);
594  }
595  poly_rings.push_back(poly_geom->getNumInteriorRings() + 1);
596  }
597  bounds.push_back(bbox.min.x);
598  bounds.push_back(bbox.min.y);
599  bounds.push_back(bbox.max.x);
600  bounds.push_back(bbox.max.y);
601 }
602 
603 GeoGeometryCollection::GeoGeometryCollection(const std::string& wkt) {
604  const auto err = GeoBase::createFromWktString(wkt, &geom_);
605  if (err != OGRERR_NONE) {
606  throw GeoTypesError("GeometryCollection", err);
607  }
608  CHECK(geom_);
609  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbGeometryCollection) {
610  throw GeoTypesError("GeometryCollection",
611  "Unexpected geometry type from WKT string: " +
612  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
613  }
614 }
615 
616 namespace {
617 
619  uint8_t table_[128];
620  constexpr HexDigitToDecimalTable() : table_{} {
621  table_['1'] = 1;
622  table_['2'] = 2;
623  table_['3'] = 3;
624  table_['4'] = 4;
625  table_['5'] = 5;
626  table_['6'] = 6;
627  table_['7'] = 7;
628  table_['8'] = 8;
629  table_['9'] = 9;
630  table_['a'] = 10;
631  table_['A'] = 10;
632  table_['b'] = 11;
633  table_['B'] = 11;
634  table_['c'] = 12;
635  table_['C'] = 12;
636  table_['d'] = 13;
637  table_['D'] = 13;
638  table_['e'] = 14;
639  table_['E'] = 14;
640  table_['f'] = 15;
641  table_['F'] = 15;
642  }
643  constexpr uint8_t operator[](const char& hex_digit) const {
644  return (hex_digit < 0) ? 0 : table_[static_cast<int>(hex_digit)];
645  }
646 } constexpr hex_digit_to_decimal_table;
647 
648 inline uint8_t hex_to_binary(const char& usb, const char& lsb) {
649  return (hex_digit_to_decimal_table[usb] << 4) | hex_digit_to_decimal_table[lsb];
650 }
651 
652 std::vector<uint8_t> hex_string_to_binary_vector(const std::string& wkb_hex) {
653  auto num_bytes = wkb_hex.size() >> 1;
654  std::vector<uint8_t> wkb(num_bytes);
655  auto* chars = wkb_hex.data();
656  auto* bytes = wkb.data();
657  for (size_t i = 0; i < num_bytes; i++) {
658  auto const& usb = *chars++;
659  auto const& lsb = *chars++;
660  *bytes++ = hex_to_binary(usb, lsb);
661  }
662  return wkb;
663 }
664 
665 } // namespace
666 
667 OGRGeometry* GeoTypesFactory::createOGRGeometry(const std::string& wkt_or_wkb_hex) {
668  OGRGeometry* geom = nullptr;
669  OGRErr err = OGRERR_NONE;
670  if (wkt_or_wkb_hex.empty()) {
671  err = OGRERR_NOT_ENOUGH_DATA;
672  } else if (wkt_or_wkb_hex[0] == '0') { // all WKB hex strings start with a 0
673  err = GeoBase::createFromWkb(hex_string_to_binary_vector(wkt_or_wkb_hex), &geom);
674  } else {
675  err = GeoBase::createFromWktString(wkt_or_wkb_hex, &geom);
676  }
677  if (err != OGRERR_NONE) {
678  throw GeoTypesError("GeoFactory", err);
679  }
680  return geom;
681 }
682 
683 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoType(
684  const std::string& wkt_or_wkb_hex) {
685  return GeoTypesFactory::createGeoTypeImpl(createOGRGeometry(wkt_or_wkb_hex));
686 }
687 
688 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoType(const std::vector<uint8_t>& wkb) {
689  OGRGeometry* geom = nullptr;
690  const auto err = GeoBase::createFromWkb(wkb, &geom);
691  if (err != OGRERR_NONE) {
692  throw GeoTypesError("GeoFactory", err);
693  }
694  return GeoTypesFactory::createGeoTypeImpl(geom);
695 }
696 
697 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoType(OGRGeometry* geom) {
698  return GeoTypesFactory::createGeoTypeImpl(geom, false);
699 }
700 
701 bool GeoTypesFactory::getGeoColumns(const std::string& wkt_or_wkb_hex,
702  SQLTypeInfo& ti,
703  std::vector<double>& coords,
704  std::vector<double>& bounds,
705  std::vector<int>& ring_sizes,
706  std::vector<int>& poly_rings,
707  const bool promote_poly_to_mpoly) {
708  try {
709  if (wkt_or_wkb_hex.empty() || wkt_or_wkb_hex == "NULL") {
710  getNullGeoColumns(
711  ti, coords, bounds, ring_sizes, poly_rings, promote_poly_to_mpoly);
712  return true;
713  }
714 
715  const auto geospatial_base = GeoTypesFactory::createGeoType(wkt_or_wkb_hex);
716 
717  int srid = 0;
718  ti.set_input_srid(srid);
719  ti.set_output_srid(srid);
720 
721  getGeoColumnsImpl(geospatial_base,
722  ti,
723  coords,
724  bounds,
725  ring_sizes,
726  poly_rings,
727  promote_poly_to_mpoly);
728 
729  } catch (const std::exception& e) {
730  LOG(ERROR) << "Geospatial Import Error: " << e.what();
731  return false;
732  }
733 
734  return true;
735 }
736 
737 bool GeoTypesFactory::getGeoColumns(const std::vector<uint8_t>& wkb,
738  SQLTypeInfo& ti,
739  std::vector<double>& coords,
740  std::vector<double>& bounds,
741  std::vector<int>& ring_sizes,
742  std::vector<int>& poly_rings,
743  const bool promote_poly_to_mpoly) {
744  try {
745  const auto geospatial_base = GeoTypesFactory::createGeoType(wkb);
746 
747  int srid = 0;
748  ti.set_input_srid(srid);
749  ti.set_output_srid(srid);
750 
751  getGeoColumnsImpl(geospatial_base,
752  ti,
753  coords,
754  bounds,
755  ring_sizes,
756  poly_rings,
757  promote_poly_to_mpoly);
758 
759  } catch (const std::exception& e) {
760  LOG(ERROR) << "Geospatial Import Error: " << e.what();
761  return false;
762  }
763 
764  return true;
765 }
766 
767 bool GeoTypesFactory::getGeoColumns(OGRGeometry* geom,
768  SQLTypeInfo& ti,
769  std::vector<double>& coords,
770  std::vector<double>& bounds,
771  std::vector<int>& ring_sizes,
772  std::vector<int>& poly_rings,
773  const bool promote_poly_to_mpoly) {
774  try {
775  const auto geospatial_base = GeoTypesFactory::createGeoType(geom);
776 
777  int srid = 0;
778  ti.set_input_srid(srid);
779  ti.set_output_srid(srid);
780 
781  getGeoColumnsImpl(geospatial_base,
782  ti,
783  coords,
784  bounds,
785  ring_sizes,
786  poly_rings,
787  promote_poly_to_mpoly);
788 
789  } catch (const std::exception& e) {
790  LOG(ERROR) << "Geospatial Import Error: " << e.what();
791  return false;
792  }
793 
794  return true;
795 }
796 
797 bool GeoTypesFactory::getGeoColumns(const std::vector<std::string>* wkt_or_wkb_hex_column,
798  SQLTypeInfo& ti,
799  std::vector<std::vector<double>>& coords_column,
800  std::vector<std::vector<double>>& bounds_column,
801  std::vector<std::vector<int>>& ring_sizes_column,
802  std::vector<std::vector<int>>& poly_rings_column,
803  const bool promote_poly_to_mpoly) {
804  try {
805  for (const auto& wkt_or_wkb_hex : *wkt_or_wkb_hex_column) {
806  std::vector<double> coords;
807  std::vector<double> bounds;
808  std::vector<int> ring_sizes;
809  std::vector<int> poly_rings;
810 
811  SQLTypeInfo row_ti;
812  getGeoColumns(wkt_or_wkb_hex,
813  row_ti,
814  coords,
815  bounds,
816  ring_sizes,
817  poly_rings,
818  promote_poly_to_mpoly);
819 
820  if (ti.get_type() != row_ti.get_type()) {
821  throw GeoTypesError("GeoFactory", "Columnar: Geometry type mismatch");
822  }
823  coords_column.push_back(coords);
824  bounds_column.push_back(bounds);
825  ring_sizes_column.push_back(ring_sizes);
826  poly_rings_column.push_back(poly_rings);
827  }
828 
829  } catch (const std::exception& e) {
830  LOG(ERROR) << "Geospatial column Import Error: " << e.what();
831  return false;
832  }
833 
834  return true;
835 }
836 
837 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoTypeImpl(OGRGeometry* geom,
838  const bool owns_geom_obj) {
839  switch (wkbFlatten(geom->getGeometryType())) {
840  case wkbPoint:
841  return std::unique_ptr<GeoPoint>(new GeoPoint(geom, owns_geom_obj));
842  case wkbLineString:
843  return std::unique_ptr<GeoLineString>(new GeoLineString(geom, owns_geom_obj));
844  case wkbPolygon:
845  return std::unique_ptr<GeoPolygon>(new GeoPolygon(geom, owns_geom_obj));
846  case wkbMultiPolygon:
847  return std::unique_ptr<GeoMultiPolygon>(new GeoMultiPolygon(geom, owns_geom_obj));
848  case wkbGeometryCollection:
849  return std::unique_ptr<GeoGeometryCollection>(
850  new GeoGeometryCollection(geom, owns_geom_obj));
851  default:
852  throw GeoTypesError(
853  "GeoTypesFactory",
854  "Unrecognized geometry type: " + std::string(geom->getGeometryName()));
855  }
856 }
857 
858 void GeoTypesFactory::getGeoColumnsImpl(const std::unique_ptr<GeoBase>& geospatial_base,
859  SQLTypeInfo& ti,
860  std::vector<double>& coords,
861  std::vector<double>& bounds,
862  std::vector<int>& ring_sizes,
863  std::vector<int>& poly_rings,
864  const bool promote_poly_to_mpoly) {
865  switch (geospatial_base->getType()) {
867  const auto geospatial_point = dynamic_cast<GeoPoint*>(geospatial_base.get());
868  CHECK(geospatial_point);
869  geospatial_point->getColumns(coords);
870  ti.set_type(kPOINT);
871  break;
872  }
874  const auto geospatial_linestring =
875  dynamic_cast<GeoLineString*>(geospatial_base.get());
876  CHECK(geospatial_linestring);
877  geospatial_linestring->getColumns(coords, bounds);
878  ti.set_type(kLINESTRING);
879  break;
880  }
882  const auto geospatial_poly = dynamic_cast<GeoPolygon*>(geospatial_base.get());
883  CHECK(geospatial_poly);
884  geospatial_poly->getColumns(coords, ring_sizes, bounds);
885  if (promote_poly_to_mpoly) {
886  if (ring_sizes.size()) {
887  CHECK_GT(coords.size(), 0u);
888  poly_rings.push_back(1 + geospatial_poly->getNumInteriorRings());
889  }
890  }
891  ti.set_type(kPOLYGON);
892  break;
893  }
895  const auto geospatial_mpoly = dynamic_cast<GeoMultiPolygon*>(geospatial_base.get());
896  CHECK(geospatial_mpoly);
897  geospatial_mpoly->getColumns(coords, ring_sizes, poly_rings, bounds);
899  break;
900  }
902  case GeoBase::GeoType::kGEOMETRYCOLLECTION:
903  default:
904  throw std::runtime_error("Unrecognized geospatial type");
905  }
906 }
907 
908 void GeoTypesFactory::getNullGeoColumns(SQLTypeInfo& ti,
909  std::vector<double>& coords,
910  std::vector<double>& bounds,
911  std::vector<int>& ring_sizes,
912  std::vector<int>& poly_rings,
913  const bool promote_poly_to_mpoly) {
914  auto t = ti.get_type();
915  switch (t) {
916  case kPOINT: {
917  // NULL fixlen coords array
918  coords.push_back(NULL_ARRAY_DOUBLE);
919  coords.push_back(NULL_DOUBLE);
920  } break;
921  case kLINESTRING:
922  case kPOLYGON:
923  case kMULTIPOLYGON: {
924  // Leaving coords array empty
925  // NULL fixlen bounds array
926  bounds.push_back(NULL_ARRAY_DOUBLE);
927  bounds.push_back(NULL_DOUBLE);
928  bounds.push_back(NULL_DOUBLE);
929  bounds.push_back(NULL_DOUBLE);
930  // Leaving ring_sizes and poly_rings arrays empty
931  } break;
932  default:
933  throw std::runtime_error("Unsupported NULL geo");
934  }
935 }
936 
937 } // namespace Geospatial
OGRGeometry * geom_
Definition: Types.h:83
constexpr uint8_t operator[](const char &hex_digit) const
Definition: Types.cpp:643
#define NULL_DOUBLE
Definition: sqltypes.h:186
constexpr auto DOUBLE_MAX
Definition: Types.cpp:38
#define NULL_ARRAY_DOUBLE
Definition: sqltypes.h:194
#define LOG(tag)
Definition: Logger.h:188
Constants for Builtin SQL Types supported by OmniSci.
struct Geospatial::anonymous_namespace{Types.cpp}::HexDigitToDecimalTable hex_digit_to_decimal_table
#define CHECK_GT(x, y)
Definition: Logger.h:209
std::string to_string(char const *&&v)
void set_input_srid(int d)
Definition: sqltypes.h:353
int process_poly_ring(OGRLinearRing *ring, std::vector< double > &coords, BoundingBox *bbox)
Definition: Types.cpp:68
void update(double x, double y)
Definition: Types.cpp:52
void set_output_srid(int s)
Definition: sqltypes.h:355
bool operator==(const SlotSize &lhs, const SlotSize &rhs)
constexpr auto DOUBLE_MIN
Definition: Types.cpp:39
#define CHECK(condition)
Definition: Logger.h:197
HOST DEVICE SQLTypes get_type() const
Definition: sqltypes.h:259
static bool run
uint8_t hex_to_binary(const char &usb, const char &lsb)
Definition: Types.cpp:648
std::vector< uint8_t > hex_string_to_binary_vector(const std::string &wkb_hex)
Definition: Types.cpp:652
HOST DEVICE void set_type(SQLTypes t)
Definition: sqltypes.h:349