Adquirir una velocidad similar a la de ArcGIS en Postgis

Question

Llevo 3/4 de año utilizando Postgis 2.0 y, aunque me gusta mucho utilizarlo, el excesivo tiempo de procesamiento de las consultas lo ha hecho básicamente inutilizable para mi caso de uso.

Suelo realizar un intenso geoprocesamiento de conjuntos de datos municipales que a menudo tienen cientos de miles de multipolígonos. Estos multipolígonos tienen a veces formas muy irregulares y pueden variar desde 4 puntos hasta 78.000 puntos por multipolígono.

Por ejemplo, cuando intersecciono un conjunto de datos de parcelas con 329.152 multipolígonos con un conjunto de datos de jurisdicción que contiene 525 multipolígonos, obtengo las siguientes estadísticas de tiempo total consumido:

ArcGIS 10.0 (on same host with windows 7 OS): 3 minutes
Postgis:56 minutes (not including geometry pre-processing queries)

En otras palabras, se requiere un 1500% más de tiempo para hacer esta intersección en Postgis que en ArcGIS --¡Y esta es una de mis consultas más sencillas!

Una de las razones por las que ArcGIS supuestamente funciona más rápido se debe a la mejora de los índices. Algunos programadores recientemente descubrí cómo funcionan estos índices y me pregunto si alguien sabe cómo construir estos índices en Postgis (o construir tablas que imiten los índices). Tal vez esto resolvería la mayoría de los problemas de velocidad en Postgis. ¡Sólo puedo esperar que debe haber alguna manera, sobre todo porque ArcGIS sólo puede utilizar 4 GB de RAM, mientras que podría utilizar hasta 4 veces que para mi servidor Postgis!

Por supuesto, hay muchas razones por las que Postgis puede funcionar lentamente, así que voy a proporcionar una versión detallada de las especificaciones de mi sistema:

Machine: Dell XPS 8300 
Processor: i7-2600 CPU @ 3.40 GHz 3.40 GHz 
Memory: Total Memory 16.0 GB (10.0 GB on virtual machine)

Platform: Ubuntu Server 12.04 Virtual Box VM

Potgres Version: 9.1.4
Postgis Version: POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER

También detallo la todo el proceso de instalación que utilicé para configurar postgis incluyendo la creación de la propia VM .

También aumenté la memoria compartida de los 24MB por defecto a 6 GB en el archivo conf y ejecuté los siguientes comandos para permitir la ejecución de postgres:

sudo sysctl -w kernel.shmmax=7516192768 (I know this setting is deleted every time you restart the OS)
sudo /etc/init.d/postgresql restart

Por lo que puedo decir esto no hace absolutamente nada notable en términos de rendimiento.

Aquí están los enlaces a los datos que he utilizado para esta prueba:

1. Parcelas: tcad_parcels_06142012.shp.zip de Ciudad de Austin, TX
2. Jurisdicciones: Límites jurisdiccionales de Ciudad de Austin, TX

Estos son los pasos que he seguido para procesar los datos:

ArcGIS

1. Añadir conjuntos de datos a ArcMap
2. Establecer el sistema de coordenadas a los pies del centro de Texas (srid 2277)
3. Utilice la herramienta de intersección del menú desplegable

Postgis

Importar paquetes con:

shp2pgsql -c -s 2277 -D -i -I -W UTF-8 "tcad_parcels_06142012.shp" "public"."tcad_parcels_06142012" |psql -d postgis_testing -U postgres -h local_ip -p 5432

Importar jurisdicciones utilizando:

shp2pgsql -c -s 2277 -D -i -I -W UTF-8 "jurisdictions.shp" "public"."jurisdictions" |psql -d postgis_testing -U postgres -h local_ip -p 5432

Limpiar la geometría no válida en las parcelas:

DROP TABLE IF EXISTS valid_parcels;
CREATE TABLE valid_parcels(
  gid serial PRIMARY KEY,
  orig_gid integer,
  geom geometry(multipolygon,2277)
);
CREATE INDEX ON valid_parcels USING gist (geom);
INSERT INTO valid_parcels(orig_gid,geom)
  SELECT 
    gid 
    orig_gid,
    st_multi(st_makevalid(geom)) 
  FROM 
    tcad_parcels_06142012;
CLUSTER valid_parcels USING valid_parcels_geom_idx;

Limpiar la geometría no válida en las jurisdicciones:

DROP TABLE IF EXISTS valid_jurisdictions;
CREATE TABLE valid_jurisdictions(
  gid serial PRIMARY KEY,
  orig_gid integer,
  geom geometry(multipolygon,2277)
);
CREATE INDEX ON valid_jurisdictions USING gist (geom);
INSERT INTO valid_jurisdictions(orig_gid,geom)
  SELECT 
    gid 
    orig_gid,
    st_multi(st_makevalid(geom)) 
  FROM 
    jurisdictions;
CLUSTER valid_jurisdictions USING valid_jurisdictions_geom_idx;

Ejecuta el clúster:

cluster;

Analiza el vacío:

vacuum analyze;

Realizar la intersección en las mesas limpias:

CREATE TABLE parcel_jurisdictions(
  gid serial primary key,
  parcel_gid integer,
  jurisdiction_gid integer,
  isect_geom geometry(multipolygon,2277)
);
CREATE INDEX ON parcel_jurisdictions using gist (isect_geom);

INSERT INTO parcel_jurisdictions(parcel_gid,jurisdiction_gid,isect_geom)
  SELECT
    a.orig_gid parcel_gid,
    b.orig_gid jurisdiction_gid,
    st_multi(st_intersection(a.geom,b.geom))
  FROM
    valid_parcels a, valid_jurisdictions b
  WHERE
    st_intersects(a.geom,b.geom);

Explicar la consulta de intersección de análisis:

Total runtime: 3446860.731 ms
        Index Cond: (geom && b.geom)
  ->  Index Scan using valid_parcels_geom_idx on valid_parcels a  (cost=0.00..11.66 rows=2 width=1592) (actual time=0.030..4.596 rows=1366 loops=525)
  ->  Seq Scan on valid_jurisdictions b  (cost=0.00..113.25 rows=525 width=22621) (actual time=0.009..0.755 rows=525 loops=1)
Nested Loop  (cost=0.00..61428.74 rows=217501 width=24213) (actual time=2.625..3445946.889 rows=329152 loops=1)
  Join Filter: _st_intersects(a.geom, b.geom)

Por todo lo que he leído, mi consulta de intersección es eficiente y no tengo ni idea de qué estoy haciendo mal para que la consulta tarde 56 minutos en una geometría limpia.

Answer 1

Un enfoque diferente. Sabiendo que el dolor está en ST_Intersection, y que las pruebas verdadero/falso son rápidas, intentar minimizar la cantidad de geometría que pasa por la intersección podría acelerar las cosas. Por ejemplo, las parcelas que están totalmente contenidas en una jurisdicción no necesitan ser recortadas, pero ST_Intersection probablemente se tomará la molestia de construir parte de la superposición de la intersección antes de darse cuenta de que no tiene que generar ninguna geometría nueva. Así que esto

INSERT INTO parcel_jurisdictions(parcel_gid,jurisdiction_gid,isect_geom)
SELECT
  a.orig_gid AS parcel_gid,
  b.orig_gid AS jurisdiction_gid,

  st_multi(st_intersection(a.geom,b.geom)) AS geom
FROM
  valid_parcels a, valid_jurisdictions b
WHERE
  st_intersects(a.geom, b.geom) and not st_within(a.geom, b.geom)
UNION ALL
SELECT
  a.orig_gid AS parcel_gid,
  b.orig_gid AS jurisdiction_gid,
  a.geom AS geom
FROM
  valid_parcels a, valid_jurisdictions b
WHERE
  st_within(a.geom, b.geom);

O incluso terser

INSERT INTO parcel_jurisdictions(parcel_gid,jurisdiction_gid,isect_geom)
SELECT
  a.orig_gid AS parcel_gid,
  b.orig_gid AS jurisdiction_gid,
  CASE 
     WHEN ST_Within(a.geom,b.geom) 
     THEN a.geom
     ELSE ST_Multi(ST_Intersection(a.geom,b.geom)) 
  END AS geom
FROM valid_parcels a
JOIN valid_jurisdictions b
ON ST_Intersects(a.geom, b.geom)

Incluso podría ser más rápido sin la UNIÓN.

Answer 2

¿Qué pasaría si se omite el "st_multi(st_intersection(a.geom,b.geom))" ¿parte?

¿No significa la siguiente consulta lo mismo sin ella? La he ejecutado con los datos que me has proporcionado.

INSERT INTO parcel_jurisdictions(parcel_gid,jurisdiction_gid,isect_geom)
  SELECT
    a.orig_gid parcel_gid,
    b.orig_gid jurisdiction_gid,
    a.geom
  FROM
    valid_parcels a, valid_jurisdictions b
  WHERE
    st_intersects(a.geom,b.geom);

Configuración

Processor: AMD Athlon II X4 635 2.9 GHz 
Memory: 4 GB
Platform: Windows 7 Professional
Potgres Version: 8.4
Postgis Version: "POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER"

Analizar los resultados

"Nested Loop  (cost=0.00..7505.18 rows=217489 width=1580) (actual time=1.994..248405.616 rows=329150 loops=1)"
"  Join Filter: _st_intersects(a.geom, b.geom)"
"  ->  Seq Scan on valid_jurisdictions b  (cost=0.00..37.25 rows=525 width=22621) (actual time=0.054..1.732 rows=525 loops=1)"
"  ->  Index Scan using valid_parcels_index on valid_parcels a  (cost=0.00..11.63 rows=2 width=1576) (actual time=0.068..6.423 rows=1366 loops=525)"
"        Index Cond: (a.geom && b.geom)"
"Total runtime: 280087.497 ms"