Spatial Filter Push-down and Clustering

WherobotsDB skips reading data files that can’t contain rows matching a spatial filter. Consider the following spatial range query:

SELECT * FROM <table_name> WHERE ST_Contains(ST_GeomFromWKT('POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))'), <geometry_column>);

WherobotsDB skips the data files that have no geometry inside the polygon. This matters most when the table is large and the filter is selective. For example, if the table holds 1 TB of data and the filter selects 1% of it, WherobotsDB ideally reads about 10 GB to answer the query. The examples on this page use the catalog name org_catalog; replace it with the name of the catalog you write to.

How spatial filter push-down works

WherobotsDB collects and records the spatial statistics of data files when writing data to the table. The spatial statistics includes the minimum bounding rectangle (MBR) of the geometries in the data file. When a spatial query is issued, WherobotsDB will first check the MBR of the geometries in the data file against the spatial filter. If the MBR doesn’t intersect with the spatial filter, WherobotsDB will skip reading the data file. The following figure shows how spatial filter push-down works. The MBRs of each data file are kept in the table metadata for pruning data files. The query window shown as the red rectangle only overlaps with the bounding box of datafile-2.parquet and datafile-4.parquet, so WherobotsDB can safely skip scanning datafile-1.parquet and datafile-3.parquet, reducing the IO cost and answer the query faster.

Spatial Filter Push-down and Data Skipping

WherobotsDB keeps these statistics for a geometry(<srid>) column in the standard Apache Iceberg column bounds. For a geometry column with no declared CRS and for a raster column, it keeps them only on format-version 1 and 2 tables, so for spatial filters on those columns it can’t skip data files on a format-version 3 table. See Spatial file skipping.

WherobotsDB skips the most files when rows that are near each other are stored in the same files. Cluster the table with CREATE SPATIAL INDEX, or partition it by a column related to location, as described in Other ways to cluster spatial data.

Raster columns

When WherobotsDB writes data files containing raster columns, it records the WGS84 bounding box of the rasters in each file. For raster spatial filters, such as RS_Intersects, RS_Contains, and RS_Within, it uses these bounds to skip unrelated data files. Push-down works for both in-db and out-db rasters. Raster bounds are recorded only in Apache Iceberg format-version 1 and 2 tables, so WherobotsDB skips data files for raster filters only on those tables. New tables in Wherobots-managed catalogs, such as org_catalog, are format-version 3, so create a raster table that you filter spatially with TBLPROPERTIES ('format-version'='2'), as in Partitioning rasters by SRID. A format-version 2 table can’t have a GEOMETRY(<srid>) or GEOGRAPHY(<srid>) column.

Clustering with CREATE SPATIAL INDEX

WherobotsDB adds the CREATE SPATIAL INDEX statement, which rewrites the table so its rows are sorted by a geometry or raster column:

CREATE SPATIAL INDEX FOR <table_name> USING hilbert(<geometry_column>, <precision> [, <minX>, <minY>, <maxX>, <maxY> ]) [ WHERE <condition> ] [ OPTIONS <options> ];

This statement will rewrite the data files of the table and cluster the data by the geometry column. Currently, WherobotsDB only supports the hilbert index strategy. hilbert sorts the data by the Hilbert index of the specified geometry field. Hilbert index is based on a space filling curve called Hilbert curve, which is very efficient at sorting geospatial data by spatial proximity. You can specify the precision of the hilbert curve using the precision parameter. The precision parameter is the number of bits used to represent the Hilbert index. The higher the precision, the more accurate the sorting will be.

The two-argument form, hilbert(<geometry_column>, <precision>), clusters over World Geodetic System 1984 (WGS84) longitude and latitude bounds, [-180, 180] and [-90, 90]. On a geometry(<srid>) column the declared CRS is checked against those bounds, so the statement is rejected when the column is in another CRS. Pass <minX>, <minY>, <maxX>, <maxY> to cluster a column in another CRS, as described in Clustering a column that is not WGS84.

For example, cluster the data in table org_catalog.test_db.test_table by the geometry column geom:

CREATE SPATIAL INDEX FOR org_catalog.test_db.test_table USING hilbert(geom, 10);

CREATE SPATIAL INDEX FOR is shorthand for calling the rewrite_data_files procedure:

CALL org_catalog.system.rewrite_data_files(
    table => 'org_catalog.test_db.test_table',
    strategy => 'sort',
    sort_order => 'hilbert(geom, 10)',
    options => map('rewrite-all', 'true')
)

The OPTIONS clause passes any option of the rewrite_data_files procedure. For example, set the target file size to 10MB:

CREATE SPATIAL INDEX FOR org_catalog.test_db.test_table USING hilbert(geom, 10)
OPTIONS map('target-file-size-bytes', '10485760')

This is equivalent to

CALL org_catalog.system.rewrite_data_files(
    table => 'org_catalog.test_db.test_table',
    strategy => 'sort',
    sort_order => 'hilbert(geom, 10)',
    options => map('target-file-size-bytes', '10485760', 'rewrite-all', 'true')
)

A WHERE clause limits the rewrite to the data files that may contain matching rows. Every row in those files is rewritten, including rows that don’t match. For example, rewrite the data files that may hold rows with state = 'CA':

CREATE SPATIAL INDEX FOR org_catalog.test_db.test_table USING hilbert(geom, 10)
WHERE state = 'CA'

This is equivalent to

CALL org_catalog.system.rewrite_data_files(
    table => 'org_catalog.test_db.test_table',
    strategy => 'sort',
    sort_order => 'hilbert(geom, 10)',
    where => "state = 'CA'",
    options => map('rewrite-all', 'true')
)

Usually, you don’t need to specify the target-file-count option, as the CREATE SPATIAL INDEX statement will automatically determine the size of output data files based on the size of the table to balance the file size and the number of data files. If the default behavior of CREATE SPATIAL INDEX is not desirable, you can specify the following options to control the number of data files:

target-file-count: The number of output data files. This option is useful when you want to rewrite the table to an expected number of data files. The exact number of data files may be different from the expected number due to various factors such as compression, heuristics for avoiding too large data files, etc.

Clustering a column that is not WGS84

Running the two-argument form on a column whose declared CRS is not WGS84 is rejected, and the message names both ways out:

Cannot build a Hilbert spatial index on geometry column 'geom' with the default
bounds: its CRS is 'EPSG:3857', but the two-argument hilbert(geom, resolution) form
clusters over WGS84 lon/lat bounds ([-180, 180] / [-90, 90]), which do not cover
this CRS's coordinates. Pass explicit bounds -- hilbert(geom, resolution, minX,
minY, maxX, maxY) -- or reproject to WGS84 (OGC:CRS84 / EPSG:4326) with
ST_Transform first.

To cluster in place, pass bounds that cover the column’s coordinates. For a Web Mercator column, that is the square Web Mercator extent, ±20037508.34 meters on both axes, which is longitude ±180 and latitude ±85.0511:

CREATE SPATIAL INDEX FOR org_catalog.test_db.test_table_3857
USING hilbert(geom, 10, -20037508.34, -20037508.34, 20037508.34, 20037508.34);

Bounds have to contain every value in the column. Widen them if your data reaches beyond that square, which EPSG:3857 permits up to latitude ±85.06, or ±20048966.10 meters.

Clustering by a raster column

CREATE SPATIAL INDEX also clusters a table by a raster column. For example, cluster the data in table org_catalog.test_db.test_table by the raster column rast:

CREATE SPATIAL INDEX FOR org_catalog.test_db.test_table USING hilbert(rast, 10);

This is only recommended for the following scenarios:

The original table does not have a good column for partitioning.
The table contains out-db rasters or very small in-db rasters.

Other ways to cluster spatial data

Besides CREATE SPATIAL INDEX, you can cluster spatial data by partitioning or sorting the table. If the data you are working with has at least one column that is closely related to the spatial property, it is usually a good idea to partition the table using that column. For example, if you are working with a table containing census data, and the data has a state column indicating the state ID of the data, you can partition the table by state:

CREATE TABLE org_catalog.test_db.test_table (
  id bigint,
  state string,
  geom geometry(4326)
) PARTITIONED BY (state)

The state partition works well for both regular queries with filters such as state = 'some_state' and spatial queries with filters such as ST_Intersects(geom, <query_window>). For the former, WherobotsDB skips data files whose state partition doesn’t match the filter. For the latter, it skips data files using their spatial statistics. The same applies to raster tables. For a table of rasters of different cities, partitioning by the city name column stores the rasters of each city in the same partition, so WherobotsDB skips more data files for raster filters.

Partitioning by spatial proximity

If your dataset does not have a field to directly partition by, you can partition the dataset by the S2 cell ID of the geometry. This requires selecting a proper S2 resolution level suited for your dataset.

CREATE TABLE org_catalog.test_db.test_table (
  id bigint,
  s2 bigint,
  state string,
  geom geometry(4326)
) PARTITIONED BY (s2)

The data inserted into the table should have the s2 column set to the S2 cell ID of the geometry, and the values in geom have to carry the SRID of the column, as described in Matching the column CRS:

Python
Scala
Java

df.withColumn("s2", expr("array_max(ST_S2CellIds(ST_Centroid(geom), 4))"))\
    .writeTo("org_catalog.test_db.test_table")\
    .append()

val dfWithS2 = df.withColumn("s2", expr("array_max(ST_S2CellIds(ST_Centroid(geom), 4))"))
dfWithS2.writeTo("org_catalog.test_db.test_table").append()

Dataset<Row> dfWithS2 = df.withColumn("s2", expr("array_max(ST_S2CellIds(ST_Centroid(geom), 4))"));
dfWithS2.writeTo("org_catalog.test_db.test_table").append();

Though the table was partitioned using S2, we don’t need to add a predicate for the s2 partitioning column to achieve data skipping. WherobotsDB will use the spatial statistics of data files and automatically skip reading data files that don’t contain data that satisfy the spatial filter.

SELECT * FROM org_catalog.test_db.test_table WHERE ST_Intersects(geom, ST_GeomFromText('...'))
-- No need to add a predicate for the s2 partitioning column

Partitioning rasters by SRID

When the rasters in a table are in Universal Transverse Mercator (UTM) coordinate reference systems spread over several UTM zones, partition the table by SRID. This example partitions the EuroSAT dataset by SRID, in a format-version 2 table so that raster bounds are recorded:

Python
Scala
Java

# Create a format-version 2 table partitioned by SRID, so raster bounds are recorded
sedona.sql("""
CREATE TABLE org_catalog.test_db.eurosat_ms_srid (
    path STRING,
    length BIGINT,
    rast RASTER,
    srid INT
) USING iceberg
PARTITIONED BY (srid)
TBLPROPERTIES ('format-version'='2')
""")

# Load rasters from GeoTIFF files
df_binary = sedona.read.format("binaryFile")\
    .option("pathGlobFilter", "*.tif")\
    .option("recursiveFileLookup", "true")\
    .load('s3://wherobots-examples/data/eurosat_small')
df_geotiff = df_binary.withColumn("rast", expr("RS_FromGeoTiff(content)"))\
    .select("path", "length", "rast")

# Append the rasters to the table, sorted by SRID
df_geotiff.withColumn("srid", expr("RS_SRID(rast)"))\
    .sort("srid")\
    .writeTo("org_catalog.test_db.eurosat_ms_srid").append()

// Create a format-version 2 table partitioned by SRID, so raster bounds are recorded
sedona.sql("""
CREATE TABLE org_catalog.test_db.eurosat_ms_srid (
    path STRING,
    length BIGINT,
    rast RASTER,
    srid INT
) USING iceberg
PARTITIONED BY (srid)
TBLPROPERTIES ('format-version'='2')
""")

// Load rasters from GeoTIFF files
val df_binary = sedona.read.format("binaryFile")
    .option("pathGlobFilter", "*.tif")
    .option("recursiveFileLookup", "true")
    .load("s3://wherobots-examples/data/eurosat_small")

val df_geotiff = df_binary.withColumn("rast", expr("RS_FromGeoTiff(content)"))
    .select("path", "length", "rast")

// Append the rasters to the table, sorted by SRID
df_geotiff.withColumn("srid", expr("RS_SRID(rast)"))
    .sort("srid")
    .writeTo("org_catalog.test_db.eurosat_ms_srid").append()

// Create a format-version 2 table partitioned by SRID, so raster bounds are recorded
sedona.sql("CREATE TABLE org_catalog.test_db.eurosat_ms_srid ("
    + "path STRING, length BIGINT, rast RASTER, srid INT) USING iceberg "
    + "PARTITIONED BY (srid) TBLPROPERTIES ('format-version'='2')");

// Load rasters from GeoTIFF files
Dataset<Row> df_binary = sedona.read().format("binaryFile")
    .option("pathGlobFilter", "*.tif")
    .option("recursiveFileLookup", "true")
    .load("s3://wherobots-examples/data/eurosat_small");

Dataset<Row> df_geotiff = df_binary.withColumn("rast", expr("RS_FromGeoTiff(content)"))
    .select("path", "length", "rast");

// Append the rasters to the table, sorted by SRID
df_geotiff.withColumn("srid", expr("RS_SRID(rast)"))
    .sort("srid")
    .writeTo("org_catalog.test_db.eurosat_ms_srid").append();

Subsequent queries such as

SELECT * FROM org_catalog.test_db.eurosat_ms_srid WHERE RS_Intersects(rast, ...)

will be much faster than the same query on a non-partitioned table.

Sorting by GeoHash

You can sort the rows by ST_GeoHash(<geometry_column>, <precision>) as you write them, to cluster the data by spatial proximity. Create the destination table first, with a geometry(<srid>) column:

CREATE TABLE org_catalog.test_db.test_table_sorted (
  id bigint,
  state string,
  geom geometry(4326)
)

Then sort the rows and append them. The example expects df to hold WGS84 geometries with SRID 4326; for geometries in another CRS, or with SRID 0, see Matching the column CRS:

Python
Scala
Java

sortedDf = df.withColumn("geohash", expr("ST_GeoHash(ST_Centroid(geom), 20)"))\
    .sort(col("geohash"))\
    .drop("geohash")
sortedDf.writeTo("org_catalog.test_db.test_table_sorted").append()

val sortedDf = df.withColumn("geohash", expr("ST_GeoHash(ST_Centroid(geom), 20)"))
    .sort(col("geohash"))
    .drop("geohash")

sortedDf.writeTo("org_catalog.test_db.test_table_sorted").append()

Dataset<Row> sortedDf = df.withColumn("geohash", expr("ST_GeoHash(ST_Centroid(geom), 20)"))
    .sort(col("geohash"))
    .drop("geohash");

sortedDf.writeTo("org_catalog.test_db.test_table_sorted").append();

It is generally recommended to use CREATE SPATIAL INDEX FOR for sorting the table over manually sorting the table. CREATE SPATIAL INDEX FOR will use the Hilbert curve to sort the table, which is more efficient than sorting by GeoHash.

Iceberg Tables Overview

Spatial file skipping on each format version.

Geometry and Geography Tables

Creating, writing, and querying geometry tables.

Raster Tables

Creating, writing, and querying raster tables.

Upgrade Existing Geometry Tables

Copy a table with no declared CRS into one whose geometry column keeps file skipping on format-version 3.

Start here

Havasu data catalogs

Vector data

Raster data

Iceberg tables

Maps and visualization

Spatial Filter Push-down and Clustering

How spatial filter push-down works

Raster columns

Clustering with CREATE SPATIAL INDEX

Clustering a column that is not WGS84

Clustering by a raster column

Other ways to cluster spatial data

Partitioning by spatial proximity

Partitioning rasters by SRID

Sorting by GeoHash

Iceberg Tables Overview

Geometry and Geography Tables

Raster Tables

Upgrade Existing Geometry Tables

How spatial filter push-down works

Raster columns

Clustering with CREATE SPATIAL INDEX

Clustering a column that is not WGS84

Clustering by a raster column

Other ways to cluster spatial data

Partitioning by columns closely related to spatial property

Partitioning by spatial proximity

Partitioning rasters by SRID

Sorting by GeoHash

Related pages

Iceberg Tables Overview

Geometry and Geography Tables

Raster Tables

Upgrade Existing Geometry Tables