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