Mejora de la velocidad de consulta espacial de SQLite para gpkg multicapa y Python

Question

Tengo un GeoPackage con dos capas de polígonos (WGS84). La primera capa es una malla de azulejos Sentinel 2 con unos 800 polígonos. La otra capa es de unos 300 000 polígonos pequeños (AOIs). Cada capa tiene un índice espacial rtree.

Tablas en el gpkg:

[('gpkg_spatial_ref_sys',), ('gpkg_contents',), ('gpkg_ogr_contents',),
 ('gpkg_geometry_columns',), ('gpkg_tile_matrix_set',), ('gpkg_tile_matrix',),
 ('aois',), ('gpkg_extensions',), ('rtree_aois_geom',),
 ('rtree_aois_geom_rowid',), ('rtree_aois_geom_node',),
 ('rtree_aois_geom_parent',), ('sentinel2_tiles',),
 ('rtree_sentinel2_tiles_geom',), ('rtree_sentinel2_tiles_geom_rowid',),
 ('rtree_sentinel2_tiles_geom_node',), ('rtree_sentinel2_tiles_geom_parent',)]

Me gustaría comprobar qué AOIs están dentro de qué mosaico de la cuadrícula de Sentinel 2 (sé de los solapamientos entre los mosaicos de Sentinel 2, pero no importa en este caso mientras cada AOI tenga asignado un ID de mosaico). El problema al que me enfrento es que la consulta tarda mucho tiempo en ejecutarse. Hacer consultas de geometría con una sola capa es bastante rápido (2 segundos o menos en promedio).

Actualmente esta es la consulta:

SELECT a.fid, b.grid
FROM aois a, sentinel2_tiles b
WHERE ST_Within(ST_envelope(a.geom), ST_envelope(b.geom));

¿Cómo puedo aumentar/optimizar la consulta? Alternativamente, ¿sería mejor/ más rápido utilizar la función OGR Within()?

Answer 1

Por el momento, obtener primero los nombres únicos de las baldosas de la cuadrícula y luego multiprocesar una consulta interna para cada baldosa de la cuadrícula parece funcionar. Se ha reducido el tiempo de consulta a menos de un minuto. La primera consulta que devuelve los IDs de las baldosas de la cuadrícula (lista llamada results ):

"""SELECT DISTINCT b.grid 
   FROM sentinel2_tiles b;"""

Segunda consulta (lista de consultas que luego se multiprocesan):

query_list = [f"""SELECT a.fid, b.grid 
                  FROM aois a,
                       sentinel2_tiles b
                  WHERE b.name = '{tile}' AND
                        ST_Within(ST_envelope(a.geom),
                        ST_envelope(b.geom));"""
              for tile in results]