Spatial join area

src/engine/SpatialJoinAlgorithms.cpp

@@ -39,71 +39,153 @@
 };
 
 // ____________________________________________________________________________
-std::string SpatialJoinAlgorithms::betweenQuotes(
-    std::string extractFrom) const {
+std::string_view SpatialJoinAlgorithms::betweenQuotes(
+    std::string_view extractFrom) const {
   size_t pos1 = extractFrom.find("\"", 0);
   size_t pos2 = extractFrom.find("\"", pos1 + 1);
   if (pos1 != std::string::npos && pos2 != std::string::npos) {
     return extractFrom.substr(pos1 + 1, pos2 - pos1 - 1);
   } else {
     return extractFrom;
   }
 }
 
-struct DistanceVisitor : public boost::static_visitor<double> {
-  template <typename Geometry1, typename Geometry2>
-  double operator()(const Geometry1& geom1, const Geometry2& geom2) const {
-    return bg::distance(geom1, geom2);
-  }
+bool SpatialJoinAlgorithms::getAnyGeometry(const IdTable* idtable, size_t row,
+                                           size_t col,
+                                           AnyGeometry& geometry) const {
+  auto printWarning = [alreadyWarned = false,
+                       &spatialJoin = spatialJoin_]() mutable {
+    if (!alreadyWarned) {
+      std::string warning =
+          "The input to a spatial join contained at least one element, "
+          "that is not a point or polygon geometry and is thus skipped. Note "
+          "that QLever currently only accepts point or polygon geometries for "
+          "the spatial joins";
+      AD_LOG_WARN << warning << std::endl;
+      alreadyWarned = true;
+      if (spatialJoin.has_value()) {
+        AD_CORRECTNESS_CHECK(spatialJoin.value() != nullptr);
+        spatialJoin.value()->addWarning(warning);
+      }
+    }
+  };
+
+  std::string str(betweenQuotes(ExportQueryExecutionTrees::idToStringAndType(
+                                    qec_->getIndex(), idtable->at(row, col), {})
+                                    .value()
+                                    .first));
+  try {
+    bg::read_wkt(str, geometry);
+  } catch (const std::exception& e) {
+    printWarning();
+    return false;
+  }
+  return true;
+}
+
+// ____________________________________________________________________________
+double SpatialJoinAlgorithms::convertDegreesToMeters(
+    AnyGeometry geometry1, AnyGeometry geometry2) const {
+  return boost::apply_visitor(DistanceVisitor(), geometry1, geometry2) * 78.630;
+};
+
+// ____________________________________________________________________________
+Point SpatialJoinAlgorithms::convertGeoPointToPoint(GeoPoint point) const {
+  return Point(point.getLng(), point.getLat());
 };
 
+// ____________________________________________________________________________
+Id SpatialJoinAlgorithms::computeDist(const rtreeEntry& geo1,
+                                      const rtreeEntry& geo2) const {
+  auto convertPoint = [&](AnyGeometry geometry, std::optional<Box> bbox) {
+    Box areaBox;
+    if (bbox.has_value()) {
+      areaBox = bbox.value();
+    } else {
+      areaBox = boost::apply_visitor(BoundingBoxVisitor(), geometry);
+    }
+    Point p = calculateMidpointOfBox(areaBox);
+    return GeoPoint(p.get<1>(), p.get<0>());
+  };
+  // use the already parsed geometries to calculate the distance
+  if (useMidpointForAreas_) {
+    auto point1 = geo1.geoPoint;
+    auto point2 = geo2.geoPoint;
+    if (!point1) {
+      point1 = convertPoint(geo1.geometry.value(), geo1.boundingBox);
+    }
+    if (!point2) {
+      point2 = convertPoint(geo2.geometry.value(), geo2.boundingBox);
+    }
+    if (!point1.has_value() || !point2.has_value()) {
+      return Id::makeUndefined();
+    }
+    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()) {
+      return Id::makeFromDouble(
+          convertDegreesToMeters(geo1.geometry.value(), geo2.geometry.value()));
+    } else {
+      // one point and one area
+      AnyGeometry geom1 = geo1.geometry.value();
+      AnyGeometry geom2 = geo2.geometry.value();
+      if (geo1.geoPoint.has_value()) {
+        geom1 = convertGeoPointToPoint(geo1.geoPoint.value());
+      } else if (geo2.geoPoint.has_value()) {
+        geom2 = convertGeoPointToPoint(geo2.geoPoint.value());
+      }
+      return Id::makeFromDouble(convertDegreesToMeters(geom1, geom2));
+    }
+  }
+}
+
 // ____________________________________________________________________________
 Id SpatialJoinAlgorithms::computeDist(const IdTable* idTableLeft,
                                       const IdTable* idTableRight,
                                       size_t rowLeft, size_t rowRight,
                                       ColumnIndex leftPointCol,
                                       ColumnIndex rightPointCol) const {
-  auto getAreaString = [&](const IdTable* idtable, size_t row, size_t col) {
-    return betweenQuotes(ExportQueryExecutionTrees::idToStringAndType(
-                             qec_->getIndex(), idtable->at(row, col), {})
-                             .value()
-                             .first);
-  };
-
   auto getAreaOrPointGeometry =
       [&](const IdTable* idtable, size_t row, size_t col,
           std::optional<GeoPoint> point) -> std::optional<AnyGeometry> {
     AnyGeometry geometry;
-    try {
-      if (!point) {
-        boost::geometry::read_wkt(getAreaString(idtable, row, col), geometry);
-      } else {
-        geometry = Point(point.value().getLng(), point.value().getLat());
+
+    if (!point) {
+      if (!getAnyGeometry(idtable, row, col, geometry)) {
+        // 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;
       }
-    } catch (...) {
-      return std::nullopt;
+    } else {
+      geometry = convertGeoPointToPoint(point.value());
     }
+
     return geometry;
   };
 
   // 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::string areastring = getAreaString(idtable, row, col);
-
-    Box areaBox;
-    try {
-      areaBox = calculateBoundingBoxOfArea(areastring);
-    } catch (...) {
+    AnyGeometry geometry;
+    if (!getAnyGeometry(idtable, row, col, geometry)) {
       // 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);
+
     Point p = calculateMidpointOfBox(areaBox);
     return GeoPoint(p.get<1>(), p.get<0>());
   };
 
+  // if the geometries have not been read already, get them first
   auto point1 = getPoint(idTableLeft, rowLeft, leftPointCol);
   auto point2 = getPoint(idTableRight, rowRight, rightPointCol);
   if (useMidpointForAreas_) {
@@ -126,14 +208,10 @@
     auto geometry2 =
         getAreaOrPointGeometry(idTableRight, rowRight, rightPointCol, point2);
     if (!geometry1.has_value() || !geometry2.has_value()) {
       return Id::makeUndefined();
     }
-    // convert to m and return. Note that the 78630 is an approximation, which
-    // does not provide accurate results for the poles.
-    return Id::makeFromDouble(boost::apply_visitor(DistanceVisitor(),
-                                                   geometry1.value(),
-                                                   geometry2.value()) *
-                              78630);
+    return Id::makeFromDouble(
+        convertDegreesToMeters(geometry1.value(), geometry2.value()));
   }
 }
 
@@ -528,23 +606,6 @@
   return std::array{northPoleReached, southPoleReached};
 }
 
-struct BoundingBoxVisitor : public boost::static_visitor<Box> {
-  template <typename Geometry>
-  Box operator()(const Geometry& geometry) const {
-    Box box;
-    bg::envelope(geometry, box);
-    return box;
-  }
-};
-
-// ____________________________________________________________________________
-Box SpatialJoinAlgorithms::calculateBoundingBoxOfArea(
-    const std::string& wktString) const {
-  AnyGeometry geometry;
-  boost::geometry::read_wkt(wktString, geometry);
-  return boost::apply_visitor(BoundingBoxVisitor(), geometry);
-}
-
 // ____________________________________________________________________________
 Point SpatialJoinAlgorithms::calculateMidpointOfBox(const Box& box) const {
   double lng = (box.min_corner().get<0>() + box.max_corner().get<0>()) / 2.0;
@@ -568,23 +629,6 @@
 
 // ____________________________________________________________________________
 Result SpatialJoinAlgorithms::BoundingBoxAlgorithm() {
-  auto printWarning = [alreadyWarned = false,
-                       &spatialJoin = spatialJoin_]() mutable {
-    if (!alreadyWarned) {
-      std::string warning =
-          "The input to a spatial join contained at least one element, "
-          "that is not a point or polygon geometry and is thus skipped. Note "
-          "that QLever currently only accepts point or polygon geometries for "
-          "the spatial joins";
-      AD_LOG_WARN << warning << std::endl;
-      alreadyWarned = true;
-      if (spatialJoin.has_value()) {
-        AD_CORRECTNESS_CHECK(spatialJoin.value() != nullptr);
-        spatialJoin.value()->addWarning(warning);
-      }
-    }
-  };
-
   const auto [idTableLeft, resultLeft, idTableRight, resultRight, leftJoinCol,
               rightJoinCol, rightSelectedCols, numColumns, maxDist,
               maxResults] = params_;
@@ -611,28 +655,30 @@
     // get point of row i
     auto geopoint = getPoint(smallerResult, i, smallerResJoinCol);
     Box bbox;
+    AnyGeometry geometry;
+    rtreeEntry entry{i, std::nullopt, std::nullopt, std::nullopt};
 
     if (!geopoint) {
-      try {
-        std::string areastring = betweenQuotes(
-            ExportQueryExecutionTrees::idToStringAndType(
-                qec_->getIndex(), smallerResult->at(i, smallerResJoinCol), {})
-                .value()
-                .first);
-        bbox = calculateBoundingBoxOfArea(areastring);
-      } catch (...) {
-        printWarning();
+      if (!getAnyGeometry(smallerResult, i, smallerResJoinCol, geometry)) {
+        // 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
         continue;
       }
+      bbox = boost::apply_visitor(BoundingBoxVisitor(), geometry);
+      entry.geometry = std::move(geometry);
+      entry.boundingBox = bbox;
+
     } else {
+      entry.geoPoint = geopoint.value();
       // create a box with a side length of at most 1mm to approximate the point
       bbox = Box(Point(geopoint.value().getLng(), geopoint.value().getLat()),
                  Point(geopoint.value().getLng() + 0.00000001,
                        geopoint.value().getLat() + 0.00000001));
     }
 
     // add every box together with the row number into the rtree
-    rtree.insert(std::make_pair(std::move(bbox), i));
+    rtree.insert(std::make_pair(std::move(bbox), std::move(entry)));
   }
 
   // query rtree with the other child
@@ -641,24 +687,26 @@
   for (size_t i = 0; i < otherResult->numRows(); i++) {
     auto geopoint = getPoint(otherResult, i, otherResJoinCol);
     std::vector<Box> bbox;
+    rtreeEntry entry{i, std::nullopt, std::nullopt, std::nullopt};
+    entry.row = i;
 
     if (!geopoint) {
-      try {
-        std::string areastring = betweenQuotes(
-            ExportQueryExecutionTrees::idToStringAndType(
-                qec_->getIndex(), otherResult->at(i, otherResJoinCol), {})
-                .value()
-                .first);
-        auto areaBox = calculateBoundingBoxOfArea(areastring);
-        auto midpoint = calculateMidpointOfBox(areaBox);
-        bbox = computeBoundingBox(
-            midpoint,
-            getMaxDistFromMidpointToAnyPointInsideTheBox(areaBox, midpoint));
-      } catch (...) {
-        printWarning();
+      AnyGeometry geometry;
+      if (!getAnyGeometry(otherResult, i, otherResJoinCol, geometry)) {
+        // 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
         continue;
       }
+      auto areaBox = boost::apply_visitor(BoundingBoxVisitor(), geometry);
+      entry.geometry = std::move(geometry);
+      auto midpoint = calculateMidpointOfBox(areaBox);
+      bbox = computeBoundingBox(
+          midpoint,
+          getMaxDistFromMidpointToAnyPointInsideTheBox(areaBox, midpoint));
+
     } else {
+      entry.geoPoint = geopoint;
       bbox = computeBoundingBox(
           Point(geopoint.value().getLng(), geopoint.value().getLat()));
     }
@@ -670,13 +718,12 @@
     });
 
     ql::ranges::for_each(results, [&](const Value& res) {
-      size_t rowLeft = res.second;
+      size_t rowLeft = res.second.row;
       size_t rowRight = i;
       if (!leftResSmaller) {
         std::swap(rowLeft, rowRight);
       }
-      auto distance = computeDist(idTableLeft, idTableRight, rowLeft, rowRight,
-                                  leftJoinCol, rightJoinCol);
+      auto distance = computeDist(res.second, entry);
       AD_CORRECTNESS_CHECK(distance.getDatatype() == Datatype::Double);
       if (distance.getDouble() * 1000 <= static_cast<double>(maxDist.value())) {
         addResultTableEntry(&result, idTableLeft, idTableRight, rowLeft,