diff --git a/.gitmodules b/.gitmodules index 4a3187d404..e7cf59b075 100644 --- a/.gitmodules +++ b/.gitmodules @@ -58,7 +58,7 @@ [submodule "extern/ueb-bmi"] path = extern/ueb-bmi url = https://github.com/NGWPC/ueb_bmi.git - branch = development + branch = idt-ueb-inclusive-end-step [submodule "extern/lstm"] path = extern/lstm url = https://github.com/NGWPC/lstm.git diff --git a/src/NGen.cpp b/src/NGen.cpp index 1333733031..c2f9045851 100644 --- a/src/NGen.cpp +++ b/src/NGen.cpp @@ -467,9 +467,10 @@ int main(int argc, char* argv[]) { #if NGEN_WITH_SQLITE3 try { nexus_collection = ngen::geopackage::read(nexusDataFile, "nexus", nexus_subset_ids); - } catch (...) { + } catch (const std::exception &e) { // Handle all exceptions - std::string msg = "Geopackage error occurred reading nexuses: " + nexusDataFile; + std::string msg = "Geopackage error occurred reading nexuses: " + nexusDataFile + + "\nError: " + e.what(); LOG(msg,LogLevel::FATAL); throw std::runtime_error(msg); } @@ -498,9 +499,10 @@ int main(int argc, char* argv[]) { try { catchment_collection = ngen::geopackage::read(catchmentDataFile, "divides", catchment_subset_ids); - } catch (...) { + } catch (const std::exception &e) { // Handle all exceptions - std::string msg = "Geopackage error occurred reading divides: " + catchmentDataFile; + std::string msg = "Geopackage error occurred reading divides: " + catchmentDataFile + + "\nError: " + e.what(); LOG(msg,LogLevel::FATAL); throw std::runtime_error(msg); } diff --git a/src/geopackage/read.cpp b/src/geopackage/read.cpp index 5cb5d9626c..c4031a00a4 100644 --- a/src/geopackage/read.cpp +++ b/src/geopackage/read.cpp @@ -26,7 +26,15 @@ std::shared_ptr ngen::geopackage::read( { // Check for malicious/invalid layer input check_table_name(layer); + std::vector features; + if (ids.size() > 0) + features.reserve(ids.size()); + double min_x = std::numeric_limits::infinity(); + double min_y = std::numeric_limits::infinity(); + double max_x = -std::numeric_limits::infinity(); + double max_y = -std::numeric_limits::infinity(); + LOG(LogLevel::DEBUG, "Establishing connection to geopackage %s.", gpkg_path.c_str()); ngen::sqlite::database db{gpkg_path}; // Check if layer exists @@ -66,80 +74,86 @@ std::shared_ptr ngen::geopackage::read( } } - // Layer exists, getting statement for it - // - // this creates a string in the form: - // WHERE id IN (?, ?, ?, ...) - // so that it can be bound by SQLite. - // This is safer than trying to concatenate - // the IDs together. - std::string joined_ids = ""; - if (!ids.empty()) { - joined_ids = " WHERE "+id_column+" IN (?"; - for (size_t i = 1; i < ids.size(); i++) { - joined_ids += ", ?"; + // execute sub-queries if the number of IDs gets too long or once if ids.size() == 0 + int bind_limit = 900; + boost::span id_span(ids); + for (int i = 0; i < ids.size() || (i == 0 && ids.size() == 0); i += bind_limit) { + int span_size = (i + bind_limit >= ids.size()) ? (ids.size() - i) : bind_limit; + boost::span sub_ids = id_span.subspan(i, span_size); + + // Layer exists, getting statement for it + // + // this creates a string in the form: + // WHERE id IN (?, ?, ?, ...) + // so that it can be bound by SQLite. + // This is safer than trying to concatenate + // the IDs together. + std::string joined_ids = ""; + if (!sub_ids.empty()) { + joined_ids = " WHERE "+id_column+" IN (?"; + for (size_t i = 1; i < sub_ids.size(); i++) { + joined_ids += ", ?"; + } + joined_ids += ")"; } - joined_ids += ")"; - } - // Get number of features - auto query_get_layer_count = db.query("SELECT COUNT(*) FROM " + layer + joined_ids, ids); - query_get_layer_count.next(); - const int layer_feature_count = query_get_layer_count.get(0); - - #ifndef NGEN_QUIET - // output debug info on what is read exactly - read_ss << "Reading " << layer_feature_count << " features from layer " << layer << " using ID column `"<< id_column << "`"; - if (!ids.empty()) { - read_ss << " (id subset:"; - for (auto& id : ids) { - read_ss << " " << id; + // Get number of features + auto query_get_layer_count = db.query("SELECT COUNT(*) FROM " + layer + joined_ids, sub_ids); + query_get_layer_count.next(); + const int layer_feature_count = query_get_layer_count.get(0); + + #ifndef NGEN_QUIET + // output debug info on what is read exactly + read_ss << "Reading " << layer_feature_count << " features from layer " << layer << " using ID column `"<< id_column << "`"; + if (!sub_ids.empty()) { + read_ss << " (id subset:"; + for (auto& id : sub_ids) { + read_ss << " " << id; + } + read_ss << ")"; } - read_ss << ")"; - } - read_ss << std::endl; - LOG(read_ss.str(), LogLevel::DEBUG); read_ss.str(""); - #endif - - // Get layer feature metadata (geometry column name + type) - auto query_get_layer_geom_meta = db.query("SELECT column_name FROM gpkg_geometry_columns WHERE table_name = ?", layer); - query_get_layer_geom_meta.next(); - const std::string layer_geometry_column = query_get_layer_geom_meta.get(0); + read_ss << std::endl; + LOG(read_ss.str(), LogLevel::DEBUG); read_ss.str(""); + #endif + + // Get layer feature metadata (geometry column name + type) + auto query_get_layer_geom_meta = db.query("SELECT column_name FROM gpkg_geometry_columns WHERE table_name = ?", layer); + query_get_layer_geom_meta.next(); + const std::string layer_geometry_column = query_get_layer_geom_meta.get(0); + + // Get layer + LOG(LogLevel::DEBUG, "Reading %d features from layer %s.", layer_feature_count, layer.c_str()); + auto query_get_layer = db.query("SELECT * FROM " + layer + joined_ids, sub_ids); + query_get_layer.next(); - // Get layer - auto query_get_layer = db.query("SELECT * FROM " + layer + joined_ids, ids); - query_get_layer.next(); + // build features out of layer query + if (ids.size() == 0) + features.reserve(layer_feature_count); + while(!query_get_layer.done()) { + geojson::Feature feature = build_feature( + query_get_layer, + id_column, + layer_geometry_column + ); + + features.push_back(feature); + query_get_layer.next(); + } - // build features out of layer query - std::vector features; - features.reserve(layer_feature_count); - while(!query_get_layer.done()) { - geojson::Feature feature = build_feature( - query_get_layer, - id_column, - layer_geometry_column - ); - - features.push_back(feature); - query_get_layer.next(); - } + // get layer bounding box from features + // + // GeoPackage contains a bounding box in the SQLite DB, + // however, it is in the SRS of the GPKG. By creating + // the bbox after the features are built, the projection + // is already done. This also should be fairly cheap to do. + for (const auto& feature : features) { + const auto& bbox = feature->get_bounding_box(); + min_x = bbox[0] < min_x ? bbox[0] : min_x; + min_y = bbox[1] < min_y ? bbox[1] : min_y; + max_x = bbox[2] > max_x ? bbox[2] : max_x; + max_y = bbox[3] > max_y ? bbox[3] : max_y; + } - // get layer bounding box from features - // - // GeoPackage contains a bounding box in the SQLite DB, - // however, it is in the SRS of the GPKG. By creating - // the bbox after the features are built, the projection - // is already done. This also should be fairly cheap to do. - double min_x = std::numeric_limits::infinity(); - double min_y = std::numeric_limits::infinity(); - double max_x = -std::numeric_limits::infinity(); - double max_y = -std::numeric_limits::infinity(); - for (const auto& feature : features) { - const auto& bbox = feature->get_bounding_box(); - min_x = bbox[0] < min_x ? bbox[0] : min_x; - min_y = bbox[1] < min_y ? bbox[1] : min_y; - max_x = bbox[2] > max_x ? bbox[2] : max_x; - max_y = bbox[3] > max_y ? bbox[3] : max_y; } auto fc = std::make_shared(