Nabycie prędkości podobnej do ArcGIS w Postgis

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!

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.

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"