Nabycie prędkości podobnej do ArcGIS w Postgis

62

Korzystam z Postgis 2.0 już od 3/4 roku i chociaż naprawdę cieszę się z jego używania, nadmierny czas przetwarzania zapytań sprawił, że jest on praktycznie bezużyteczny w moim przypadku użycia.

Często wykonuję ciężkie geoprzetwarzanie w miejskich zestawach danych, które często mają setki tysięcy wieloboków. Te multipolygony są czasami kształtowane bardzo nieregularnie i mogą różnić się od 4 punktów do 78 000 punktów na multipolygon.

Na przykład, gdy przecinam zestaw danych działki z 329 152 wielobokami z zestawem danych jurysdykcji zawierającym 525 wieloboków, otrzymuję następujące statystyki dotyczące całkowitego czasu pracy:

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

Innymi słowy, potrzeba 1500% więcej czasu na wykonanie tego skrzyżowania w Postgis niż w ArcGIS - i to jest jedno z moich prostszych zapytań!

Jednym z powodów, dla których ArcGIS prawdopodobnie działa szybciej, są lepsze indeksy. Niektórzy programiści niedawno odkryli, jak działają te indeksy i zastanawiam się, czy ktoś wie, jak zbudować te indeksy w Postgis (lub zbudować tabele, które naśladują indeksy). Być może rozwiązałoby to większość problemów z prędkością w Postgis. Mogę tylko mieć nadzieję, że musi być jakiś sposób, zwłaszcza, że ​​ArcGIS może zużywać tylko 4 GB pamięci RAM, podczas gdy ja mogłem użyć nawet 4 razy więcej niż na moim serwerze Postgis!

Oczywiście istnieje wiele powodów, dla których postgis może działać wolno, dlatego przedstawię szczegółową wersję specyfikacji mojego systemu:

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

Szczegółowo opisuję również cały proces instalacji, którego użyłem do skonfigurowania postgis, w tym stworzenie samej maszyny wirtualnej .

Zwiększyłem również pamięć współdzieloną z domyślnych 24 MB do 6 GB w pliku conf i uruchomiłem następujące polecenia, aby umożliwić uruchomienie 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

O ile mogę powiedzieć, nie robi to absolutnie nic zauważalnego pod względem wydajności.

Oto linki do danych, które wykorzystałem do tego testu:

  1. Paczki: tcad_parcels_06142012.shp.zip z City of Austin, Teksas
  2. Jurysdykcje: Granice jurysdykcji z miasta Austin, Teksas

Oto kroki, które podjąłem, aby przetworzyć dane:

ArcGIS

  1. Dodaj zestawy danych do ArcMap
  2. Ustaw układ współrzędnych na środkowe stopy w Teksasie (srid 2277)
  3. Użyj narzędzia przecięcia z menu rozwijanego

Postgis

Importuj paczki za pomocą:

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

Importuj jurysdykcje, używając:

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

Wyczyść niepoprawną geometrię w paczkach:

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;

Wyczyść nieprawidłową geometrię w jurysdykcjach:

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;

Uruchom klaster:

cluster;

Przeprowadź analizę próżni:

vacuum analyze;

Wykonaj skrzyżowanie na oczyszczonych tabelach:

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

Wyjaśnij Analizuj zapytanie o przecięcie:

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)

Ze wszystkiego, co przeczytałem, moje zapytanie o przecięcie jest skuteczne i nie mam absolutnie pojęcia, co robię źle, aby zapytanie zajęło 56 minut na czystej geometrii!

THX1138
źródło
2
W PostGIS powszechny jest idiom polegający na dodawaniu sprawdzania przecięcia obwiedni w celu przyspieszenia. Spróbuj dodać „AND a.geom && b.geom” do klauzuli WHERE i przekonaj się, ile to robi różnicy.
Sean
2
st_intersects () zawiera zapytanie pola ograniczającego przed wykonaniem jakichkolwiek testów skrzyżowań w Postgis 2.x, więc niestety nie zaoszczędzi to czasu.
THX1138,
1
Czy możesz uruchomić zapytanie za pomocą EXPLAIN ANALYZE i opublikować wyniki
Nathan W
1
powinieneś również pamiętać, że korzystasz z różnych zestawów danych na postgis vs arcgis, ponieważ mówisz, że musisz je zatwierdzić, aby zostały zaakceptowane przez postgis.
Nicklas Avén,
2
Czy można uzyskać zbiory danych, aby na nie spojrzeć?
Nicklas Avén,

Odpowiedzi:

87

Odmienne podejście. Wiedząc, że ból występuje w ST_Intersection i że testy prawda / fałsz są szybkie, próba zminimalizowania ilości geometrii przechodzącej przez skrzyżowanie może przyspieszyć. Na przykład paczki, które są całkowicie objęte jurysdykcją, nie muszą być przycinane, ale ST_Intersection prawdopodobnie nadal będzie miał problem z budowaniem części nakładki przecięcia, zanim zda sobie sprawę, że nie musi generować żadnej nowej geometrii. Więc to

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

A nawet 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)

Może nawet być szybszy bez UNII.

Paul Ramsey
źródło
13
Dzięki, że sprowadziłem mnie do 3,63 minuty! Nigdy bym nie pomyślał, że związek będzie szybszy. Ta odpowiedź naprawdę sprawi, że ponownie pomyślę o sposobie, w jaki teraz zadaję zapytania.
THX1138,
2
To jest bardzo fajne. Miałem przypadek w pracy, w którym moje zapytanie st_intersection trwało ponad 30 minut i teraz wiem, jak mogę tego uniknąć :)
Nathan W
1
to pytanie skłoniło mnie do nauki Postgis! dziś dobrze spać, widząc Postgisa biegającego ramię w ramię z Arcgisem :-)
vinayan
2
Jeszcze jedno ulepszenie od Martina Davisa, możesz wpisać „in or out?” pytaj do SELECT za pomocą instrukcji CASE i unikaj w ten sposób UNION.
Paul Ramsey,
2
UNIONeliminuje zduplikowane wiersze z dwóch zapytań, ale te dwa zapytania nie mogą mieć tego samego wiersza w zestawie wyników. Tak więc UNION ALL, które pomija sprawdzanie duplikatów, byłoby tutaj właściwe. (Nie używam UNIONzbyt wiele, ale generalnie stwierdzam, że poza czasami używam, znacznie częściej używam UNION ALL.)
jpmc26,
4

Co by się stało, gdybyś pominął tę "st_multi(st_intersection(a.geom,b.geom))"część?

Czy poniższe zapytanie nie oznacza tego samego bez niego? Uruchomiłem to na danych, które podałeś.

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

Konfiguracja

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"

Analizuj wyniki

"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"
vinayan
źródło
Nie, chce wynikowych wielokątów przecięcia, ale twoje zapytanie bardzo ładnie pokazuje, że cały ból powstaje podczas generowania skrzyżowań, a nie w binarnej części zapytania typu prawda / fałsz. Jest to dość oczekiwane, ponieważ kod prawda / fałsz jest wysoce zoptymalizowany, podczas gdy generowanie przecięcia nie jest.
Paul Ramsey,