¿Cómo configurar correctamente los índices para las consultas de distancia PostGIS?

Estoy creando una aplicación que se supone que consulta y devuelve cada Recorduna en una tabla que está a Xkilómetros de distancia PointX. Recordsy PointXlas posiciones se determinan a partir de la (long/lat)información proporcionada por Google Geocode API.

Soy nuevo en PostGIS. Después de una investigación rápida, encontré esta pregunta . La respuesta parece estar en la línea de:

SELECT *
FROM your_table
WHERE ST_Distance_Sphere(the_geom, ST_MakePoint(your_lon,your_lat)) <= radius_mi * 1609.34

El problema es: aunque solo estoy comenzando en SIG, cuando miro la consulta anterior, no puedo imaginar cómo esto puede usar un índice. Hay 2 llamadas a funciones. Me imagino que la mesa está siendo escaneada para cada uno Record. Quiero estar equivocado :)

Pregunta: ¿PostGIS tiene algún tipo de índice capaz de hacer que la consulta anterior sea eficaz? Si no, ¿cuál sería el enfoque recomendado para hacer lo que necesito?

Hay dos claves para obtener un buen rendimiento de consultas geodésicas con tablas grandes con geometrycolumnas utilizando datos geográficos WGS 1984 (SRID 4326):

1. Use la ST_DWithinfunción, que busca utilizando un índice espacial disponible, y encontrará características de geografía con una distancia cartesiana
2. Cree un índice adicional en el reparto de geografía, por lo que ST_DWithinpuede usarlo

Así que echemos un vistazo a lo que sucede en el mundo real. Primero necesitamos crear y llenar una tabla de un millón de puntos aleatorios:

DROP TABLE IF EXISTS example1
;

CREATE TABLE example1 (
    idcol   serial      NOT NULL,
    geomcol geometry        NULL,
    CONSTRAINT  example1_pk PRIMARY KEY (idcol),
    CONSTRAINT  enforce_srid CHECK (st_srid(geomcol) = 4326)
)
with (
    OIDS=FALSE
);

INSERT INTO example1(geomcol)
SELECT  ST_SetSRID(
            ST_MakePoint(
            (random()*360.0) - 180.0,
            (acos(1.0 - 2.0 * random()) * 2.0 - pi()) * 90.0 / pi()),
            4326) as geomcol
FROM  generate_series(1, 1000000) vtab;

CREATE INDEX example1_spx ON example1 USING GIST (geomcol);
-- (took about 22 sec)

Si ejecutamos la consulta ST_Distance, obtenemos el escaneo completo esperado de la tabla:

EXPLAIN ANALYZE VERBOSE
SELECT  count(*)
FROM    example1
WHERE   ST_Distance(geomcol::geography,ST_SetSRID(ST_MakePoint(6.9333,46.8167),4326)::geography) < 30 * 1609.34
;

Aggregate  (cost=274167.33..274167.34 rows=1 width=0) (actual time=4940.531..4940.532 rows=1 loops=1)
  Output: count(*)
  ->  Seq Scan on bob.example1  (cost=0.00..273334.00 rows=333333 width=0) (actual time=592.766..4940.509 rows=11 loops=1)
        Output: idcol, geomcol
        Filter: (_st_distance((example1.geomcol)::geography, '0101000020E61000005D6DC5FEB2BB1B40545227A089684740'::geography, 0::double precision, true) < 48280.2::double precision)
        Rows Removed by Filter: 999989
Planning time: 2.137 ms
Execution time: 4940.568 ms

Ahora, si usamos ST_DWithin, todavía obtenemos un escaneo completo de la tabla (aunque más rápido):

EXPLAIN ANALYZE VERBOSE
SELECT  count(*)
FROM    example1
WHERE   ST_DWithin(geomcol::geography,ST_SetSRID(ST_MakePoint(6.9333,46.8167),4326)::geography,30 * 1609.34)
;

Aggregate  (cost=405867.33..405867.34 rows=1 width=0) (actual time=908.716..908.716 rows=1 loops=1)
  Output: count(*)
  ->  Seq Scan on bob.example1  (cost=0.00..405834.00 rows=13333 width=0) (actual time=38.449..908.700 rows=7 loops=1)
        Output: idcol, geomcol
        Filter: (((example1.geomcol)::geography && '0101000020E61000005D6DC5FEB2BB1B40545227A089684740'::geography) AND ('0101000020E61000005D6DC5FEB2BB1B40545227A089684740'::geography && _st_expand((example1.geomcol)::geography, 48280.2::double precision) (...)
        Rows Removed by Filter: 999993
Planning time: 2.017 ms
Execution time: 908.763 ms

Y esta es la última pieza: construir el índice de cobertura (geografía de reparto):

CREATE INDEX example1_gpx ON example1 USING GIST (geography(geomcol));
-- (Takes an extra 13 sec)

EXPLAIN ANALYZE VERBOSE
SELECT  count(*)
FROM    example1
WHERE   ST_DWithin(geomcol::geography,ST_SetSRID(ST_MakePoint(6.9333,46.8167),4326)::geography,30 * 1609.34)
;

Aggregate  (cost=96538.95..96538.96 rows=1 width=0) (actual time=0.775..0.775 rows=1 loops=1)
  Output: count(*)
  ->  Bitmap Heap Scan on bob.example1  (cost=8671.62..96505.62 rows=13333 width=0) (actual time=0.586..0.769 rows=19 loops=1)
        Output: idcol, geomcol
        Recheck Cond: ((example1.geomcol)::geography && '0101000020E61000005D6DC5FEB2BB1B40545227A089684740'::geography)
        Filter: (('0101000020E61000005D6DC5FEB2BB1B40545227A089684740'::geography && _st_expand((example1.geomcol)::geography, 48280.2::double precision)) AND _st_dwithin((example1.geomcol)::geography, '0101000020E61000005D6DC5FEB2BB1B40545227A089684740':: (...)
        Rows Removed by Filter: 14
        Heap Blocks: exact=33
        ->  Bitmap Index Scan on example1_gpx  (cost=0.00..8668.29 rows=200000 width=0) (actual time=0.384..0.384 rows=33 loops=1)
              Index Cond: ((example1.geomcol)::geography && '0101000020E61000005D6DC5FEB2BB1B40545227A089684740'::geography)
Planning time: 2.572 ms
Execution time: 0.820 ms

Finalmente, el optimizador está usando el índice espacial, y se muestra, pero ¿cuáles son los tres órdenes de magnitud entre amigos?

Algunas advertencias:

• Soy un nerd de la base de datos, por lo que la PC de mi casa tiene 16 Gb de RAM, seis núcleos de 3.3 Ghz y un SSD de 256 Gb para el espacio de tabla predeterminado de la base de datos; Su experiencia puede ser diferente

• Volví a ejecutar el SQL de creación antes de cada consulta, para nivelar el campo de juego con respecto a las páginas "activas" en el caché, pero esto podría producir resultados ligeramente diferentes porque la misma semilla aleatoria no se usó para diferentes ejecuciones

Y una nota:

• Ajusté el rango de latitud original {-90, + 90} para usar el arco-coseno para una distribución de área igual (menos sesgada hacia los polos)