Spatial join area

src/engine/SpatialJoinAlgorithms.cpp

@@ -46,23 +46,22 @@
   if (pos1 != std::string::npos && pos2 != std::string::npos) {
     return extractFrom.substr(pos1 + 1, pos2 - pos1 - 1);
   } else {
     return extractFrom;
   }
 }
 
-std::optional<AnyGeometry> SpatialJoinAlgorithms::getAnyGeometry(
-    const IdTable* idtable, size_t row, size_t col) const {
-  auto printWarning = [alreadyWarned = false,
-                       &spatialJoin = spatialJoin_]() mutable {
-    if (!alreadyWarned) {
+std::optional<size_t> SpatialJoinAlgorithms::getAnyGeometry(
+    const IdTable* idtable, size_t row, size_t col) {
+  auto printWarning = [this, &spatialJoin = spatialJoin_]() mutable {
+    if (this->nrOfFailedParsedGeometries_ == 0) {
       std::string warning =
           "The input to a spatial join contained at least one element, "
           "that is not a Point, Linestring, Polygon, MultiPoint, "
           "MultiLinestring or MultiPolygon geometry and is thus skipped. Note "
           "that QLever currently only accepts those geometries for "
           "the spatial joins";
       AD_LOG_WARN << warning << std::endl;
-      alreadyWarned = true;
+      this->nrOfFailedParsedGeometries_ += 1;
       if (spatialJoin.has_value()) {
         AD_CORRECTNESS_CHECK(spatialJoin.value() != nullptr);
         spatialJoin.value()->addWarning(warning);
@@ -82,32 +81,37 @@
   AnyGeometry geometry;
   try {
     bg::read_wkt(str, geometry);
+    geometries_.push_back(std::move(geometry));
   } catch (...) {
     printWarning();
     return std::nullopt;
   }
-  return geometry;
+  return geometries_.size() - 1;  // index of the last element
 }
 
 // ____________________________________________________________________________
-double SpatialJoinAlgorithms::computeDist(const AnyGeometry& geometry1,
-                                          const AnyGeometry& geometry2) const {
-  return boost::apply_visitor(DistanceVisitor(), geometry1, geometry2) * 78.630;
+double SpatialJoinAlgorithms::computeDist(
+    const std::optional<size_t> geometryIndex1,
+    const std::optional<size_t> geometryIndex2) const {
+  return boost::apply_visitor(ClosestPointVisitor(),
+                              geometries_.at(geometryIndex1.value()),
+                              geometries_.at(geometryIndex2.value()));
 };
 
 // ____________________________________________________________________________
-Point SpatialJoinAlgorithms::convertGeoPointToPoint(GeoPoint point) const {
-  return Point(point.getLng(), point.getLat());
+size_t SpatialJoinAlgorithms::convertGeoPointToPoint(GeoPoint point) {
+  geometries_.emplace_back(Point(point.getLng(), point.getLat()));
+  return geometries_.size() - 1;  // index of the last element
 };
 
 // ____________________________________________________________________________
-Id SpatialJoinAlgorithms::computeDist(const rtreeEntry& geo1,
-                                      const rtreeEntry& geo2) const {
-  auto convertPoint = [&](const AnyGeometry& geometry,
-                          const std::optional<Box>& bbox) {
-    Box areaBox;
-    areaBox =
-        bbox.value_or(boost::apply_visitor(BoundingBoxVisitor(), geometry));
+Id SpatialJoinAlgorithms::computeDist(rtreeEntry& geo1, rtreeEntry& geo2) {
+  auto convertPoint = [&](size_t geometryIndex, std::optional<Box> bbox) {
+    if (!bbox.has_value()) {
+      bbox = boost::apply_visitor(BoundingBoxVisitor(),
+                                  geometries_.at(geometryIndex));
+    }
+    const auto& areaBox = bbox.value();
     Point p = calculateMidpointOfBox(areaBox);
     return GeoPoint(p.get<1>(), p.get<0>());
   };
@@ -116,30 +120,32 @@
     auto point1 = geo1.geoPoint_;
     auto point2 = geo2.geoPoint_;
     if (!point1) {
-      point1 = convertPoint(geo1.geometry_.value(), geo1.boundingBox_);
+      point1 =
+          convertPoint(geo1.geometryIndex_.value(), geo1.boundingBox_.value());
     }
     if (!point2) {
-      point2 = convertPoint(geo2.geometry_.value(), geo2.boundingBox_);
+      point2 =
+          convertPoint(geo2.geometryIndex_.value(), geo2.boundingBox_.value());
     }
     return Id::makeFromDouble(
         ad_utility::detail::wktDistImpl(point1.value(), point2.value()));
   } else {
     if (geo1.geoPoint_.has_value() && geo2.geoPoint_.has_value()) {
       return Id::makeFromDouble(ad_utility::detail::wktDistImpl(
           geo1.geoPoint_.value(), geo2.geoPoint_.value()));
-    } else if (geo1.geometry_.has_value() && geo2.geometry_.has_value()) {
+    } else if (geo1.geometryIndex_.has_value() &&
+               geo2.geometryIndex_.has_value()) {
       return Id::makeFromDouble(
-          computeDist(geo1.geometry_.value(), geo2.geometry_.value()));
+          computeDist(geo1.geometryIndex_, geo2.geometryIndex_));
     } else {
       // one point and one area
-      std::optional<AnyGeometry> geom1 = geo1.geometry_;
-      std::optional<AnyGeometry> geom2 = geo2.geometry_;
-      if (geo1.geoPoint_.has_value()) {
-        geom1 = convertGeoPointToPoint(geo1.geoPoint_.value());
-      } else if (geo2.geoPoint_.has_value()) {
-        geom2 = convertGeoPointToPoint(geo2.geoPoint_.value());
+      if (!geo1.geometryIndex_) {
+        geo1.geometryIndex_ = convertGeoPointToPoint(geo1.geoPoint_.value());
+      } else if (!geo2.geometryIndex_) {
+        geo2.geometryIndex_ = convertGeoPointToPoint(geo2.geoPoint_.value());
       }
-      return Id::makeFromDouble(computeDist(geom1.value(), geom2.value()));
+      return Id::makeFromDouble(
+          computeDist(geo1.geometryIndex_, geo2.geometryIndex_));
     }
   }
 }
@@ -149,35 +155,32 @@
                                       const IdTable* idTableRight,
                                       size_t rowLeft, size_t rowRight,
                                       ColumnIndex leftPointCol,
-                                      ColumnIndex rightPointCol) const {
-  auto getAreaOrPointGeometry = [&](const IdTable* idtable, size_t row,
-                                    size_t col, std::optional<GeoPoint> point) {
-    std::optional<AnyGeometry> geometry;
-
+                                      ColumnIndex rightPointCol) {
+  auto getAreaOrPointGeometry =
+      [&](const IdTable* idtable, size_t row, size_t col,
+          std::optional<GeoPoint> point) -> std::optional<size_t> {
     if (!point) {
       // nothing to do. When parsing a point or an area fails, a warning
       // message gets printed at another place and the point/area just gets
       // skipped
-      geometry = getAnyGeometry(idtable, row, col);
-    } else {
-      geometry = convertGeoPointToPoint(point.value());
+      return getAnyGeometry(idtable, row, col);
     }
-
-    return geometry;
+    return convertGeoPointToPoint(point.value());
   };
 
   // for now we need to get the data from the disk, but in the future, this
   // information will be stored in an ID, just like GeoPoint
   auto getAreaPoint = [&](const IdTable* idtable, size_t row,
                           size_t col) -> std::optional<GeoPoint> {
-    std::optional<AnyGeometry> geometry = getAnyGeometry(idtable, row, col);
-    if (!geometry.has_value()) {
+    std::optional<size_t> geometryIndex = getAnyGeometry(idtable, row, col);
+    if (!geometryIndex) {
       // nothing to do. When parsing a point or an area fails, a warning message
       // gets printed at another place and the point/area just gets skipped
       return std::nullopt;
     }
 
-    Box areaBox = boost::apply_visitor(BoundingBoxVisitor(), geometry.value());
+    Box areaBox = boost::apply_visitor(BoundingBoxVisitor(),
+                                       geometries_.at(geometryIndex.value()));
 
     Point p = calculateMidpointOfBox(areaBox);
     return GeoPoint(p.get<1>(), p.get<0>());
@@ -203,15 +206,16 @@
       return Id::makeUndefined();
     }
   } else {
-    entryLeft.geometry_ = getAreaOrPointGeometry(
+    entryLeft.geometryIndex_ = getAreaOrPointGeometry(
         idTableLeft, rowLeft, leftPointCol, entryLeft.geoPoint_);
-    entryRight.geometry_ = getAreaOrPointGeometry(
+    entryRight.geometryIndex_ = getAreaOrPointGeometry(
         idTableRight, rowRight, rightPointCol, entryRight.geoPoint_);
-    if (entryLeft.geometry_ && entryRight.geometry_) {
+    if (entryLeft.geometryIndex_.has_value() &&
+        entryRight.geometryIndex_.has_value()) {
       return computeDist(entryLeft, entryRight);
     } else {
       return Id::makeUndefined();
     }
   }
 }
 
@@ -403,7 +407,7 @@
 }
 
 // ____________________________________________________________________________
-std::vector<Box> SpatialJoinAlgorithms::computeBoundingBox(
+std::vector<Box> SpatialJoinAlgorithms::computeQueryBox(
     const Point& startPoint, double additionalDist) const {
   const auto [idTableLeft, resultLeft, idTableRight, resultRight, leftJoinCol,
               rightJoinCol, rightSelectedCols, numColumns, maxDist,
@@ -434,7 +438,7 @@
   // a single bounding box for the whole planet, do an optimized version
   if (static_cast<double>(maxDist.value()) > circumferenceMax_ / 4.0 &&
       static_cast<double>(maxDist.value()) < circumferenceMax_ / 2.01) {
-    return computeBoundingBoxForLargeDistances(startPoint);
+    return computeQueryBoxForLargeDistances(startPoint);
   }
 
   // compute latitude bound
@@ -487,7 +491,7 @@
 }
 
 // ____________________________________________________________________________
-std::vector<Box> SpatialJoinAlgorithms::computeBoundingBoxForLargeDistances(
+std::vector<Box> SpatialJoinAlgorithms::computeQueryBoxForLargeDistances(
     const Point& startPoint) const {
   const auto [idTableLeft, resultLeft, idTableRight, resultRight, leftJoinCol,
               rightJoinCol, rightSelectedCols, numColumns, maxDist,
@@ -628,29 +632,44 @@
 }
 
 // ____________________________________________________________________________
-Result SpatialJoinAlgorithms::BoundingBoxAlgorithm() {
-  auto getRtreeEntry = [&](const IdTable* idTable, const size_t row,
-                           const ColumnIndex col) -> std::optional<rtreeEntry> {
-    rtreeEntry entry{row, std::nullopt, std::nullopt, std::nullopt};
-    entry.geoPoint_ = getPoint(idTable, row, col);
-
-    if (!entry.geoPoint_) {
-      entry.geometry_ = getAnyGeometry(idTable, row, col);
-      if (!entry.geometry_.has_value()) {
-        return std::nullopt;
-      }
-      entry.boundingBox_ =
-          boost::apply_visitor(BoundingBoxVisitor(), entry.geometry_.value());
-    } else {
-      entry.boundingBox_ =
-          Box(Point(entry.geoPoint_.value().getLng(),
-                    entry.geoPoint_.value().getLat()),
-              Point(entry.geoPoint_.value().getLng() + 0.00000001,
-                    entry.geoPoint_.value().getLat() + 0.00000001));
+std::optional<rtreeEntry> SpatialJoinAlgorithms::getRtreeEntry(
+    const IdTable* idTable, const size_t row, const ColumnIndex col) {
+  rtreeEntry entry{row, std::nullopt, std::nullopt, std::nullopt};
+  entry.geoPoint_ = getPoint(idTable, row, col);
+
+  if (!entry.geoPoint_) {
+    entry.geometryIndex_ = getAnyGeometry(idTable, row, col);
+    if (!entry.geometryIndex_) {
+      return std::nullopt;
     }
-    return entry;
-  };
+    entry.boundingBox_ = boost::apply_visitor(
+        BoundingBoxVisitor(), geometries_.at(entry.geometryIndex_.value()));
+  } else {
+    entry.boundingBox_ =
+        Box(Point(entry.geoPoint_.value().getLng(),
+                  entry.geoPoint_.value().getLat()),
+            Point(entry.geoPoint_.value().getLng() + 0.00000001,
+                  entry.geoPoint_.value().getLat() + 0.00000001));
+  }
+  return entry;
+}
 
+// ____________________________________________________________________________
+std::vector<Box> SpatialJoinAlgorithms::getQueryBox(
+    const std::optional<rtreeEntry>& entry) const {
+  if (!entry.value().geoPoint_) {
+    auto midpoint = calculateMidpointOfBox(entry.value().boundingBox_.value());
+    return computeQueryBox(midpoint,
+                           getMaxDistFromMidpointToAnyPointInsideTheBox(
+                               entry.value().boundingBox_.value(), midpoint));
+  } else {
+    return computeQueryBox(Point(entry.value().geoPoint_.value().getLng(),
+                                 entry.value().geoPoint_.value().getLat()));
+  }
+}
+
+// ____________________________________________________________________________
+Result SpatialJoinAlgorithms::BoundingBoxAlgorithm() {
   const auto [idTableLeft, resultLeft, idTableRight, resultRight, leftJoinCol,
               rightJoinCol, rightSelectedCols, numColumns, maxDist,
               maxResults] = params_;
@@ -683,8 +702,8 @@
       // skipped
       continue;
     }
-    rtree.insert(std::make_pair(entry.value().boundingBox_.value(),
-                                std::move(entry.value())));
+    rtree.insert(std::pair(entry.value().boundingBox_.value(),
+                           std::move(entry.value())));
   }
 
   // query rtree with the other child
@@ -699,26 +718,15 @@
       // skipped
       continue;
     }
-    std::vector<Box> queryBox;
-    if (!entry.value().geoPoint_) {
-      auto midpoint =
-          calculateMidpointOfBox(entry.value().boundingBox_.value());
-      queryBox = computeBoundingBox(
-          midpoint, getMaxDistFromMidpointToAnyPointInsideTheBox(
-                        entry.value().boundingBox_.value(), midpoint));
-    } else {
-      queryBox =
-          computeBoundingBox(Point(entry.value().geoPoint_.value().getLng(),
-                                   entry.value().geoPoint_.value().getLat()));
-    }
+    std::vector<Box> queryBox = getQueryBox(entry);
 
     results.clear();
 
     ql::ranges::for_each(queryBox, [&](const Box& bbox) {
       rtree.query(bgi::intersects(bbox), std::back_inserter(results));
     });
 
-    ql::ranges::for_each(results, [&](const Value& res) {
+    ql::ranges::for_each(results, [&](Value& res) {
       size_t rowLeft = res.second.row_;
       size_t rowRight = i;
       if (!leftResSmaller) {