Adquirir velocidad similar a ArcGIS en Postgis

61

He estado usando Postgis 2.0 durante 3/4 de un año y aunque realmente disfruto usarlo, el tiempo de procesamiento de consultas excesivo lo ha hecho prácticamente inutilizable para mi caso de uso.

Tiendo a realizar geoprocesos pesados en conjuntos de datos municipales que a menudo tienen cientos de miles de multipolígonos. Estos multipolígonos a veces tienen una forma muy irregular y pueden variar de 4 puntos a 78,000 puntos por multipolígono.

Por ejemplo, cuando intersecto 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 para el 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, requiere 1500% más de tiempo para hacer esta intersección en Postgis que en ArcGIS . ¡Y esta es una de mis consultas más simples!

Una de las razones por las que ArcGIS se ejecuta más rápido se debe a los mejores índices. Algunos programadores descubrieron recientemente cómo funcionan estos índices y me pregunto si alguien sabe cómo crear estos índices en Postgis (o crear tablas que imitaría los índices). Quizás esto resolvería la mayoría de los problemas de velocidad en Postgis. Solo puedo esperar que haya alguna manera, especialmente porque ArcGIS solo puede usar 4 GB de RAM, mientras que yo podría usar hasta 4 veces más que para mi servidor de Postgis

Por supuesto, hay muchas razones por las que postgis puede ejecutarse lentamente, por lo que 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 el proceso de instalación completo que utilicé para configurar Postgis, incluida la creación de la máquina virtual en sí .

También aumenté la memoria compartida de los 24 MB predeterminados 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í hay enlaces a los datos que utilicé para esta prueba:

  1. Parcelas: tcad_parcels_06142012.shp.zip desde < a href="ftp://ftp.ci.austin.tx.us/GIS-Data/planning/temp/"> Ciudad de Austin, TX
  2. Jurisdicciones: Límites jurisdiccionales de Ciudad de Austin, TX

Estos son los pasos que tomé para procesar los datos:

ArcGIS

  1. Agregar conjuntos de datos a ArcMap
  2. Establezca el sistema de coordenadas en pies centrales de texas (srid 2277)
  3. Use la herramienta de intersección del menú desplegable

Postgis

Importar parcelas usando:

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 usando:

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 geometría no válida en 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 geometría inválida en 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;

Ejecutar clúster:

cluster;

Ejecutar análisis de vacío:

vacuum analyze;

Realizar intersección en tablas 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 Analizar consulta de intersección:

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 lo que he leído, mi consulta de intersección es eficiente y no tengo ni idea de lo que estoy haciendo mal para que la consulta tome 56 minutos en geometría limpia.

    
pregunta THX1138 11.08.2012 - 23:59

2 respuestas

86

Enfoque diferente. Sabiendo que el dolor está en ST_Intersection, y que las pruebas de verdadero / falso son rápidas, tratar de minimizar la cantidad de geometría que pasa a través de la intersección podría acelerar las cosas. Por ejemplo, las parcelas que están totalmente contenidas en una jurisdicción no necesitan recortarse, pero ST_Intersection probablemente se tomará la molestia de construir parte de la superposición de intersección antes de darse cuenta de que no tiene que generar ninguna nueva geometría. 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.

    
respondido por el Paul Ramsey 15.08.2012 - 22:43
4

¿Qué pasaría si omites la parte "st_multi(st_intersection(a.geom,b.geom))" ?

¿La siguiente consulta no significa lo mismo sin ella? Lo ejecuté en los datos que proporcionaste.

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 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"
    
respondido por el vinayan 15.08.2012 - 07:38

Lea otras preguntas en las etiquetas