bz2 and gz input file support, options --stats and --verbose for more output, update util

Author: patrickbr

CMakeLists.txt

@@ -27,10 +27,14 @@ find_package(ZLIB)
 
 if (NOT BZIP2_FOUND)
 	add_definitions( -DSPATIALJOIN_NO_BZIP2=True )
+else()
+	add_definitions( -DPBUTIL_BZIP2_FOUND=True )
 endif()
 
 if (NOT ZLIB_FOUND)
 	add_definitions( -DSPATIALJOIN_NO_ZLIB=True )
+else()
+	add_definitions( -DPBUTIL_ZLIB_FOUND=True )
 endif()
 
 # export compile commands to tools like clang

src/spatialjoin/CMakeLists.txt

@@ -1,9 +1,6 @@
 file(GLOB_RECURSE SPATIALJOIN_SRC *.cpp)
 list(REMOVE_ITEM SPATIALJOIN_SRC TestMain.cpp)
-
-if(NOT CMAKE_TESTING_ENABLED)
-	list(FILTER util_SRC EXCLUDE REGEX ".*/tests/.*")
-endif()
+list(FILTER SPATIALJOIN_SRC EXCLUDE REGEX ".*/tests/.*")
 
 set(spatialjoin_main SpatialJoinMain.cpp)
 
@@ -20,7 +17,7 @@ target_include_directories(spatialjoin-dev PUBLIC
 	$<BUILD_INTERFACE:${SPATIALJOIN_INCLUDE_DIR}>
 )
 
-target_link_libraries(spatialjoin spatialjoin-dev pb_util -lpthread)
+target_link_libraries(spatialjoin spatialjoin-dev pb_util pb_util_geo -lpthread)
 
 if (BZIP2_FOUND)
 	target_include_directories(spatialjoin PUBLIC ${BZIP2_INCLUDE_DIR} )

src/spatialjoin/SpatialJoinMain.cpp

@@ -41,7 +41,11 @@ void printHelp(int argc, char** argv) {
       << "\n"
       << "(C) 2023-" << YEAR << " " << COPY << "\n"
       << "Authors: " << AUTHORS << "\n\n"
-      << "Usage: " << argv[0] << " [--help] [-h] [input]\n\n"
+      << "Usage: " << argv[0] << " [OPTIONS] [INPUT1] [INPUT1]\n\n"
+      << "With input from either stdin, or from file(s) [INPUT1] and [INPUT2] "
+         "(.bz2 or .gz\nsupported). If "
+      << "both [INPUT1] and [INPUT2] are given, compute a non-self join with\n"
+      << "[INPUT1] on the left side, [INPUT2] on the right side.\n\n"
       << "Allowed options:\n\n"
       << std::setfill(' ') << std::left << "General:\n"
       << std::setw(42) << "  -h [ --help ]"
@@ -52,6 +56,11 @@ void printHelp(int argc, char** argv) {
       << "cache directory for intermediate files\n"
       << std::setw(42) << "  --de9im"
       << "output DE-9IM relationships\n"
+      << std::setw(42) << "  --within-distance (default: '')"
+      << "if set to non-negative value, only compute for each object\n"
+      << std::setw(42) << " "
+      << "the objects within the given distance\n\n"
+      << std::setfill(' ') << std::left << "Formatting:\n"
       << std::setw(42) << "  --prefix (default: '')"
       << "prefix added at the beginning of every relation\n"
       << std::setw(42) << "  --intersects (default: ' intersects ')"
@@ -70,12 +79,11 @@ void printHelp(int argc, char** argv) {
       << "separator between crossing geometry IDs\n"
       << std::setw(42) << "  --suffix (default: '\\n')"
       << "suffix added at the beginning of every relation\n\n"
-      << std::setw(42) << "  --within-distance (default: '')"
-      << "if set to non-negative value, only compute for each object the "
-         "objects within the given distance\n\n"
       << std::setfill(' ') << std::left << "Geometric computation:\n"
       << std::setw(42) << "  --no-box-ids"
-      << "disable box id criteria for contains/covers/intersect computation\n"
+      << "disable box id criteria for contains/covers/intersect\n"
+      << std::setw(42) << " "
+      << "computation\n"
       << std::setw(42) << "  --no-surface-area"
       << "disable surface area criteria for polygon contains/covers\n"
       << std::setw(42) << "  --no-oriented-envelope"
@@ -94,11 +102,17 @@ void printHelp(int argc, char** argv) {
       << "  --num-caches (default: " + std::to_string(NUM_THREADS) + ")"
       << "number of geometry caches (if < --num-threads, syncing\n"
       << std::setw(42)
-      << "  --cache-max-size (default: " + std::to_string(DEFAULT_CACHE_SIZE) + ")"
+      << "  --cache-max-size (default: " + std::to_string(DEFAULT_CACHE_SIZE) +
+             ")"
       << "maximum approx. size in bytes of cache per type and thread\n"
       << std::setw(42) << "  --no-geometry-checks"
-      << "do not compute geometric relations, only report number of "
-         "candidates\n"
+      << "do not compute geometric relations, only report number of\n"
+      << std::setw(42) << " "
+      << "candidates\n"
+      << std::setw(42) << "  --stats"
+      << "output stats\n"
+      << std::setw(42) << "  -v [ --verbose ]"
+      << "verbose logging\n"
       << std::endl;
 }
 
@@ -134,6 +148,9 @@ int main(int argc, char** argv) {
   bool noGeometryChecks = false;
   bool computeDE9IM = false;
 
+  bool printStats = false;
+  bool verbose = false;
+
   size_t numThreads = NUM_THREADS;
   size_t numCaches = NUM_THREADS;
   size_t geomCacheMaxSizeBytes = DEFAULT_CACHE_SIZE;
@@ -194,6 +211,10 @@ int main(int argc, char** argv) {
           useFastSweepSkip = false;
         } else if (cur == "--use-inner-outer") {
           useInnerOuter = true;
+        } else if (cur == "--stats") {
+          printStats = true;
+        } else if (cur == "--verbose" || cur == "-v") {
+          verbose = true;
         } else {
           inputFiles.push_back(cur);
         }
@@ -265,53 +286,116 @@ int main(int argc, char** argv) {
   unsigned char* buf = new unsigned char[CACHE_SIZE];
   size_t len;
 
-  Sweeper sweeper({numThreads,
-                   numCaches,
-                   geomCacheMaxSizeBytes / (numThreads * 3),
-                   prefix,
-                   intersects,
-                   contains,
-                   covers,
-                   touches,
-                   equals,
-                   overlaps,
-                   crosses,
-                   suffix,
-                   useBoxIds,
-                   useArea,
-                   useOBB,
-                   useDiagBox,
-                   useFastSweepSkip,
-                   useInnerOuter,
-                   noGeometryChecks,
-                   withinDist,
-                   computeDE9IM,
-                   {},
-                   [](const std::string& s) { LOGTO(INFO, std::cerr) << s; },
-                   [](const std::string& s) { std::cerr << s; },
-                   {},
-                   {}},
-                  cache, output);
-
-  LOGTO(INFO, std::cerr) << "Parsing input geometries...";
+  sj::SweeperCfg sweeperCfg{numThreads,
+                            numCaches,
+                            geomCacheMaxSizeBytes / (numThreads * 3),
+                            prefix,
+                            intersects,
+                            contains,
+                            covers,
+                            touches,
+                            equals,
+                            overlaps,
+                            crosses,
+                            suffix,
+                            useBoxIds,
+                            useArea,
+                            useOBB,
+                            useDiagBox,
+                            useFastSweepSkip,
+                            useInnerOuter,
+                            noGeometryChecks,
+                            withinDist,
+                            computeDE9IM,
+                            {},
+                            {},
+                            {},
+                            {},
+                            {}};
+
+  if (printStats)
+    sweeperCfg.statsCb = [](const std::string& s) { std::cerr << s; };
+
+  if (verbose)
+    sweeperCfg.logCb = [](const std::string& s) {
+      LOGTO(INFO, std::cerr) << s;
+    };
+
+  Sweeper sweeper(sweeperCfg, cache, output);
+
+  sweeper.log("Parsing input geometries...");
   auto ts = TIME();
 
   sj::WKTParser parser(&sweeper, NUM_THREADS);
 
   if (!inputFiles.empty()) {
     if (inputFiles.size() > 2) {
-      std::cerr << "Either 1 input files (for self join), or 2 input files (for non-self join) can be provided." << std::endl;
+      std::cerr << "Either 1 input files (for self join), or 2 input files "
+                   "(for non-self join) can be provided."
+                << std::endl;
       exit(1);
     }
     for (size_t i = 0; i < inputFiles.size(); i++) {
-      int f = open(inputFiles[i].c_str(), O_RDONLY);
+      if (util::endsWith(inputFiles[i], ".bz2")) {
+#ifndef SPATIALJOIN_NO_BZIP2
+        auto fh = fopen(inputFiles[i].c_str(), "r");
+        if (!fh) {
+          std::cerr << "Could not open input file " << inputFiles[i]
+                    << std::endl;
+          exit(1);
+        }
+        int err;
+        BZFILE* f = BZ2_bzReadOpen(&err, fh, 0, 0, NULL, 0);
+        if (!f || err != BZ_OK) {
+          std::cerr << "Could not open input file " << inputFiles[i]
+                    << std::endl;
+          exit(1);
+        }
+        while ((len = util::bz2readAll(f, buf, CACHE_SIZE)) > 0) {
+          parser.parse(reinterpret_cast<char*>(buf), len, i != 0);
+        }
 
-      if (f < 0) {
-        throw std::runtime_error("Could not open input file " + inputFiles[i]);
-      }
+        BZ2_bzReadClose(&err, f);
+        fclose(fh);
+#else
+        std::cerr << "Could not open input file " << inputFiles[i]
+                  << ", spatialjoin was compiled without BZip2 support"
+                  << std::endl;
+        exit(1);
+#endif
+      } else if (util::endsWith(inputFiles[i], ".gz")) {
+#ifndef SPATIALJOIN_NO_ZLIB
+        gzFile f = gzopen(inputFiles[i].c_str(), "r");
+        if (f == Z_NULL) {
+          std::cerr << "Could not open input file " << inputFiles[i]
+                    << std::endl;
+          exit(1);
+        }
+        while ((len = util::zreadAll(f, buf, CACHE_SIZE)) > 0) {
+          parser.parse(reinterpret_cast<char*>(buf), len, i != 0);
+        }
+
+        gzclose(f);
+#else
+        std::cerr << "Could not open input file " << inputFiles[i]
+                  << ", spatialjoin was compiled without gzip support"
+                  << std::endl;
+        exit(1);
+#endif
+      } else {
+        int f = open(inputFiles[i].c_str(), O_RDONLY);
+
+        if (f < 0) {
+          std::cerr << "Could not open input file " << inputFiles[i]
+                    << std::endl;
+          exit(1);
+        }
+
+        while ((len = util::readAll(f, buf, CACHE_SIZE)) > 0) {
+          parser.parse(reinterpret_cast<char*>(buf), len, i != 0);
+        }
 
-      while ((len = util::readAll(f, buf, CACHE_SIZE)) > 0) {
-        parser.parse(reinterpret_cast<char*>(buf), len, i != 0);
+        close(f);
       }
     }
   } else {
@@ -322,20 +406,20 @@ int main(int argc, char** argv) {
 
   parser.done();
 
-  LOGTO(INFO, std::cerr) << "Done parsing ("
-                         << TOOK(ts) / 1000000000.0 << "s).";
+  sweeper.log("Done parsing (" + std::to_string(TOOK(ts) / 1000000000.0) +
+              "s).");
   ts = TIME();
 
-  LOGTO(INFO, std::cerr) << "Sorting sweep events...";
+  sweeper.log("Sorting sweep events...");
 
   sweeper.flush();
 
-  LOGTO(INFO, std::cerr) << "done (" << TOOK(ts) / 1000000000.0 << "s).";
+  sweeper.log("done (" + std::to_string(TOOK(ts) / 1000000000.0) + "s).");
 
-  LOGTO(INFO, std::cerr) << "Sweeping...";
+  sweeper.log("Sweeping...");
   ts = TIME();
   sweeper.sweep();
-  LOGTO(INFO, std::cerr) << "done (" << TOOK(ts) / 1000000000.0 << "s).";
+  sweeper.log("done (" + std::to_string(TOOK(ts) / 1000000000.0) + "s).");
 
   delete[] buf;
 }

src/spatialjoin/Sweeper.h

@@ -252,6 +252,8 @@ class Sweeper {
     }
   }
 
+  void log(const std::string& msg);
+
   util::geo::I32Box add(const util::geo::I32MultiPolygon& a,
                         const std::string& gid, bool side,
                         WriteBatch& batch) const;
@@ -541,8 +543,6 @@ class Sweeper {
   void prepareOutputFiles();
   void flushOutputFiles();
 
-  void log(const std::string& msg);
-
   bool notOverlaps(const std::string& a, const std::string& b);
   bool notTouches(const std::string& a, const std::string& b);
   bool notCrosses(const std::string& a, const std::string& b);

src/spatialjoin/WKTParse.cpp

@@ -19,6 +19,7 @@ WKTParser::WKTParser(sj::Sweeper *sweeper, size_t numThreads)
   }
 }
 
+// _____________________________________________________________________________
 void WKTParser::processQueue(size_t t) {
   ParseBatch batch;
   while ((batch = _jobs.get()).size()) {

src/spatialjoin/WKTParse.h

@@ -137,35 +137,35 @@ class WKTParserBase {
     } else {
       auto wktType = getWKTType(c, &c);
       if (wktType == util::geo::WKTType::POINT) {
-        const auto &point = pointFromWKT<int32_t>(c, 0, &projFunc);
+        const auto &point = pointFromWKTProj<int32_t>(c, 0, &projFunc);
         _bboxes[t] = util::geo::extendBox(_sweeper->add(point, id, side, batch),
                                           _bboxes[t]);
       } else if (wktType == util::geo::WKTType::MULTIPOINT) {
-        const auto &mp = multiPointFromWKT<int32_t>(c, 0, &projFunc);
+        const auto &mp = multiPointFromWKTProj<int32_t>(c, 0, &projFunc);
         if (mp.size() != 0)
           _bboxes[t] = util::geo::extendBox(_sweeper->add(mp, id, side, batch),
                                             _bboxes[t]);
       } else if (wktType == util::geo::WKTType::LINESTRING) {
-        const auto &line = lineFromWKT<int32_t>(c, 0, &projFunc);
+        const auto &line = lineFromWKTProj<int32_t>(c, 0, &projFunc);
         if (line.size() > 1)
           _bboxes[t] = util::geo::extendBox(
               _sweeper->add(line, id, side, batch), _bboxes[t]);
       } else if (wktType == util::geo::WKTType::MULTILINESTRING) {
-        const auto &ml = multiLineFromWKT<int32_t>(c, 0, &projFunc);
+        const auto &ml = multiLineFromWKTProj<int32_t>(c, 0, &projFunc);
         _bboxes[t] = util::geo::extendBox(_sweeper->add(ml, id, side, batch),
                                           _bboxes[t]);
       } else if (wktType == util::geo::WKTType::POLYGON) {
-        const auto &poly = polygonFromWKT<int32_t>(c, 0, &projFunc);
+        const auto &poly = polygonFromWKTProj<int32_t>(c, 0, &projFunc);
         if (poly.getOuter().size() > 1)
           _bboxes[t] = util::geo::extendBox(
               _sweeper->add(poly, id, side, batch), _bboxes[t]);
       } else if (wktType == util::geo::WKTType::MULTIPOLYGON) {
-        const auto &mp = multiPolygonFromWKT<int32_t>(c, 0, &projFunc);
+        const auto &mp = multiPolygonFromWKTProj<int32_t>(c, 0, &projFunc);
         if (mp.size())
           _bboxes[t] = util::geo::extendBox(_sweeper->add(mp, id, side, batch),
                                             _bboxes[t]);
       } else if (wktType == util::geo::WKTType::COLLECTION) {
-        const auto &col = collectionFromWKT<int32_t>(c, 0, &projFunc);
+        const auto &col = collectionFromWKTProj<int32_t>(c, 0, &projFunc);
 
         size_t numGeoms = 0;
         for (const auto &a : col) {

src/spatialjoin/tests/CMakeLists.txt

@@ -1,2 +1,2 @@
 add_executable(spatialjoinTest TestMain.cpp)
-target_link_libraries(spatialjoinTest spatialjoin-dev pb_util ${BZIP2_LIBRARIES} -lpthread)
+target_link_libraries(spatialjoinTest spatialjoin-dev pb_util pb_util_geo ${BZIP2_LIBRARIES} -lpthread)