Memperoleh kecepatan seperti ArcGIS di Postgis

Saya telah menggunakan Postgis 2.0 selama 3/4 tahun sekarang dan sementara saya benar-benar menikmati menggunakannya, waktu pemrosesan kueri yang berlebihan pada dasarnya membuatnya tidak dapat digunakan untuk kasus penggunaan saya.

Saya cenderung melakukan geoproses berat pada dataset kota yang sering memiliki ratusan ribu multipoligon. Multipoligon ini terkadang berbentuk sangat tidak teratur dan dapat bervariasi dari 4 poin hingga 78.000 poin per multipoligon.

Misalnya, ketika saya memotong dataset parcel dengan 329.152 multipolygons dengan dataset yurisdiksi yang mengandung 525 multipolygons, saya mendapatkan statistik berikut untuk total waktu yang digunakan:

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

Dengan kata lain, ini membutuhkan 1500% lebih banyak waktu untuk melakukan persimpangan di Postgis ini daripada di ArcGIS - dan ini adalah salah satu pertanyaan saya yang lebih sederhana!

Salah satu alasan mengapa ArcGIS berjalan lebih cepat adalah karena indeks yang lebih baik. Beberapa programmer baru-baru ini menemukan cara kerja indeks ini dan saya bertanya-tanya apakah ada yang tahu bagaimana membangun indeks ini di Postgis (atau membuat tabel yang akan meniru indeks). Mungkin ini akan menyelesaikan sebagian besar masalah kecepatan di Postgis. Saya hanya bisa berharap harus ada beberapa cara, terutama karena ArcGIS hanya dapat menggunakan 4 GB RAM sementara saya bisa menggunakan hingga 4 kali lipat dari server postgis saya!

Tentu saja ada banyak alasan postgis dapat berjalan lambat, jadi saya akan memberikan versi detail dari spesifikasi sistem saya:

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

Saya juga merinci seluruh proses instalasi yang saya gunakan untuk mengatur postgis termasuk pembuatan VM itu sendiri .

Saya juga meningkatkan memori bersama dari 24MB default ke 6 GB di file conf dan menjalankan perintah berikut untuk memungkinkan postgres berjalan:

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

Sejauh yang saya tahu ini, tidak ada yang terlihat dalam hal kinerja.

Berikut ini tautan ke data yang saya gunakan untuk tes ini:

1. Paket: tcad_parcels_06142012.shp.zip dari City of Austin, TX
2. Yurisdiksi: Batas Yurisdiksi dari City of Austin, TX

Berikut langkah-langkah yang saya ambil untuk memproses data:

ArcGIS

1. Tambahkan dataset ke ArcMap
2. Atur sistem koordinat ke kaki texas pusat (srid 2277)
3. Gunakan alat persimpangan dari menu dropdown

Postgis

Impor paket menggunakan:

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

Impor yurisdiksi menggunakan:

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

Bersihkan geometri yang tidak valid dalam paket:

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;

Bersihkan geometri yang tidak valid di yurisdiksi:

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;

Jalankan cluster:

cluster;

Jalankan analisis vakum:

vacuum analyze;

Lakukan persimpangan pada meja yang dibersihkan:

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);

Jelaskan Menganalisis kueri persimpangan:

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)

Dari semua yang saya baca, kueri persimpangan saya efisien dan saya sama sekali tidak tahu apa yang saya lakukan salah untuk kueri membutuhkan waktu 56 menit pada geometri bersih!

Pendekatan yang berbeda. Mengetahui bahwa rasa sakit ada di ST_Intersection, dan bahwa tes benar / salah cepat, mencoba untuk meminimalkan jumlah geometri yang melewati persimpangan mungkin mempercepat segalanya. Misalnya, parsel yang benar-benar terkandung dalam yurisdiksi tidak perlu dipotong, tetapi ST_Intersection mungkin masih akan mengalami kesulitan membangun bagian dari hamparan persimpangan sebelum menyadari bahwa itu tidak harus menghasilkan geometri baru. Jadi ini

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);

Atau bahkan 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)

Bahkan bisa lebih cepat tanpa UNION.

Apa yang akan terjadi jika Anda menghilangkan "st_multi(st_intersection(a.geom,b.geom))"bagian itu?

Bukankah permintaan di bawah ini memiliki arti yang sama tanpanya? Saya menjalankannya pada data yang Anda berikan.

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);

Konfigurasi

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"

Analisis Hasil

"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"