Faster distance calculation between multigeometries

src/spatialjoin/BoxIds.h

@@ -68,8 +68,8 @@ inline void getBoxIds(const util::geo::I32XSortedLine& line,
 
       const util::geo::I32XSortedPolygon boxPoly{util::geo::I32Polygon(box)};
 
-      auto check = util::geo::intersectsContainsCovers(
-          line, envelope, boxPoly, box, &firstInA, &firstInB);
+      auto check = util::geo::intersectsContainsCovers(line, boxPoly, &firstInA,
+                                                       &firstInB);
 
       if (std::get<0>(check)) {
         if (localXWidth == 1 && localYHeight == 1) {
@@ -118,10 +118,8 @@ inline void getBoxIds(const util::geo::I32XSortedPolygon& poly,
 
       const util::geo::I32XSortedPolygon boxPoly{util::geo::I32Polygon(box)};
 
-      double boxArea = GRID_AREA * (localXWidth) * (localYHeight);
-
-      auto check = util::geo::intersectsContainsCovers(
-          boxPoly, box, boxArea, poly, envelope, area, &firstInA, &firstInB);
+      auto check = util::geo::intersectsContainsCovers(boxPoly, poly, &firstInA,
+                                                       &firstInB);
 
       if (std::get<1>(check)) {
         // we can insert all at once

src/spatialjoin/GeometryCache.cpp

@@ -9,7 +9,7 @@
 #include "GeometryCache.h"
 #include "util/geo/Geo.h"
 
-const static size_t MAX_MEM_CACHE_SIZE = 1 * 1024 * 1024 * 20l;
+const static size_t MAX_MEM_CACHE_SIZE = 3 * 1024 * 1024 * 20l;
 
 // ____________________________________________________________________________
 template <typename W>
@@ -18,7 +18,7 @@ std::shared_ptr<W> sj::GeometryCache<W>::get(size_t off, ssize_t desTid) const {
 
   if (_inMemory) {
     // completely circumvent cache system
-    return std::make_shared<W>(_memStore.at(off));
+    return _memStore.at(off);
   } else if (desTid == -1) {
     // special cache for large geometries
     tid = _numThreads;
@@ -164,10 +164,6 @@ std::pair<size_t, sj::Line> sj::GeometryCache<sj::Line>::getFrom(
   // geom
   estSize += readLine(str, ret.geom);
 
-  // envelope
-  str.read(reinterpret_cast<char*>(&ret.box), sizeof(util::geo::I32Box));
-  estSize += sizeof(util::geo::I32Box);
-
   // id
   uint16_t len;
   str.read(reinterpret_cast<char*>(&len), sizeof(uint16_t));
@@ -182,10 +178,6 @@ std::pair<size_t, sj::Line> sj::GeometryCache<sj::Line>::getFrom(
   str.read(reinterpret_cast<char*>(&ret.subId), sizeof(size_t));
   estSize += sizeof(size_t);
 
-  // length
-  str.read(reinterpret_cast<char*>(&ret.length), sizeof(double));
-  estSize += sizeof(double);
-
   // boxIds
   uint32_t numBoxIds;
   str.read(reinterpret_cast<char*>(&numBoxIds), sizeof(uint32_t));
@@ -216,10 +208,6 @@ std::pair<size_t, sj::Area> sj::GeometryCache<sj::Area>::getFrom(
   // geom
   estSize += readPoly(str, ret.geom);
 
-  // envelope
-  str.read(reinterpret_cast<char*>(&ret.box), sizeof(util::geo::I32Box));
-  estSize += sizeof(util::geo::I32Box);
-
   // id
   uint16_t len;
   str.read(reinterpret_cast<char*>(&len), sizeof(uint16_t));
@@ -234,14 +222,6 @@ std::pair<size_t, sj::Area> sj::GeometryCache<sj::Area>::getFrom(
   str.read(reinterpret_cast<char*>(&ret.subId), sizeof(size_t));
   estSize += sizeof(size_t);
 
-  // area
-  str.read(reinterpret_cast<char*>(&ret.area), sizeof(double));
-  estSize += sizeof(double);
-
-  // outer area
-  str.read(reinterpret_cast<char*>(&ret.outerArea), sizeof(double));
-  estSize += sizeof(double);
-
   // boxIds
   uint32_t numBoxIds;
   str.read(reinterpret_cast<char*>(&numBoxIds), sizeof(uint32_t));
@@ -261,22 +241,8 @@ std::pair<size_t, sj::Area> sj::GeometryCache<sj::Area>::getFrom(
     // simplified inner
     estSize += readPoly(str, ret.inner);
 
-    if (!ret.inner.empty()) {
-      str.read(reinterpret_cast<char*>(&ret.innerBox),
-               sizeof(util::geo::I32Box));
-      str.read(reinterpret_cast<char*>(&ret.innerOuterArea), sizeof(double));
-      estSize += sizeof(double) + sizeof(util::geo::I32Box);
-    }
-
     // simplified outer
     estSize += readPoly(str, ret.outer);
-
-    if (!ret.outer.empty()) {
-      str.read(reinterpret_cast<char*>(&ret.outerBox),
-               sizeof(util::geo::I32Box));
-      str.read(reinterpret_cast<char*>(&ret.outerOuterArea), sizeof(double));
-      estSize += sizeof(double) + sizeof(util::geo::I32Box);
-    }
   }
 
   return {estSize, ret};
@@ -318,13 +284,13 @@ size_t sj::GeometryCache<W>::add(const std::string& raw) {
     // cache for later use
     std::istringstream ss(raw);
     const auto& val = getFrom(0, ss);
-    _memStore[ret] = std::move(val.second);
+    _memStore[ret] = std::make_shared<W>(std::move(val.second));
 
     if (_geomsOffset > MAX_MEM_CACHE_SIZE) {
       _inMemory = false;
 
       for (const auto& val : _memStore) {
-        writeTo(val.second, _geomsF);
+        writeTo(*val.second.get(), _geomsF);
       }
 
       // clear mem store
@@ -404,10 +370,6 @@ size_t sj::GeometryCache<sj::Line>::writeTo(const sj::Line& val,
   // geoms
   ret += writeLine(val.geom, str);
 
-  // envelopes
-  str.write(reinterpret_cast<const char*>(&val.box), sizeof(util::geo::I32Box));
-  ret += sizeof(util::geo::I32Box);
-
   // id
   if (val.id.size() > std::numeric_limits<uint16_t>::max()) {
     throw std::out_of_range(
@@ -425,10 +387,6 @@ size_t sj::GeometryCache<sj::Line>::writeTo(const sj::Line& val,
   str.write(reinterpret_cast<const char*>(&val.subId), sizeof(size_t));
   ret += sizeof(size_t);
 
-  // length
-  str.write(reinterpret_cast<const char*>(&val.length), sizeof(double));
-  ret += sizeof(double);
-
   // boxIds
   uint32_t size = val.boxIds.size();
   str.write(reinterpret_cast<const char*>(&size), sizeof(uint32_t));
@@ -456,10 +414,6 @@ size_t sj::GeometryCache<sj::Area>::writeTo(const sj::Area& val,
   // geoms
   ret += writePoly(val.geom, str);
 
-  // envelope
-  str.write(reinterpret_cast<const char*>(&val.box), sizeof(util::geo::I32Box));
-  ret += sizeof(util::geo::I32Box);
-
   // id
   if (val.id.size() > std::numeric_limits<uint16_t>::max()) {
     throw std::out_of_range(
@@ -469,6 +423,7 @@ size_t sj::GeometryCache<sj::Area>::writeTo(const sj::Area& val,
   uint16_t s = val.id.size();
   str.write(reinterpret_cast<const char*>(&s), sizeof(uint16_t));
   ret += sizeof(uint16_t);
+
   str.write(reinterpret_cast<const char*>(val.id.c_str()),
             val.id.size() * sizeof(char));
   ret += sizeof(char) * val.id.size();
@@ -477,14 +432,6 @@ size_t sj::GeometryCache<sj::Area>::writeTo(const sj::Area& val,
   str.write(reinterpret_cast<const char*>(&val.subId), sizeof(size_t));
   ret += sizeof(size_t);
 
-  // area
-  str.write(reinterpret_cast<const char*>(&val.area), sizeof(double));
-  ret += sizeof(double);
-
-  // outer area
-  str.write(reinterpret_cast<const char*>(&val.outerArea), sizeof(double));
-  ret += sizeof(double);
-
   // boxIds
   uint32_t size = val.boxIds.size();
   str.write(reinterpret_cast<const char*>(&size), sizeof(uint32_t));
@@ -504,30 +451,8 @@ size_t sj::GeometryCache<sj::Area>::writeTo(const sj::Area& val,
     // innerGeom
     ret += writePoly(val.inner, str);
 
-    if (!val.inner.empty()) {
-      str.write(reinterpret_cast<const char*>(&val.innerBox),
-                sizeof(util::geo::I32Box));
-      ret += sizeof(util::geo::I32Box);
-
-      // inner area
-      str.write(reinterpret_cast<const char*>(&val.innerOuterArea),
-                sizeof(double));
-      ret += sizeof(double);
-    }
-
     // outerGeom
     ret += writePoly(val.outer, str);
-
-    if (!val.outer.empty()) {
-      str.write(reinterpret_cast<const char*>(&val.outerBox),
-                sizeof(util::geo::I32Box));
-      ret += sizeof(util::geo::I32Box);
-
-      // outer area
-      str.write(reinterpret_cast<const char*>(&val.outerOuterArea),
-                sizeof(double));
-      ret += sizeof(double);
-    }
   }
 
   return ret;
@@ -547,11 +472,22 @@ template <typename W>
 size_t sj::GeometryCache<W>::readPoly(std::istream& str,
                                       util::geo::I32XSortedPolygon& ret) const {
   size_t estSize = 0;
+
+  util::geo::I32Box box;
+  str.read(reinterpret_cast<char*>(&box), sizeof(util::geo::I32Box));
+  estSize += sizeof(util::geo::I32Box);
+
   double maxSegLen;
   str.read(reinterpret_cast<char*>(&maxSegLen), sizeof(double));
   estSize += sizeof(double);
   ret.getOuter().setMaxSegLen(maxSegLen);
 
+  // outer area
+  double area;
+  str.read(reinterpret_cast<char*>(&area), sizeof(double));
+  estSize += sizeof(double);
+  ret.getOuter().setArea(area);
+
   uint32_t sizeOuter;
   str.read(reinterpret_cast<char*>(&sizeOuter), sizeof(uint32_t));
 
@@ -562,6 +498,8 @@ size_t sj::GeometryCache<W>::readPoly(std::istream& str,
     str.read(reinterpret_cast<char*>(&ret.getOuter().rawRing()[0]),
              sizeof(util::geo::XSortedTuple<int32_t>) * sizeOuter);
     estSize += sizeof(util::geo::XSortedTuple<int32_t>) * sizeOuter;
+
+    ret.getOuter().setBoundingBox(box);
   }
 
   uint32_t numInners;
@@ -609,6 +547,8 @@ size_t sj::GeometryCache<W>::readPoly(std::istream& str,
       str.read(reinterpret_cast<char*>(&ret.getInners()[j].rawRing()[0]),
                sizeof(util::geo::XSortedTuple<int32_t>) * sizeInner);
       estSize += sizeof(util::geo::XSortedTuple<int32_t>) * sizeInner;
+      ret.getInners()[j].setBoundingBox(ret.getInnerBoxes()[j]);
+      ret.getInners()[j].setArea(ret.getInnerAreas()[j]);
     }
   }
 
@@ -621,11 +561,21 @@ size_t sj::GeometryCache<W>::writePoly(const util::geo::I32XSortedPolygon& geom,
                                        std::ostream& str) {
   size_t ret = 0;
 
-  // geom, outer
+  // outer envelope
+  const auto box = geom.boundingBox();
+  str.write(reinterpret_cast<const char*>(&box), sizeof(util::geo::I32Box));
+  ret += sizeof(util::geo::I32Box);
+
+  // max seg len
   double maxSegLen = geom.getOuter().getMaxSegLen();
   str.write(reinterpret_cast<const char*>(&maxSegLen), sizeof(double));
   ret += sizeof(double);
 
+  // outer area
+  double area = geom.getOuter().area();
+  str.write(reinterpret_cast<const char*>(&area), sizeof(double));
+  ret += sizeof(double);
+
   uint32_t locSize = geom.getOuter().rawRing().size();
   str.write(reinterpret_cast<const char*>(&locSize), sizeof(uint32_t));
   if (locSize) {
@@ -685,6 +635,18 @@ template <typename W>
 size_t sj::GeometryCache<W>::readLine(std::istream& str,
                                       util::geo::I32XSortedLine& ret) const {
   size_t estSize = 0;
+
+  util::geo::I32Box box;
+  str.read(reinterpret_cast<char*>(&box), sizeof(util::geo::I32Box));
+  estSize += sizeof(util::geo::I32Box);
+  ret.setBoundingBox(box);
+
+  // length
+  double length;
+  str.read(reinterpret_cast<char*>(&length), sizeof(double));
+  estSize += sizeof(double);
+  ret.setLength(length);
+
   double maxSegLen;
   str.read(reinterpret_cast<char*>(&maxSegLen), sizeof(double));
   estSize += sizeof(double);
@@ -721,6 +683,16 @@ size_t sj::GeometryCache<W>::writeLine(const util::geo::I32XSortedLine& geom,
                                        std::ostream& str) {
   size_t ret = 0;
 
+  // outer envelope
+  const auto box = geom.boundingBox();
+  str.write(reinterpret_cast<const char*>(&box), sizeof(util::geo::I32Box));
+  ret += sizeof(util::geo::I32Box);
+
+  // length
+  double length = geom.length();
+  str.write(reinterpret_cast<const char*>(&length), sizeof(double));
+  ret += sizeof(double);
+
   double maxSegLen = geom.getMaxSegLen();
   str.write(reinterpret_cast<const char*>(&maxSegLen), sizeof(double));
   ret += sizeof(double);

src/spatialjoin/GeometryCache.h

@@ -29,21 +29,12 @@ struct Area {
   // polygons
   util::geo::I32XSortedPolygon geom;
 
-  // envelope
-  util::geo::I32Box box;
-
   // id
   std::string id;
 
   // sub id (for multipolygons)
   size_t subId;
 
-  // area
-  double area;
-
-  // outer area
-  double outerArea;
-
   // box ids
   std::vector<sj::boxids::BoxId> boxIds;
 
@@ -53,20 +44,8 @@ struct Area {
   // inner geom
   util::geo::I32XSortedPolygon inner;
 
-  // inner polygon envelope
-  util::geo::I32Box innerBox;
-
-  // outer area for inner polygon
-  double innerOuterArea;
-
   // outer geom
   util::geo::I32XSortedPolygon outer;
-
-  // outer polygon envelope
-  util::geo::I32Box outerBox;
-
-  // outer area for outer polygon
-  double outerOuterArea;
 };
 
 struct SimpleLine {
@@ -78,18 +57,12 @@ struct Line {
   // line
   util::geo::I32XSortedLine geom;
 
-  // envelope
-  util::geo::I32Box box;
-
   // id
   std::string id;
 
   // sub id (for multilines)
   size_t subId;
 
-  // length
-  double length;
-
   // box ids
   std::vector<sj::boxids::BoxId> boxIds;
 
@@ -217,7 +190,8 @@ class GeometryCache {
   std::string _dir, _tmpPrefix;
   std::string _fName;
 
-  std::map<size_t, W> _memStore;
+  // must be ordered map to ensure correct writing order when flushing to disk!
+  std::map<size_t, std::shared_ptr<W>> _memStore;
   bool _inMemory = true;
 
   char* _writeBuffer = 0;