features: use the GeoPackage spatial index for spatial predicates

xtraplatform-features-sql/src/main/java/de/ii/xtraplatform/features/sql/domain/FeatureProviderSql.java

@@ -114,6 +114,7 @@
 import de.ii.xtraplatform.streams.domain.Reactive.Stream;
 import de.ii.xtraplatform.streams.domain.Reactive.Transformer;
 import de.ii.xtraplatform.values.domain.ValueStore;
+import java.sql.SQLException;
 import java.time.ZoneId;
 import java.util.AbstractMap.SimpleImmutableEntry;
 import java.util.ArrayList;
@@ -613,7 +614,8 @@ protected boolean onStartup() throws InterruptedException {
             type ->
                 Optional.ofNullable(queryMappings.get(type))
                     .filter(mappings -> mappings.size() == 1)
-                    .map(mappings -> mappings.get(0)));
+                    .map(mappings -> mappings.get(0)),
+            getSpatialIndexes());
     AggregateStatsQueryGenerator queryGeneratorSql =
         new AggregateStatsQueryGenerator(sqlDialect, filterEncoder);
 
@@ -667,6 +669,24 @@ protected boolean onStartup() throws InterruptedException {
     return true;
   }
 
+  /**
+   * The geometry columns with a spatial index that the query generator has to name explicitly to
+   * make use of. Without it, spatial queries keep working, only without the index, so a failure to
+   * determine it is logged and does not keep the provider from starting.
+   */
+  private Map<String, String> getSpatialIndexes() {
+    try {
+      return getSqlClient().getSpatialIndexes();
+    } catch (SQLException | RuntimeException e) {
+      LogContext.errorAsWarn(
+          LOGGER,
+          e,
+          "Could not determine the spatial indexes of the feature provider with id '{}'. Spatial queries will not use a spatial index.",
+          getId());
+      return Map.of();
+    }
+  }
+
   @Override
   protected void onStarted() {
     changes()

xtraplatform-features-sql/src/main/java/de/ii/xtraplatform/features/sql/domain/SqlClientBasic.java

@@ -27,4 +27,8 @@ default DbInfo getDbInfo() throws SQLException {
   default Map<String, GeoInfo> getGeoInfo() throws SQLException {
     return getDbmsAdapter().getGeoInfo(getConnection(), getDbInfo());
   }
+
+  default Map<String, String> getSpatialIndexes() throws SQLException {
+    return getDbmsAdapter().getSpatialIndexes(this);
+  }
 }

xtraplatform-features-sql/src/main/java/de/ii/xtraplatform/features/sql/domain/SqlDbmsAdapter.java

@@ -40,6 +40,20 @@ public interface SqlDbmsAdapter {
 
   Map<String, GeoInfo> getGeoInfo(Connection connection, DbInfo dbInfo) throws SQLException;
 
+  /**
+   * Returns the geometry columns that have a spatial index which the query generator has to name
+   * explicitly to make use of, keyed by {@code table.column} in lower case (SQL identifiers are
+   * compared case-insensitively here, as SQLite does for ASCII). The value is the column of the
+   * table that the index entries are keyed on, which a query has to join the index on.
+   *
+   * <p>The default is an empty map, which is correct for every DBMS whose spatial operators consult
+   * the spatial index by themselves. It takes the client rather than a connection so that those
+   * adapters do not have to be handed one they will not use.
+   */
+  default Map<String, String> getSpatialIndexes(SqlClientBasic sqlClient) throws SQLException {
+    return Map.of();
+  }
+
   DbInfo getDbInfo(Connection connection) throws SQLException;
 
   Collator getRowSortingCollator(Optional<String> defaultCollation);

xtraplatform-features-sql/src/main/java/de/ii/xtraplatform/features/sql/domain/SqlDialect.java

@@ -138,6 +138,28 @@ default String getSpatialOperatorMatch(SpatialFunction spatialFunction) {
     return "";
   }
 
+  /**
+   * Returns a predicate that restricts the rows to those whose geometry bounding box intersects the
+   * bounding box given as {@code min}/{@code max}, evaluated with a spatial index on {@code
+   * table}.{@code column}, or an empty optional if the dialect cannot use a spatial index that way.
+   *
+   * <p>The predicate is a conjunct that is added to the exact spatial predicate, never a
+   * replacement for it. Dialects whose spatial operators already consult the spatial index by
+   * themselves (such as PostGIS, where {@code ST_Intersects} includes the {@code &&} bounding box
+   * operator) do not implement this.
+   *
+   * @param table the table holding the geometry column
+   * @param column the geometry column
+   * @param alias the alias the table is referenced by in the query
+   * @param keyColumn the column of the table that the index entries are keyed on
+   * @param min the lower corner of the bounding box, in the native CRS
+   * @param max the upper corner of the bounding box, in the native CRS
+   */
+  default Optional<String> getSpatialIndexPredicate(
+      String table, String column, String alias, String keyColumn, double[] min, double[] max) {
+    return Optional.empty();
+  }
+
   default String getTemporalOperator(TemporalFunction temporalFunction) {
     // this is implementation specific
     return null;
@@ -166,4 +188,23 @@ default String applyToExpression(
 
   Map<SpatialFunction, String> SPATIAL_OPERATORS_3D =
       new ImmutableMap.Builder<SpatialFunction, String>().build();
+
+  /**
+   * The spatial operators for which a match implies that the bounding boxes of both operands
+   * intersect. Only for these may a bounding box predicate from a spatial index be added as a
+   * conjunct: the conjunct is then implied by the exact predicate, so it changes neither the result
+   * nor the meaning of the predicate under negation.
+   *
+   * <p>{@code S_DISJOINT} is absent because two geometries can be disjoint while being arbitrarily
+   * far apart, so a bounding box predicate would discard matching rows.
+   */
+  Set<SpatialFunction> SPATIAL_OPERATORS_IMPLYING_BBOX_INTERSECTION =
+      ImmutableSet.of(
+          SpatialFunction.S_EQUALS,
+          SpatialFunction.S_TOUCHES,
+          SpatialFunction.S_WITHIN,
+          SpatialFunction.S_OVERLAPS,
+          SpatialFunction.S_CROSSES,
+          SpatialFunction.S_INTERSECTS,
+          SpatialFunction.S_CONTAINS);
 }

xtraplatform-features-sql/src/main/java/de/ii/xtraplatform/features/sql/domain/SqlDialectGpkg.java

@@ -206,6 +206,43 @@ public String escapeString(String value) {
     return value.replaceAll("'", "''");
   }
 
+  /**
+   * GeoPackage keeps the bounding boxes of a geometry column in an R-Tree, a virtual table named
+   * {@code rtree_<table>_<column>}. SQLite cannot infer from a spatial operator that such a table
+   * is relevant, so unless the query names it, the operator is evaluated for every row of the
+   * table. The predicate returned here names it, which turns the table scan into an index lookup.
+   *
+   * <p>The R-Tree stores the bounding boxes as 32-bit floats, rounded outwards, so the lookup never
+   * discards a geometry that the exact predicate would have matched.
+   */
+  @Override
+  public Optional<String> getSpatialIndexPredicate(
+      String table, String column, String alias, String keyColumn, double[] min, double[] max) {
+    if (min.length < 2 || max.length < 2) {
+      return Optional.empty();
+    }
+
+    return Optional.of(
+        String.format(
+            "%s.%s IN (SELECT id FROM %s WHERE maxx >= %s AND minx <= %s AND maxy >= %s AND miny <= %s)",
+            alias,
+            quoteIdentifier(keyColumn),
+            quoteIdentifier(String.format("rtree_%s_%s", table, column)),
+            toNumericLiteral(min[0]),
+            toNumericLiteral(max[0]),
+            toNumericLiteral(min[1]),
+            toNumericLiteral(max[1])));
+  }
+
+  private static String quoteIdentifier(String identifier) {
+    return String.format("\"%s\"", identifier.replace("\"", "\"\""));
+  }
+
+  private static String toNumericLiteral(double value) {
+    // Double.toString is independent of the default locale, unlike String.format("%f", ...)
+    return Double.toString(value);
+  }
+
   @Override
   public String applyToExpression(
       String table, String name, Map<String, String> subDecoderPaths, boolean spatial) {

xtraplatform-features-sql/src/main/java/de/ii/xtraplatform/features/sql/infra/db/SqlDbmsAdapterGpkg.java

@@ -9,12 +9,14 @@
 
 import com.github.azahnen.dagger.annotations.AutoBind;
 import com.google.common.collect.ImmutableList;
+import com.google.common.collect.ImmutableMap;
 import de.ii.xtraplatform.base.domain.AppContext;
 import de.ii.xtraplatform.blobs.domain.ResourceStore;
 import de.ii.xtraplatform.crs.domain.EpsgCrs;
 import de.ii.xtraplatform.crs.domain.EpsgCrs.Force;
 import de.ii.xtraplatform.features.sql.domain.ConnectionInfoSql;
 import de.ii.xtraplatform.features.sql.domain.ImmutableGeoInfo;
+import de.ii.xtraplatform.features.sql.domain.SqlClientBasic;
 import de.ii.xtraplatform.features.sql.domain.SqlDbmsAdapter;
 import de.ii.xtraplatform.features.sql.domain.SqlDialect;
 import de.ii.xtraplatform.features.sql.domain.SqlDialectGpkg;
@@ -172,6 +174,46 @@ public List<String> getSystemTables() {
         "SpatialIndex");
   }
 
+  /**
+   * Reports the geometry columns that have a GeoPackage R-Tree, mapped to the column that the
+   * R-Tree entries are keyed on.
+   *
+   * <p>The R-Tree virtual table itself is looked up instead of the {@code gpkg_rtree_index} row in
+   * {@code gpkg_extensions}, because it is the table that the generated query has to name, and the
+   * two can disagree: the extension may be registered for a column whose R-Tree was dropped.
+   *
+   * <p>The key column is the primary key of the feature table, which is what the R-Tree triggers
+   * store. It is reported instead of assuming {@code rowid}: a GeoPackage feature table is required
+   * to have an {@code INTEGER PRIMARY KEY}, which SQLite makes an alias of {@code rowid}, but in a
+   * file that does not follow that requirement the two differ and joining on {@code rowid} would
+   * match the wrong rows. Feature tables with a composite primary key are left out, as their rows
+   * cannot be addressed by the single id column of an R-Tree.
+   */
+  @Override
+  public Map<String, String> getSpatialIndexes(SqlClientBasic sqlClient) throws SQLException {
+    String query =
+        "SELECT gc.table_name, gc.column_name, ti.name FROM gpkg_geometry_columns gc"
+            + " JOIN sqlite_master m ON m.type = 'table' AND lower(m.name) = lower('rtree_' || gc.table_name || '_' || gc.column_name)"
+            + " JOIN pragma_table_info(gc.table_name) ti ON ti.pk = 1"
+            + " GROUP BY gc.table_name, gc.column_name HAVING count(*) = 1;";
+
+    ImmutableMap.Builder<String, String> spatialIndexes = ImmutableMap.builder();
+
+    // the connection comes from the pool, closing it hands it back
+    try (Connection connection = sqlClient.getConnection();
+        Statement statement = connection.createStatement();
+        ResultSet resultSet = statement.executeQuery(query)) {
+      while (resultSet.next()) {
+        spatialIndexes.put(
+            String.format("%s.%s", resultSet.getString(1), resultSet.getString(2))
+                .toLowerCase(Locale.ROOT),
+            resultSet.getString(3));
+      }
+    }
+
+    return spatialIndexes.build();
+  }
+
   @Override
   public Map<String, GeoInfo> getGeoInfo(Connection connection, DbInfo dbInfo) throws SQLException {
     if (!(dbInfo instanceof DbInfoGpkg)