OmniSciDB  ca0c39ec8f
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Types.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2022 HEAVY.AI, 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 #include <mutex>
21 
22 #include <gdal.h>
23 #include <ogr_geometry.h>
24 #include <ogrsf_frmts.h>
25 
26 #include "Geospatial/GDAL.h"
27 #include "Logger/Logger.h"
28 #include "Shared/sqltypes.h"
29 
39 namespace {
40 constexpr auto DOUBLE_MAX = std::numeric_limits<double>::max();
41 constexpr auto DOUBLE_MIN = std::numeric_limits<double>::lowest();
42 
43 struct Coords {
44  double x;
45  double y;
46  Coords(double x, double y) : x(x), y(y) {}
47 };
48 
49 struct BoundingBox {
52  BoundingBox() : min(DOUBLE_MAX, DOUBLE_MAX), max(DOUBLE_MIN, DOUBLE_MIN) {}
53 
54  void update(double x, double y) {
55  if (x < min.x) {
56  min.x = x;
57  }
58  if (y < min.y) {
59  min.y = y;
60  }
61  if (x > max.x) {
62  max.x = x;
63  }
64  if (y > max.y) {
65  max.y = y;
66  }
67  }
68 };
69 
70 int process_poly_ring(OGRLinearRing* ring,
71  std::vector<double>& coords,
72  BoundingBox* bbox) {
73  double last_x = DOUBLE_MAX, last_y = DOUBLE_MAX;
74  size_t first_index = coords.size();
75  int num_points_added = 0;
76  int num_points_in_ring = ring->getNumPoints();
77  if (num_points_in_ring < 3) {
79  "PolyRing",
80  "All poly rings must have more than 3 points. Found ring with " +
81  std::to_string(num_points_in_ring) + " points.");
82  }
83  for (auto i = 0; i < num_points_in_ring; i++) {
84  OGRPoint point;
85  ring->getPoint(i, &point);
86  last_x = point.getX();
87  last_y = point.getY();
88  coords.push_back(last_x);
89  coords.push_back(last_y);
90  if (bbox) {
91  bbox->update(last_x, last_y);
92  }
93  num_points_added++;
94  }
95  // Store all rings as open rings (implicitly assumes all rings are closed)
96  if ((coords[first_index] == last_x) && (coords[first_index + 1] == last_y)) {
97  coords.pop_back();
98  coords.pop_back();
99  num_points_added--;
100  if (num_points_added < 3) {
102  "PolyRing",
103  "All exterior rings must have more than 3 points. Found ring with " +
104  std::to_string(num_points_added) + " points.");
105  }
106  }
107  return num_points_added;
108 }
109 
110 } // namespace
111 
112 namespace Geospatial {
113 
115 std::map<std::tuple<int32_t, int32_t>, std::shared_ptr<OGRCoordinateTransformation>>
117 
118 std::string GeoTypesError::OGRErrorToStr(const int ogr_err) {
119  switch (ogr_err) {
120  case OGRERR_NOT_ENOUGH_DATA:
121  return std::string("not enough input data");
122  case OGRERR_NOT_ENOUGH_MEMORY:
123  return std::string("not enough memory");
124  case OGRERR_UNSUPPORTED_GEOMETRY_TYPE:
125  return std::string("unsupported geometry type");
126  case OGRERR_UNSUPPORTED_OPERATION:
127  return std::string("unsupported operation");
128  case OGRERR_CORRUPT_DATA:
129  return std::string("corrupt input data");
130  case OGRERR_FAILURE:
131  return std::string("ogr failure");
132  case OGRERR_UNSUPPORTED_SRS:
133  return std::string("unsupported spatial reference system");
134  case OGRERR_INVALID_HANDLE:
135  return std::string("invalid file handle");
136 #if (GDAL_VERSION_MAJOR > 1)
137  case OGRERR_NON_EXISTING_FEATURE:
138  return std::string("feature does not exist in input geometry");
139 #endif
140  default:
141  return std::string("Unknown OGOR error encountered: ") + std::to_string(ogr_err);
142  }
143 }
144 
146  // Note: Removing the geometry object that was pulled from an OGRFeature results in a
147  // segfault. If we are wrapping around a pre-existing OGRGeometry object, we let the
148  // caller manage the memory.
149  if (geom_ && owns_geom_obj_) {
150  OGRGeometryFactory::destroyGeometry(geom_);
151  }
152 }
153 
154 OGRErr GeoBase::createFromWktString(const std::string& wkt, OGRGeometry** geom) {
155 #if (GDAL_VERSION_MAJOR > 2) || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 3)
156  OGRErr ogr_status = OGRGeometryFactory::createFromWkt(wkt.c_str(), nullptr, geom);
157 #else
158  auto data = (char*)wkt.c_str();
159  OGRErr ogr_status = OGRGeometryFactory::createFromWkt(&data, NULL, geom);
160 #endif
161  return ogr_status;
162 }
163 
164 std::string GeoBase::getWktString() const {
165  char* wkt = nullptr;
166  geom_->exportToWkt(&wkt);
167  CHECK(wkt);
168  std::string wkt_str(wkt);
169  CPLFree(wkt);
170  return wkt_str;
171 }
172 
173 OGRErr GeoBase::createFromWkbView(OGRGeometry** geom, WkbView const wkb_view) {
174 #if (GDAL_VERSION_MAJOR > 2) || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 3)
175  OGRErr ogr_status =
176  OGRGeometryFactory::createFromWkb(wkb_view.ptr_, nullptr, geom, wkb_view.size_);
177  return ogr_status;
178 #else
179  CHECK(false);
180 #endif
181 }
182 
183 // GeoBase could also generate GEOS geometries through geom_->exportToGEOS(),
184 // conversion to WKB and subsequent load of WKB into GEOS could be eliminated.
185 
186 bool GeoBase::getWkb(std::vector<uint8_t>& wkb) const {
187  auto size = geom_->WkbSize();
188  if (size > 0) {
189  wkb.resize(size);
190  geom_->exportToWkb(wkbNDR, wkb.data());
191  return true;
192  }
193  return false;
194 }
195 
196 bool GeoBase::isEmpty() const {
197  return geom_ && geom_->IsEmpty();
198 }
199 
200 bool GeoBase::operator==(const GeoBase& other) const {
201  if (!this->geom_ || !other.geom_) {
202  return false;
203  }
204  return this->geom_->Equals(other.geom_);
205 }
206 
208 #define SRID_WORLD_MERCATOR 999000
209 
210 #define SRID_NORTH_UTM_START 999001
211 
212 #define SRID_NORTH_UTM_END 999060
213 
214 #define SRID_NORTH_LAMBERT 999061
215 
216 #define SRID_SOUTH_UTM_START 999101
217 
218 #define SRID_SOUTH_UTM_END 999160
219 
220 #define SRID_SOUTH_LAMBERT 999161
221 
222 #define SRID_LAEA_START 999163
223 
224 #define SRID_LAEA_END 999283
225 
226 int32_t GeoBase::getBestPlanarSRID() const {
227  if (!this->geom_) {
228  return 0;
229  }
230  double cx, cy, xwidth, ywidth;
231  OGREnvelope envelope;
232  geom_->getEnvelope(&envelope);
233  // Can't use GDAL's Centroid geom_->Centroid(OGRPoint*): requires GEOS
234  // Use center of the bounding box for now.
235  // TODO: hook up our own Centroid implementation
236  cx = (envelope.MaxX + envelope.MinX) / 2.0;
237  cy = (envelope.MaxY + envelope.MinY) / 2.0;
238  xwidth = envelope.MaxX - envelope.MinX;
239  ywidth = envelope.MaxY - envelope.MinY;
240 
241  // Arctic coords: Lambert Azimuthal Equal Area North
242  if (cy > 70.0 && ywidth < 45.0) {
243  return SRID_NORTH_LAMBERT;
244  }
245  // Antarctic coords: Lambert Azimuthal Equal Area South
246  if (cy < -70.0 && ywidth < 45.0) {
247  return SRID_SOUTH_LAMBERT;
248  }
249 
250  // Can geometry fit into a single UTM zone?
251  if (xwidth < 6.0) {
252  int zone = floor((cx + 180.0) / 6.0);
253  if (zone > 59) {
254  zone = 59;
255  }
256  // Below the equator: UTM South
257  // Above the equator: UTM North
258  if (cy < 0.0) {
259  return SRID_SOUTH_UTM_START + zone;
260  } else {
261  return SRID_NORTH_UTM_START + zone;
262  }
263  }
264 
265  // TODO: to be removed once we add custom LAEA zone transforms
266  // Can geometry still fit into 5 consecutive UTM zones?
267  // Then go for the mid-UTM zone, tolerating some limited distortion
268  // in the left and right corners. That's still better than Mercator.
269  if (xwidth < 30.0) {
270  int zone = floor((cx + 180.0) / 6.0);
271  if (zone > 59) {
272  zone = 59;
273  }
274  // Below the equator: UTM South
275  // Above the equator: UTM North
276  if (cy < 0.0) {
277  return SRID_SOUTH_UTM_START + zone;
278  } else {
279  return SRID_NORTH_UTM_START + zone;
280  }
281  }
282 
283  // Can geometry fit into a custom LAEA area 30 degrees high? Allow some overlap.
284  if (ywidth < 25.0) {
285  int xzone = -1;
286  int yzone = 3 + floor(cy / 30.0); // range of 0-5
287  if ((yzone == 2 || yzone == 3) && xwidth < 30.0) {
288  // Equatorial band, 12 zones, 30 degrees wide
289  xzone = 6 + floor(cx / 30.0);
290  } else if ((yzone == 1 || yzone == 4) && xwidth < 45.0) {
291  // Temperate band, 8 zones, 45 degrees wide
292  xzone = 4 + floor(cx / 45.0);
293  } else if ((yzone == 0 || yzone == 5) && xwidth < 90.0) {
294  // Arctic band, 4 zones, 90 degrees wide
295  xzone = 2 + floor(cx / 90.0);
296  }
297  // Found an appropriate xzone to fit in?
298  if (xzone != -1) {
299  return SRID_LAEA_START + 20 * yzone + xzone;
300  }
301  }
302 
303  // Fall-back to Mercator
304  return SRID_WORLD_MERCATOR;
305 }
306 
307 std::shared_ptr<OGRCoordinateTransformation> GeoBase::getTransformation(int32_t srid0,
308  int32_t srid1) {
309  std::lock_guard<std::mutex> guard(transformation_map_mutex_);
310  std::tuple<int32_t, int32_t> key{srid0, srid1};
311  auto it = transformation_map_.find(key);
312  if (it != transformation_map_.end()) {
313  return it->second;
314  }
315  auto setSpatialReference = [&](OGRSpatialReference* sr, int32_t srid) -> bool {
316  OGRErr status = OGRERR_NONE;
317  if (srid == 4326) {
318  status = sr->importFromEPSG(4326);
319  } else if (srid == SRID_NORTH_LAMBERT) {
320  // +proj=laea +lat_0=90 +lon_0=-40 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m
321  // +no_defs
322  status = sr->importFromEPSG(3574);
323  } else if (srid == SRID_SOUTH_LAMBERT) {
324  // +proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m
325  // +no_defs
326  status = sr->importFromEPSG(3409);
327  } else if (SRID_SOUTH_UTM_START <= srid && srid <= SRID_SOUTH_UTM_END) {
328  // +proj=utm +zone=%d +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs
329  int32_t zone = srid - SRID_SOUTH_UTM_START;
330  status = sr->importFromEPSG(32701 + zone);
331  } else if (SRID_NORTH_UTM_START <= srid && srid <= SRID_NORTH_UTM_END) {
332  // +proj=utm +zone=%d +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
333  int32_t zone = srid - SRID_NORTH_UTM_START;
334  status = sr->importFromEPSG(32601 + zone);
335  } else if (SRID_LAEA_START <= srid && srid <= SRID_LAEA_END) {
336  // TODO: add support and coordinate operations for custom Lambert zones,
337  // need to calculate lat/lon for the zone, SetCoordinateOperation in options.
338  // +proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=%g +lon_0=%g +units=m +no_defs
339  // Go with Mercator for now
340  status = sr->importFromEPSG(3395);
341  } else if (srid == SRID_WORLD_MERCATOR) {
342  // +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m
343  // +no_defs
344  status = sr->importFromEPSG(3395);
345  } else if (srid > 0) {
346  // Attempt to import from srid directly
347  status = sr->importFromEPSG(srid);
348  } else {
349  return false;
350  }
351 #if GDAL_VERSION_MAJOR >= 3
352  // GDAL 3.x (really Proj.4 6.x) now enforces lat, lon order
353  // this results in X and Y being transposed for angle-based
354  // coordinate systems. This restores the previous behavior.
355  if (status == OGRERR_NONE) {
356  sr->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
357  }
358 #endif
359  return (status == OGRERR_NONE);
360  };
361 
362  // lazy init GDAL
363  GDAL::init();
364 
365  OGRSpatialReference sr0;
366  if (!setSpatialReference(&sr0, srid0)) {
367  return nullptr;
368  }
369  OGRSpatialReference sr1;
370  if (!setSpatialReference(&sr1, srid1)) {
371  return nullptr;
372  }
373  // GDAL 3 allows specification of advanced transformations in
374  // OGRCoordinateTransformationOptions, including multi-step pipelines.
375  // GDAL 3 would be required to handle Lambert zone proj4 strings.
376  // Using a simple transform for now.
377  std::shared_ptr<OGRCoordinateTransformation> new_transformation;
378  new_transformation.reset(OGRCreateCoordinateTransformation(&sr0, &sr1));
379  transformation_map_[key] = new_transformation;
380  return new_transformation;
381 }
382 
383 bool GeoBase::transform(int32_t srid0, int32_t srid1) {
384  auto coordinate_transformation = getTransformation(srid0, srid1);
385  if (!coordinate_transformation) {
386  return false;
387  }
388  auto ogr_status = geom_->transform(coordinate_transformation.get());
389  return (ogr_status == OGRERR_NONE);
390 }
391 
393  auto srid1 = ti.get_output_srid();
394  if (srid1 == 4326) {
395  auto srid0 = ti.get_input_srid();
396  if (srid0 > 0 && srid0 != 4326) {
397  if (!transform(srid0, srid1)) {
398  return false;
399  }
400  }
401  }
402  return true;
403 }
404 
405 // Run a specific geo operation on this and other geometries
406 std::unique_ptr<GeoBase> GeoBase::run(GeoBase::GeoOp op, const GeoBase& other) const {
407  OGRGeometry* result = nullptr;
408  // Invalid geometries are derailing geo op performance,
409  // Checking validity before running an operation doesn't have lesser penalty.
410  // DOES NOT HELP: if (geom_->IsValid() && other.geom_->IsValid()) {
411  switch (op) {
413  result = geom_->Intersection(other.geom_);
414  break;
416  result = geom_->Difference(other.geom_);
417  break;
419  result = geom_->Union(other.geom_);
420  break;
421  default:
422  break;
423  }
424  // TODO: Need to handle empty/non-POLYGON result
425  if (!result || result->IsEmpty() ||
426  !(result->getGeometryType() == wkbPolygon ||
427  result->getGeometryType() == wkbMultiPolygon)) {
428  // throw GeoTypesError(std::string(OGRGeometryTypeToName(geom_->getGeometryType())),
429  // "Currenly don't support invalid or empty result");
430  // return GeoTypesFactory::createGeoType("POLYGON EMPTY");
431  // Not supporting EMPTY polygons, return a dot polygon
433  "MULTIPOLYGON(((0 0,0.0000001 0,0 0.0000001)))");
434  }
435  return GeoTypesFactory::createGeoType(result);
436 }
437 
438 // Run a specific geo operation on this and other geometries
439 std::unique_ptr<GeoBase> GeoBase::optimized_run(GeoBase::GeoOp op,
440  const GeoBase& other) const {
441  OGRGeometry* result = nullptr;
442  // Loop through polys combinations, check validity, do intersections
443  // where needed, return a union of all intersections
444  auto gc1 = geom_->toGeometryCollection();
445  auto gc2 = other.geom_->toGeometryCollection();
446  if (!gc1 || !gc2 || gc1->IsEmpty() || gc2->IsEmpty()) {
447  return nullptr;
448  }
449  for (int i1 = 0; i1 < gc1->getNumGeometries(); i1++) {
450  auto g1 = gc1->getGeometryRef(i1);
451  // Validity check is very slow
452  if (!g1 || g1->IsEmpty() /*|| !g1->IsValid()*/) {
453  continue;
454  }
455  OGREnvelope ge1;
456  g1->getEnvelope(&ge1);
457  for (int i2 = 0; i2 < gc2->getNumGeometries(); i2++) {
458  auto g2 = gc2->getGeometryRef(i2);
459  // Validity check is very slow
460  if (!g2 || g2->IsEmpty() /*|| !g2->IsValid()*/) {
461  continue;
462  }
463  // Check for bounding box overlap
464  OGREnvelope ge2;
465  g2->getEnvelope(&ge2);
466  if (!ge1.Intersects(ge2)) {
467  continue;
468  }
469  // Do intersection
470  auto intermediate_result = g1->Intersection(g2);
471  if (!intermediate_result || intermediate_result->IsEmpty()) {
472  continue;
473  }
474  if (!result) {
475  result = intermediate_result;
476  } else {
477  result = result->Union(intermediate_result);
478  }
479  }
480  }
481 
482  // TODO: Need to handle empty/non-POLYGON result
483  if (!result || result->IsEmpty() ||
484  !(result->getGeometryType() == wkbPolygon ||
485  result->getGeometryType() == wkbMultiPolygon)) {
486  // throw GeoTypesError(std::string(OGRGeometryTypeToName(geom_->getGeometryType())),
487  // "Currenly don't support invalid or empty result");
488  // return GeoTypesFactory::createGeoType("POLYGON EMPTY");
489  // Not supporting EMPTY polygons, return a dot polygon
491  "MULTIPOLYGON(((0 0,0.0000001 0,0 0.0000001)))");
492  }
493  return GeoTypesFactory::createGeoType(result);
494 }
495 
496 // Run a specific geo operation on this geometry and a double param
497 std::unique_ptr<GeoBase> GeoBase::run(GeoBase::GeoOp op, double param) const {
498  OGRGeometry* result = nullptr;
499  switch (op) {
501  result = geom_->Buffer(param);
502  break;
503  default:
504  break;
505  }
506  // TODO: Need to handle empty/non-POLYGON result
507  if (!result || result->IsEmpty() ||
508  !(result->getGeometryType() == wkbPolygon ||
509  result->getGeometryType() == wkbMultiPolygon)) {
510  // throw GeoTypesError(std::string(OGRGeometryTypeToName(geom_->getGeometryType())),
511  // "Currenly don't support invalid or empty result");
512  // return GeoTypesFactory::createGeoType("POLYGON EMPTY");
513  // Not supporting EMPTY polygons, return a dot polygon
515  "MULTIPOLYGON(((0 0,0.0000001 0,0 0.0000001)))");
516  }
517  return GeoTypesFactory::createGeoType(result);
518 }
519 
520 // Run a specific predicate operation on this geometry
521 bool GeoBase::run(GeoBase::GeoOp op) const {
522  auto result = false;
523  switch (op) {
525  result = geom_->IsValid();
526  break;
528  result = isEmpty();
529  break;
530  default:
531  break;
532  }
533  return result;
534 }
535 
536 std::unique_ptr<GeoBase> GeoPoint::clone() const {
537  CHECK(geom_);
538  return std::unique_ptr<GeoBase>(new GeoPoint(geom_->clone(), true));
539 }
540 
541 GeoPoint::GeoPoint(const std::vector<double>& coords) {
542  if (coords.size() != 2) {
543  throw GeoTypesError("Point",
544  "Incorrect coord size of " + std::to_string(coords.size()) +
545  " supplied. Expected 2.");
546  }
547  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbPoint);
548  OGRPoint* point = dynamic_cast<OGRPoint*>(geom_);
549  CHECK(point);
550  point->setX(coords[0]);
551  point->setY(coords[1]);
552 }
553 
554 GeoPoint::GeoPoint(const std::string& wkt) {
555  const auto err = GeoBase::createFromWktString(wkt, &geom_);
556  if (err != OGRERR_NONE) {
557  throw GeoTypesError("Point", err);
558  }
559  CHECK(geom_);
560  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbPoint) {
561  throw GeoTypesError("Point",
562  "Unexpected geometry type from WKT string: " +
563  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
564  }
565 }
566 
567 void GeoPoint::getColumns(std::vector<double>& coords) const {
568  const auto point_geom = dynamic_cast<OGRPoint*>(geom_);
569  CHECK(point_geom);
570 
571  if (point_geom->IsEmpty()) {
572  // until the run-time can handle empties
573  throw GeoTypesError("Point", "'EMPTY' not supported");
574  // we cannot yet handle NULL fixed-length array
575  // so we have to store sentinel values instead
576  coords.push_back(NULL_DOUBLE);
577  coords.push_back(NULL_DOUBLE);
578  return;
579  }
580 
581  coords.push_back(point_geom->getX());
582  coords.push_back(point_geom->getY());
583 }
584 
585 std::unique_ptr<GeoBase> GeoMultiPoint::clone() const {
586  CHECK(geom_);
587  return std::unique_ptr<GeoBase>(new GeoMultiPoint(geom_->clone(), true));
588 }
589 
590 GeoMultiPoint::GeoMultiPoint(const std::vector<double>& coords) {
591  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbMultiPoint);
592  OGRMultiPoint* multipoint = dynamic_cast<OGRMultiPoint*>(geom_);
593  CHECK(multipoint);
594  for (size_t i = 0; i < coords.size(); i += 2) {
595  OGRPoint pt(coords[i], coords[i + 1]);
596  multipoint->addGeometry(&pt);
597  }
598 }
599 
600 GeoMultiPoint::GeoMultiPoint(const std::string& wkt) {
601  const auto err = GeoBase::createFromWktString(wkt, &geom_);
602  if (err != OGRERR_NONE) {
603  throw GeoTypesError("MultiPoint", err);
604  }
605  CHECK(geom_);
606  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbMultiPoint) {
607  throw GeoTypesError("MultiPoint",
608  "Unexpected geometry type from WKT string: " +
609  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
610  }
611 }
612 
613 void GeoMultiPoint::getColumns(std::vector<double>& coords,
614  std::vector<double>& bounds) const {
615  auto multipoint_geom = dynamic_cast<OGRMultiPoint*>(geom_);
616  CHECK(multipoint_geom);
617 
618  if (multipoint_geom->IsEmpty()) {
619  // until the run-time can handle empties
620  throw GeoTypesError("MultiPoint", "'EMPTY' not supported");
621  // return null bounds
622  bounds.push_back(NULL_DOUBLE);
623  bounds.push_back(NULL_DOUBLE);
624  bounds.push_back(NULL_DOUBLE);
625  bounds.push_back(NULL_DOUBLE);
626  return;
627  }
628 
629  BoundingBox bbox;
630  for (auto i = 0; i < multipoint_geom->getNumGeometries(); i++) {
631  OGRPoint* point = dynamic_cast<OGRPoint*>(multipoint_geom->getGeometryRef(i));
632  CHECK(point);
633  double x = point->getX();
634  double y = point->getY();
635  coords.push_back(x);
636  coords.push_back(y);
637  bbox.update(x, y);
638  }
639  bounds.push_back(bbox.min.x);
640  bounds.push_back(bbox.min.y);
641  bounds.push_back(bbox.max.x);
642  bounds.push_back(bbox.max.y);
643 }
644 
645 std::unique_ptr<GeoBase> GeoLineString::clone() const {
646  CHECK(geom_);
647  return std::unique_ptr<GeoBase>(new GeoLineString(geom_->clone(), true));
648 }
649 
650 GeoLineString::GeoLineString(const std::vector<double>& coords) {
651  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbLineString);
652  OGRLineString* line = dynamic_cast<OGRLineString*>(geom_);
653  CHECK(line);
654  for (size_t i = 0; i < coords.size(); i += 2) {
655  line->addPoint(coords[i], coords[i + 1]);
656  }
657 }
658 
659 GeoLineString::GeoLineString(const std::string& wkt) {
660  const auto err = GeoBase::createFromWktString(wkt, &geom_);
661  if (err != OGRERR_NONE) {
662  throw GeoTypesError("LineString", err);
663  }
664  CHECK(geom_);
665  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbLineString) {
666  throw GeoTypesError("LineString",
667  "Unexpected geometry type from WKT string: " +
668  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
669  }
670 }
671 
672 void GeoLineString::getColumns(std::vector<double>& coords,
673  std::vector<double>& bounds) const {
674  auto linestring_geom = dynamic_cast<OGRLineString*>(geom_);
675  CHECK(linestring_geom);
676 
677  if (linestring_geom->IsEmpty()) {
678  // until the run-time can handle empties
679  throw GeoTypesError("LineString", "'EMPTY' not supported");
680  // return null bounds
681  bounds.push_back(NULL_DOUBLE);
682  bounds.push_back(NULL_DOUBLE);
683  bounds.push_back(NULL_DOUBLE);
684  bounds.push_back(NULL_DOUBLE);
685  return;
686  }
687 
688  BoundingBox bbox;
689  for (auto i = 0; i < linestring_geom->getNumPoints(); i++) {
690  OGRPoint point;
691  linestring_geom->getPoint(i, &point);
692  double x = point.getX();
693  double y = point.getY();
694  coords.push_back(x);
695  coords.push_back(y);
696  bbox.update(x, y);
697  }
698  bounds.push_back(bbox.min.x);
699  bounds.push_back(bbox.min.y);
700  bounds.push_back(bbox.max.x);
701  bounds.push_back(bbox.max.y);
702 }
703 
704 std::unique_ptr<GeoBase> GeoMultiLineString::clone() const {
705  CHECK(geom_);
706  return std::unique_ptr<GeoBase>(new GeoMultiLineString(geom_->clone(), true));
707 }
708 
709 GeoMultiLineString::GeoMultiLineString(const std::vector<double>& coords,
710  const std::vector<int32_t>& linestring_sizes) {
711  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbMultiLineString);
712  OGRMultiLineString* multilinestring = dynamic_cast<OGRMultiLineString*>(geom_);
713  CHECK(multilinestring);
714 
715  size_t coords_ctr = 0;
716  for (const auto linestring_sz : linestring_sizes) {
717  OGRLineString linestring;
718  auto next_coords_ctr = coords_ctr + 2 * linestring_sz;
719  CHECK(next_coords_ctr <= coords.size());
720  for (auto i = coords_ctr; i < next_coords_ctr; i += 2) {
721  linestring.addPoint(coords[i], coords[i + 1]);
722  }
723  coords_ctr = next_coords_ctr;
724  multilinestring->addGeometry(&linestring);
725  }
726 }
727 
728 GeoMultiLineString::GeoMultiLineString(const std::string& wkt) {
729  const auto err = GeoBase::createFromWktString(wkt, &geom_);
730  if (err != OGRERR_NONE) {
731  throw GeoTypesError("MultiLineString", err);
732  }
733  CHECK(geom_);
734  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbMultiLineString) {
735  throw GeoTypesError("MultiLineString",
736  "Unexpected geometry type from WKT string: " +
737  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
738  }
739 }
740 
741 void GeoMultiLineString::getColumns(std::vector<double>& coords,
742  std::vector<int32_t>& linestring_sizes,
743  std::vector<double>& bounds) const {
744  auto multilinestring = dynamic_cast<OGRMultiLineString*>(geom_);
745  CHECK(multilinestring);
746 
747  if (multilinestring->IsEmpty()) {
748  // until the run-time can handle empties
749  throw GeoTypesError("MultiLineString", "'EMPTY' not supported");
750  // return null bounds
751  bounds.push_back(NULL_DOUBLE);
752  bounds.push_back(NULL_DOUBLE);
753  bounds.push_back(NULL_DOUBLE);
754  bounds.push_back(NULL_DOUBLE);
755  return;
756  }
757 
758  BoundingBox bbox;
759  for (auto l = 0; l < multilinestring->getNumGeometries(); l++) {
760  const auto geom = multilinestring->getGeometryRef(l);
761  CHECK(geom);
762  const auto linestring = dynamic_cast<OGRLineString*>(geom);
763  if (!linestring) {
764  throw GeoTypesError("MultiLineString",
765  "Failed to read linestring geometry from multilinestring");
766  }
767  auto linestring_sz = linestring->getNumPoints();
768  linestring_sizes.push_back(linestring_sz);
769  for (auto i = 0; i < linestring_sz; i++) {
770  OGRPoint point;
771  linestring->getPoint(i, &point);
772  double x = point.getX();
773  double y = point.getY();
774  coords.push_back(x);
775  coords.push_back(y);
776  bbox.update(x, y);
777  }
778  }
779  bounds.push_back(bbox.min.x);
780  bounds.push_back(bbox.min.y);
781  bounds.push_back(bbox.max.x);
782  bounds.push_back(bbox.max.y);
783 }
784 
785 std::unique_ptr<GeoBase> GeoPolygon::clone() const {
786  CHECK(geom_);
787  return std::unique_ptr<GeoBase>(new GeoPolygon(geom_->clone(), true));
788 }
789 
790 GeoPolygon::GeoPolygon(const std::vector<double>& coords,
791  const std::vector<int32_t>& ring_sizes) {
792  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbPolygon);
793  OGRPolygon* poly = dynamic_cast<OGRPolygon*>(geom_);
794  CHECK(poly);
795 
796  size_t coords_ctr = 0;
797  for (size_t r = 0; r < ring_sizes.size(); r++) {
798  OGRLinearRing ring;
799  const auto ring_sz = ring_sizes[r];
800  for (auto i = 0; i < 2 * ring_sz; i += 2) {
801  ring.addPoint(coords[coords_ctr + i], coords[coords_ctr + i + 1]);
802  }
803  ring.addPoint(coords[coords_ctr], coords[coords_ctr + 1]);
804  coords_ctr += 2 * ring_sz;
805  poly->addRing(&ring);
806  }
807 }
808 
809 GeoPolygon::GeoPolygon(const std::string& wkt) {
810  const auto err = GeoBase::createFromWktString(wkt, &geom_);
811  if (err != OGRERR_NONE) {
812  throw GeoTypesError("Polygon", err);
813  }
814  CHECK(geom_);
815  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbPolygon) {
816  throw GeoTypesError("Polygon",
817  "Unexpected geometry type from WKT string: " +
818  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
819  }
820 }
821 
822 void GeoPolygon::getColumns(std::vector<double>& coords,
823  std::vector<int32_t>& ring_sizes,
824  std::vector<double>& bounds) const {
825  const auto poly_geom = dynamic_cast<OGRPolygon*>(geom_);
826  CHECK(poly_geom);
827 
828  if (poly_geom->IsEmpty()) {
829  // until the run-time can handle empties
830  throw GeoTypesError("Polygon", "'EMPTY' not supported");
831  // return null bounds
832  bounds.push_back(NULL_DOUBLE);
833  bounds.push_back(NULL_DOUBLE);
834  bounds.push_back(NULL_DOUBLE);
835  bounds.push_back(NULL_DOUBLE);
836  return;
837  }
838 
839  BoundingBox bbox;
840  const auto exterior_ring = poly_geom->getExteriorRing();
841  CHECK(exterior_ring);
842  // All exterior rings are imported CCW
843  if (exterior_ring->isClockwise()) {
844  exterior_ring->reverseWindingOrder();
845  }
846  const auto num_points_added = process_poly_ring(exterior_ring, coords, &bbox);
847  ring_sizes.push_back(num_points_added);
848  for (auto r = 0; r < poly_geom->getNumInteriorRings(); r++) {
849  auto interior_ring = poly_geom->getInteriorRing(r);
850  CHECK(interior_ring);
851  // All interior rings are imported CW
852  if (!interior_ring->isClockwise()) {
853  interior_ring->reverseWindingOrder();
854  }
855  const auto num_points_added = process_poly_ring(interior_ring, coords, nullptr);
856  ring_sizes.push_back(num_points_added);
857  }
858  bounds.push_back(bbox.min.x);
859  bounds.push_back(bbox.min.y);
860  bounds.push_back(bbox.max.x);
861  bounds.push_back(bbox.max.y);
862 }
863 
865  const auto poly_geom = dynamic_cast<OGRPolygon*>(geom_);
866  CHECK(poly_geom);
867  return poly_geom->getNumInteriorRings();
868 }
869 
870 std::unique_ptr<GeoBase> GeoMultiPolygon::clone() const {
871  CHECK(geom_);
872  return std::unique_ptr<GeoBase>(new GeoMultiPolygon(geom_->clone(), true));
873 }
874 
875 GeoMultiPolygon::GeoMultiPolygon(const std::vector<double>& coords,
876  const std::vector<int32_t>& ring_sizes,
877  const std::vector<int32_t>& poly_rings) {
878  geom_ = OGRGeometryFactory::createGeometry(OGRwkbGeometryType::wkbMultiPolygon);
879  OGRMultiPolygon* multipoly = dynamic_cast<OGRMultiPolygon*>(geom_);
880  CHECK(multipoly);
881 
882  size_t ring_ctr = 0;
883  size_t coords_ctr = 0;
884  for (const auto& rings_in_poly : poly_rings) {
885  OGRPolygon poly;
886  for (auto r = 0; r < rings_in_poly; r++) {
887  OGRLinearRing ring;
888  const auto ring_sz = ring_sizes[ring_ctr];
889  for (auto i = 0; i < 2 * ring_sz; i += 2) {
890  ring.addPoint(coords[coords_ctr + i], coords[coords_ctr + i + 1]);
891  }
892  ring.addPoint(coords[coords_ctr], coords[coords_ctr + 1]);
893  coords_ctr += 2 * ring_sz;
894  poly.addRing(&ring);
895  ring_ctr++;
896  }
897  multipoly->addGeometry(&poly);
898  }
899 }
900 
901 GeoMultiPolygon::GeoMultiPolygon(const std::string& wkt) {
902  const auto err = GeoBase::createFromWktString(wkt, &geom_);
903  if (err != OGRERR_NONE) {
904  throw GeoTypesError("MultiPolygon", err);
905  }
906  CHECK(geom_);
907  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbMultiPolygon) {
908  throw GeoTypesError("MultiPolygon",
909  "Unexpected geometry type from WKT string: " +
910  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
911  }
912 }
913 
914 void GeoMultiPolygon::getColumns(std::vector<double>& coords,
915  std::vector<int32_t>& ring_sizes,
916  std::vector<int32_t>& poly_rings,
917  std::vector<double>& bounds) const {
918  const auto mpoly = dynamic_cast<OGRMultiPolygon*>(geom_);
919  CHECK(mpoly);
920 
921  if (mpoly->IsEmpty()) {
922  // until the run-time can handle empties
923  throw GeoTypesError("MultiPolygon", "'EMPTY' not supported");
924  // return null bounds
925  bounds.push_back(NULL_DOUBLE);
926  bounds.push_back(NULL_DOUBLE);
927  bounds.push_back(NULL_DOUBLE);
928  bounds.push_back(NULL_DOUBLE);
929  return;
930  }
931 
932  BoundingBox bbox;
933  for (auto p = 0; p < mpoly->getNumGeometries(); p++) {
934  const auto mpoly_geom = mpoly->getGeometryRef(p);
935  CHECK(mpoly_geom);
936  const auto poly_geom = dynamic_cast<OGRPolygon*>(mpoly_geom);
937  if (!poly_geom) {
938  throw GeoTypesError("MultiPolygon",
939  "Failed to read polygon geometry from multipolygon");
940  }
941  const auto exterior_ring = poly_geom->getExteriorRing();
942  CHECK(exterior_ring);
943  // All exterior rings are imported CCW
944  if (exterior_ring->isClockwise()) {
945  exterior_ring->reverseWindingOrder();
946  }
947  const auto num_points_added = process_poly_ring(exterior_ring, coords, &bbox);
948  ring_sizes.push_back(num_points_added);
949 
950  for (auto r = 0; r < poly_geom->getNumInteriorRings(); r++) {
951  auto interior_ring = poly_geom->getInteriorRing(r);
952  CHECK(interior_ring);
953  // All interior rings are imported CW
954  if (!interior_ring->isClockwise()) {
955  interior_ring->reverseWindingOrder();
956  }
957  const auto num_points_added = process_poly_ring(interior_ring, coords, nullptr);
958  ring_sizes.push_back(num_points_added);
959  }
960  poly_rings.push_back(poly_geom->getNumInteriorRings() + 1);
961  }
962  bounds.push_back(bbox.min.x);
963  bounds.push_back(bbox.min.y);
964  bounds.push_back(bbox.max.x);
965  bounds.push_back(bbox.max.y);
966 }
967 
968 std::unique_ptr<GeoBase> GeoGeometry::clone() const {
969  CHECK(geom_);
970  return std::unique_ptr<GeoBase>(new GeoGeometry(geom_->clone(), true));
971 }
972 
973 std::unique_ptr<GeoBase> GeoGeometryCollection::clone() const {
974  CHECK(geom_);
975  return std::unique_ptr<GeoBase>(new GeoGeometryCollection(geom_->clone(), true));
976 }
977 
979  const auto err = GeoBase::createFromWktString(wkt, &geom_);
980  if (err != OGRERR_NONE) {
981  throw GeoTypesError("GeometryCollection", err);
982  }
983  CHECK(geom_);
984  if (wkbFlatten(geom_->getGeometryType()) != OGRwkbGeometryType::wkbGeometryCollection) {
985  throw GeoTypesError("GeometryCollection",
986  "Unexpected geometry type from WKT string: " +
987  std::string(OGRGeometryTypeToName(geom_->getGeometryType())));
988  }
989 }
990 
991 namespace {
992 
994  uint8_t table_[128];
995  constexpr HexDigitToDecimalTable() : table_{} {
996  table_[static_cast<int>('1')] = 1;
997  table_[static_cast<int>('2')] = 2;
998  table_[static_cast<int>('3')] = 3;
999  table_[static_cast<int>('4')] = 4;
1000  table_[static_cast<int>('5')] = 5;
1001  table_[static_cast<int>('6')] = 6;
1002  table_[static_cast<int>('7')] = 7;
1003  table_[static_cast<int>('8')] = 8;
1004  table_[static_cast<int>('9')] = 9;
1005  table_[static_cast<int>('a')] = 10;
1006  table_[static_cast<int>('A')] = 10;
1007  table_[static_cast<int>('b')] = 11;
1008  table_[static_cast<int>('B')] = 11;
1009  table_[static_cast<int>('c')] = 12;
1010  table_[static_cast<int>('C')] = 12;
1011  table_[static_cast<int>('d')] = 13;
1012  table_[static_cast<int>('D')] = 13;
1013  table_[static_cast<int>('e')] = 14;
1014  table_[static_cast<int>('E')] = 14;
1015  table_[static_cast<int>('f')] = 15;
1016  table_[static_cast<int>('F')] = 15;
1017  }
1018  constexpr uint8_t operator[](const char& hex_digit) const {
1019  return (hex_digit < 0) ? 0 : table_[static_cast<int>(hex_digit)];
1020  }
1021 };
1022 
1024 
1025 inline uint8_t hex_to_binary(const char& usb, const char& lsb) {
1026  return (hex_digit_to_decimal_table[usb] << 4) | hex_digit_to_decimal_table[lsb];
1027 }
1028 
1029 std::vector<uint8_t> hex_string_to_binary_vector(const std::string& wkb_hex) {
1030  auto num_bytes = wkb_hex.size() >> 1;
1031  std::vector<uint8_t> wkb(num_bytes);
1032  auto* chars = wkb_hex.data();
1033  auto* bytes = wkb.data();
1034  for (size_t i = 0; i < num_bytes; i++) {
1035  auto const& usb = *chars++;
1036  auto const& lsb = *chars++;
1037  *bytes++ = hex_to_binary(usb, lsb);
1038  }
1039  return wkb;
1040 }
1041 
1042 } // namespace
1043 
1044 OGRGeometry* GeoTypesFactory::createOGRGeometry(const std::string& wkt_or_wkb_hex) {
1045  OGRGeometry* geom = nullptr;
1046  OGRErr err = OGRERR_NONE;
1047  if (wkt_or_wkb_hex.empty()) {
1048  err = OGRERR_NOT_ENOUGH_DATA;
1049  } else if (wkt_or_wkb_hex[0] == '0') { // all WKB hex strings start with a 0
1050  auto const wkb = hex_string_to_binary_vector(wkt_or_wkb_hex);
1051  err = GeoBase::createFromWkbView(&geom, {wkb.data(), wkb.size()});
1052  } else {
1053  err = GeoBase::createFromWktString(wkt_or_wkb_hex, &geom);
1054  }
1055  if (err != OGRERR_NONE) {
1056  throw GeoTypesError("GeoFactory", err);
1057  }
1058  return geom;
1059 }
1060 
1061 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoType(
1062  const std::string& wkt_or_wkb_hex) {
1063  return GeoTypesFactory::createGeoTypeImpl(createOGRGeometry(wkt_or_wkb_hex));
1064 }
1065 
1066 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoType(const WkbView wkb_view) {
1067  OGRGeometry* geom = nullptr;
1068  const auto err = GeoBase::createFromWkbView(&geom, wkb_view);
1069  if (err != OGRERR_NONE) {
1070  throw GeoTypesError("GeoFactory", err);
1071  }
1073 }
1074 
1075 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoType(OGRGeometry* geom) {
1076  return GeoTypesFactory::createGeoTypeImpl(geom, false);
1077 }
1078 
1079 bool GeoTypesFactory::getGeoColumns(const std::string& wkt_or_wkb_hex,
1080  SQLTypeInfo& ti,
1081  std::vector<double>& coords,
1082  std::vector<double>& bounds,
1083  std::vector<int>& ring_sizes,
1084  std::vector<int>& poly_rings,
1085  const bool promote_poly_to_mpoly) {
1086  try {
1087  if (wkt_or_wkb_hex.empty() || wkt_or_wkb_hex == "NULL") {
1089  ti, coords, bounds, ring_sizes, poly_rings, promote_poly_to_mpoly);
1090  return true;
1091  }
1092 
1093  const auto geospatial_base = GeoTypesFactory::createGeoType(wkt_or_wkb_hex);
1094 
1095  if (!geospatial_base || !geospatial_base->transform(ti)) {
1096  return false;
1097  }
1098 
1099  getGeoColumnsImpl(geospatial_base,
1100  ti,
1101  coords,
1102  bounds,
1103  ring_sizes,
1104  poly_rings,
1105  promote_poly_to_mpoly);
1106 
1107  } catch (const std::exception& e) {
1108  LOG(ERROR) << "Geospatial Import Error: " << e.what();
1109  return false;
1110  }
1111 
1112  return true;
1113 }
1114 
1115 bool GeoTypesFactory::getGeoColumns(const std::vector<uint8_t>& wkb,
1116  SQLTypeInfo& ti,
1117  std::vector<double>& coords,
1118  std::vector<double>& bounds,
1119  std::vector<int>& ring_sizes,
1120  std::vector<int>& poly_rings,
1121  const bool promote_poly_to_mpoly) {
1122  try {
1123  const auto geospatial_base = GeoTypesFactory::createGeoType({wkb.data(), wkb.size()});
1124 
1125  if (!geospatial_base || !geospatial_base->transform(ti)) {
1126  return false;
1127  }
1128 
1129  getGeoColumnsImpl(geospatial_base,
1130  ti,
1131  coords,
1132  bounds,
1133  ring_sizes,
1134  poly_rings,
1135  promote_poly_to_mpoly);
1136 
1137  } catch (const std::exception& e) {
1138  LOG(ERROR) << "Geospatial Import Error: " << e.what();
1139  return false;
1140  }
1141 
1142  return true;
1143 }
1144 
1145 bool GeoTypesFactory::getGeoColumns(OGRGeometry* geom,
1146  SQLTypeInfo& ti,
1147  std::vector<double>& coords,
1148  std::vector<double>& bounds,
1149  std::vector<int>& ring_sizes,
1150  std::vector<int>& poly_rings,
1151  const bool promote_poly_to_mpoly) {
1152  try {
1153  const auto geospatial_base = GeoTypesFactory::createGeoType(geom);
1154 
1155  if (!geospatial_base || !geospatial_base->transform(ti)) {
1156  return false;
1157  }
1158 
1159  getGeoColumnsImpl(geospatial_base,
1160  ti,
1161  coords,
1162  bounds,
1163  ring_sizes,
1164  poly_rings,
1165  promote_poly_to_mpoly);
1166 
1167  } catch (const std::exception& e) {
1168  LOG(ERROR) << "Geospatial Import Error: " << e.what();
1169  return false;
1170  }
1171 
1172  return true;
1173 }
1174 
1175 bool GeoTypesFactory::getGeoColumns(const std::vector<std::string>* wkt_or_wkb_hex_column,
1176  SQLTypeInfo& ti,
1177  std::vector<std::vector<double>>& coords_column,
1178  std::vector<std::vector<double>>& bounds_column,
1179  std::vector<std::vector<int>>& ring_sizes_column,
1180  std::vector<std::vector<int>>& poly_rings_column,
1181  const bool promote_poly_to_mpoly) {
1182  try {
1183  for (const auto& wkt_or_wkb_hex : *wkt_or_wkb_hex_column) {
1184  std::vector<double> coords;
1185  std::vector<double> bounds;
1186  std::vector<int> ring_sizes;
1187  std::vector<int> poly_rings;
1188 
1189  SQLTypeInfo row_ti{ti};
1190  getGeoColumns(wkt_or_wkb_hex,
1191  row_ti,
1192  coords,
1193  bounds,
1194  ring_sizes,
1195  poly_rings,
1196  promote_poly_to_mpoly);
1197  if (ti.get_type() != row_ti.get_type()) {
1198  if (!promote_poly_to_mpoly || !(row_ti.get_type() == SQLTypes::kPOLYGON &&
1199  ti.get_type() == SQLTypes::kMULTIPOLYGON)) {
1200  throw GeoTypesError("GeoFactory", "Columnar: Geometry type mismatch");
1201  }
1202  }
1203  coords_column.push_back(coords);
1204  bounds_column.push_back(bounds);
1205  ring_sizes_column.push_back(ring_sizes);
1206  poly_rings_column.push_back(poly_rings);
1207  }
1208 
1209  } catch (const std::exception& e) {
1210  LOG(ERROR) << "Geospatial column Import Error: " << e.what();
1211  return false;
1212  }
1213 
1214  return true;
1215 }
1216 
1217 std::unique_ptr<GeoBase> GeoTypesFactory::createGeoTypeImpl(OGRGeometry* geom,
1218  const bool owns_geom_obj) {
1219  switch (wkbFlatten(geom->getGeometryType())) {
1220  case wkbPoint:
1221  return std::unique_ptr<GeoPoint>(new GeoPoint(geom, owns_geom_obj));
1222  case wkbMultiPoint:
1223  return std::unique_ptr<GeoMultiPoint>(new GeoMultiPoint(geom, owns_geom_obj));
1224  case wkbLineString:
1225  return std::unique_ptr<GeoLineString>(new GeoLineString(geom, owns_geom_obj));
1226  case wkbMultiLineString:
1227  return std::unique_ptr<GeoMultiLineString>(
1228  new GeoMultiLineString(geom, owns_geom_obj));
1229  case wkbPolygon:
1230  return std::unique_ptr<GeoPolygon>(new GeoPolygon(geom, owns_geom_obj));
1231  case wkbMultiPolygon:
1232  return std::unique_ptr<GeoMultiPolygon>(new GeoMultiPolygon(geom, owns_geom_obj));
1233  case wkbGeometryCollection:
1234  return std::unique_ptr<GeoGeometryCollection>(
1235  new GeoGeometryCollection(geom, owns_geom_obj));
1236  default:
1237  throw GeoTypesError(
1238  "GeoTypesFactory",
1239  "Unrecognized geometry type: " + std::string(geom->getGeometryName()));
1240  }
1241 }
1242 
1243 void GeoTypesFactory::getGeoColumnsImpl(const std::unique_ptr<GeoBase>& geospatial_base,
1244  SQLTypeInfo& ti,
1245  std::vector<double>& coords,
1246  std::vector<double>& bounds,
1247  std::vector<int>& ring_sizes,
1248  std::vector<int>& poly_rings,
1249  const bool promote_poly_to_mpoly) {
1250  switch (geospatial_base->getType()) {
1251  case GeoBase::GeoType::kPOINT: {
1252  const auto geospatial_point = dynamic_cast<GeoPoint*>(geospatial_base.get());
1253  CHECK(geospatial_point);
1254  geospatial_point->getColumns(coords);
1255  ti.set_type(kPOINT);
1256  break;
1257  }
1259  const auto geospatial_multipoint =
1260  dynamic_cast<GeoMultiPoint*>(geospatial_base.get());
1261  CHECK(geospatial_multipoint);
1262  geospatial_multipoint->getColumns(coords, bounds);
1263  ti.set_type(kMULTIPOINT);
1264  break;
1265  }
1267  const auto geospatial_linestring =
1268  dynamic_cast<GeoLineString*>(geospatial_base.get());
1269  CHECK(geospatial_linestring);
1270  geospatial_linestring->getColumns(coords, bounds);
1271  ti.set_type(kLINESTRING);
1272  break;
1273  }
1275  const auto geospatial_multilinestring =
1276  dynamic_cast<GeoMultiLineString*>(geospatial_base.get());
1277  CHECK(geospatial_multilinestring);
1278  geospatial_multilinestring->getColumns(coords, ring_sizes, bounds);
1280  break;
1281  }
1283  const auto geospatial_poly = dynamic_cast<GeoPolygon*>(geospatial_base.get());
1284  CHECK(geospatial_poly);
1285  geospatial_poly->getColumns(coords, ring_sizes, bounds);
1286  if (promote_poly_to_mpoly) {
1287  if (ring_sizes.size()) {
1288  CHECK_GT(coords.size(), 0u);
1289  poly_rings.push_back(1 + geospatial_poly->getNumInteriorRings());
1290  }
1291  }
1292  ti.set_type(kPOLYGON);
1293  break;
1294  }
1296  const auto geospatial_mpoly = dynamic_cast<GeoMultiPolygon*>(geospatial_base.get());
1297  CHECK(geospatial_mpoly);
1298  geospatial_mpoly->getColumns(coords, ring_sizes, poly_rings, bounds);
1299  ti.set_type(kMULTIPOLYGON);
1300  break;
1301  }
1304  default:
1305  throw std::runtime_error("Unrecognized geospatial type");
1306  }
1307 }
1308 
1310  std::vector<double>& coords,
1311  std::vector<double>& bounds,
1312  std::vector<int>& ring_sizes,
1313  std::vector<int>& poly_rings,
1314  const bool promote_poly_to_mpoly) {
1315  auto t = ti.get_type();
1316  switch (t) {
1317  case kPOINT: {
1318  // NULL fixlen coords array
1319  coords.push_back(NULL_ARRAY_DOUBLE);
1320  coords.push_back(NULL_DOUBLE);
1321  } break;
1322  case kMULTIPOINT:
1323  case kLINESTRING:
1324  case kMULTILINESTRING:
1325  case kPOLYGON:
1326  case kMULTIPOLYGON: {
1327  // Leaving coords array empty
1328  // NULL fixlen bounds array
1329  bounds.push_back(NULL_ARRAY_DOUBLE);
1330  bounds.push_back(NULL_DOUBLE);
1331  bounds.push_back(NULL_DOUBLE);
1332  bounds.push_back(NULL_DOUBLE);
1333  // Leaving ring_sizes and poly_rings arrays empty
1334  } break;
1335  default:
1336  throw std::runtime_error("Unsupported NULL geo");
1337  }
1338 }
1339 
1340 } // namespace Geospatial
OGRGeometry * geom_
Definition: Types.h:102
std::unique_ptr< GeoBase > optimized_run(GeoOp op, const GeoBase &other) const
Definition: Types.cpp:439
static std::unique_ptr< GeoBase > createGeoType(const std::string &wkt_or_wkb_hex)
Definition: Types.cpp:1061
GeoPoint(const std::vector< double > &coords)
Definition: Types.cpp:541
#define NULL_DOUBLE
void getColumns(std::vector< double > &coords, std::vector< double > &bounds) const
Definition: Types.cpp:613
void getColumns(std::vector< double > &coords) const
Definition: Types.cpp:567
#define SRID_LAEA_START
Definition: Types.cpp:222
int32_t getNumInteriorRings() const
Definition: Types.cpp:864
#define SRID_NORTH_UTM_END
Definition: Types.cpp:212
std::unique_ptr< GeoBase > clone() const final
Definition: Types.cpp:585
uint8_t const * ptr_
Definition: Types.h:32
GeoMultiPolygon(const std::vector< double > &coords, const std::vector< int32_t > &ring_sizes, const std::vector< int32_t > &poly_rings)
Definition: Types.cpp:875
bool owns_geom_obj_
Definition: Types.h:103
constexpr auto DOUBLE_MAX
Definition: Types.cpp:40
#define LOG(tag)
Definition: Logger.h:216
constexpr HexDigitToDecimalTable hex_digit_to_decimal_table
Definition: Types.cpp:1023
GeoLineString(const std::vector< double > &coords)
Definition: Types.cpp:650
GeoMultiLineString(const std::vector< double > &coords, const std::vector< int32_t > &linestring_sizes)
Definition: Types.cpp:709
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: Types.h:265
static void init()
Definition: GDAL.cpp:67
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: Types.cpp:1243
Constants for Builtin SQL Types supported by HEAVY.AI.
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: Types.cpp:1309
HOST DEVICE SQLTypes get_type() const
Definition: sqltypes.h:379
static std::unique_ptr< Geospatial::GeoBase > createGeoTypeImpl(OGRGeometry *geom, const bool owns_geom_obj=true)
Definition: Types.cpp:1217
bool getWkb(std::vector< uint8_t > &) const
Definition: Types.cpp:186
#define CHECK_GT(x, y)
Definition: Logger.h:234
constexpr uint8_t operator[](const char &hex_digit) const
Definition: Types.cpp:1018
static std::shared_ptr< OGRCoordinateTransformation > getTransformation(int32_t srid0, int32_t srid1)
Definition: Types.cpp:307
std::string to_string(char const *&&v)
std::unique_ptr< GeoBase > clone() const final
Definition: Types.cpp:870
bool isEmpty() const
Definition: Types.cpp:196
int process_poly_ring(OGRLinearRing *ring, std::vector< double > &coords, BoundingBox *bbox)
Definition: Types.cpp:70
virtual ~GeoBase()
Definition: Types.cpp:145
static int createFromWkbView(OGRGeometry **geom, WkbView const)
Definition: Types.cpp:173
void update(double x, double y)
Definition: Types.cpp:54
static std::string OGRErrorToStr(const int ogr_err)
Definition: Types.cpp:118
GeoGeometry(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: Types.h:237
#define SRID_NORTH_UTM_START
Definition: Types.cpp:210
std::unique_ptr< GeoBase > clone() const override
Definition: Types.cpp:536
void getColumns(std::vector< double > &coords, std::vector< int32_t > &linestring_sizes, std::vector< double > &bounds) const
Definition: Types.cpp:741
void getColumns(std::vector< double > &coords, std::vector< int32_t > &ring_sizes, std::vector< int32_t > &poly_rings, std::vector< double > &bounds) const
Definition: Types.cpp:914
static bool getGeoColumns(const std::string &wkt_or_wkb_hex, 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: Types.cpp:1079
virtual bool operator==(const GeoBase &other) const
Definition: Types.cpp:200
static OGRGeometry * createOGRGeometry(const std::string &wkt_or_wkb_hex)
Definition: Types.cpp:1044
std::unique_ptr< GeoBase > clone() const final
Definition: Types.cpp:785
static int createFromWktString(const std::string &wkt, OGRGeometry **geom)
Definition: Types.cpp:154
tuple line
Definition: parse_ast.py:10
void getColumns(std::vector< double > &coords, std::vector< double > &bounds) const
Definition: Types.cpp:672
std::mutex transformation_map_mutex_
Definition: Types.cpp:114
std::unique_ptr< GeoBase > clone() const final
Definition: Types.cpp:968
constexpr auto DOUBLE_MIN
Definition: Types.cpp:41
GeoMultiPoint(const std::vector< double > &coords)
Definition: Types.cpp:590
HOST DEVICE int get_input_srid() const
Definition: sqltypes.h:383
bool transform(int32_t srid0, int32_t srid1)
Definition: Types.cpp:383
#define NULL_ARRAY_DOUBLE
#define SRID_SOUTH_UTM_START
Definition: Types.cpp:216
#define SRID_SOUTH_UTM_END
Definition: Types.cpp:218
#define CHECK(condition)
Definition: Logger.h:222
std::string getWktString() const
Definition: Types.cpp:164
#define SRID_SOUTH_LAMBERT
Definition: Types.cpp:220
#define SRID_WORLD_MERCATOR
Definition: Types.cpp:208
#define SRID_LAEA_END
Definition: Types.cpp:224
std::unique_ptr< GeoBase > clone() const final
Definition: Types.cpp:645
uint8_t hex_to_binary(const char &usb, const char &lsb)
Definition: Types.cpp:1025
std::unique_ptr< GeoBase > run(GeoOp op, const GeoBase &other) const
Definition: Types.cpp:406
#define SRID_NORTH_LAMBERT
Definition: Types.cpp:214
std::map< std::tuple< int32_t, int32_t >, std::shared_ptr< OGRCoordinateTransformation > > transformation_map_
Definition: Types.cpp:116
std::unique_ptr< GeoBase > clone() const final
Definition: Types.cpp:973
size_t size_
Definition: Types.h:33
std::unique_ptr< GeoBase > clone() const final
Definition: Types.cpp:704
std::vector< uint8_t > hex_string_to_binary_vector(const std::string &wkb_hex)
Definition: Types.cpp:1029
void getColumns(std::vector< double > &coords, std::vector< int32_t > &ring_sizes, std::vector< double > &bounds) const
Definition: Types.cpp:822
int32_t getBestPlanarSRID() const
Definition: Types.cpp:226
GeoPolygon(const std::vector< double > &coords, const std::vector< int32_t > &ring_sizes)
Definition: Types.cpp:790
HOST DEVICE int get_output_srid() const
Definition: sqltypes.h:385
HOST DEVICE void set_type(SQLTypes t)
Definition: sqltypes.h:489