Klasteryzacja zbioru punktów w PostGIS

Celem projektu jest zbadanie możliwości zaimplementowania algorytmu klasteryzacji punktów pod stronie bazy danych. Problem wiele aplikacji jest wyświetlanie zbyt wielu informacji na wybranym poziomie szczegółowości mapy. Jednym z rozwiązań tego zagadnienia jest grupowanie sąsiednich punktów w jeden i wyświetlanie ich dopiero na wyższym poziomie szczegółowości. Rozwiązanie to pozwala na przejrzystsza prezentacje, wyświetlanie zbyt wielu szczegółów wprowadza w zakłopotanie i nie pozwala skupić wzroku skupiającego.

Istnieje wiele algorytmów klasteryzacji. Trudno wskazać optymalny algorytm dla wszystkich zastosować, wybóru należy dokonać indywidualnie po porzedzającej analizie plusów i minusów odpowiednich algorytmów. Większość o ile nie wszystkie istniejące algorytmy realizowane są po stronie aplikacji - pobierane zostają dane z bazy danych a następnie odpowiednio przetwarzane. Zadanie postawione przed projektem, polega na stworzeniu implementacji po stronie serweru bazodanowego. Podłożem i platformą bazodanową która wykorzystywana jest przy realizacji projektu jest Postgres wraz z platformą Postgis - odpowiedzialną za udostępnianie struktór i funkcji dla danych przestrzennych.

• klasteryzaja punktów
• określenie warunku stopu
• określenie warunków brzegowych, punków które zawierają się w zadanym prostokącie
• określenie maksymalnej ilości punktów

Projekt posiada jeden przypadek użycia. Zaimplementowany mechanizm pozwala na pobranie punktów dla danego obszaru. Pozwala on na ustalenie gęstości punktów które zostaną wyświetlone na mapie. Wszelkie obszary które cechują się za dużą gęstością ( odległości pomiędzy istniejącymi punktami są mniejsze od maksymalnie zadanej) grupowanę są w elementy zwane klustrami. Każdy kluster posiada swoją współrzędną geograficzną ( będącą średnią zawieranych punktów) i liczbę prezentującą ilość zawieranych punktów.

Wykorzystywane będą funkcje Postgis i możliwość procedur składowych języka Postgres. Udostępnione zostaną procedury pozwalające na wykonanie klasteryzacji dla zadanych parametrów.

• Przechowywanie danych geograficznych w strukturze geometry, należy pamiętać iż przechowywane są one w sposób binarny. Tak więc przeglądając bazę danych poprzez interfejs bazodanowy, zawartość kolumn geometry nie może zostać odczytana przez człowieka
• wypisywanie danych w postaci możliwej dla odczytu dla człowieka za pomocą metody ST_AsText
• porównywanie punktów za pomocą operatorów dwuargumentowych &< ( znajduje punkty o mniejszej szerokości geograficznej) |&> ( znajduje punkty o większej wysokości geograficznej) , &> ( znajduje punkty o większej szerokości geograficznej) AND coords &<| ( znajduje punkty o mniejszej wysokości geograficznej)
• Pobieranie szerokości punktu xmin, pobieranie wysokości punku ymin
• Ujednolicanie SRID punktów, konwersja dokonywana jest za pomocą ST_Transform
• Znajdowanie odległości między punktami, jednostka zależna od SRID punktu ST_Distance
• Tworzenie nowych punktów i geometrii - geometryFromText
• Znajdywanie punktu ciężkości dla dwóch punktów - ST_Line_Interpolate_Point, gdzie pierwszym argumentem jest odcinek posiadający jako wierzchołki dwa punty, a drugim argumentem 0.5 czyli połowa odległości
• Dla projektu wybrany został SRID 900913, jest to SRID używany przez Google Maps, podający odległości w metrach

Aplikacja napisana została w języku /pgSQL Działanie aplikacji przebiega w następującym porządku:

1. Pobranie z bazy punktów dla zadanych ograniczeń (lewy górny wierzchołek prostokąta, prawy dolny wierzchołek prostokąta, limit)
2. Wpisanie punktów do tablicy
3. Wykonanie klasteryzacji
4. Zwrócenie wyniku

Całość obsługiwana jest przez procedurę która jest kontrolerem

Struktura ta zawiera w sobie współrzędne i liczbę zawieranych punktów

CREATE TYPE _cluster AS (count_ int, point_ geometry);

-- kontroller, zarządza wszystkim

CREATE OR REPLACE FUNCTION _cluster_controller(lpoint geometry, rpoint geometry, stop_dist int, limit_ int) RETURNS _cluster[] AS $$
DECLARE
 _clusters _cluster[];
 tmp integer;
BEGIN
 _clusters := _get_clusters(lpoint, rpoint,limit_);
--PERFORM show_clusters(_clusters);
 _clusters := calculate_distance(_clusters, stop_dist);
--PERFORM show_clusters(_clusters);
 RETURN _clusters;

END
$$
LANGUAGE plpgsql;

Pobrane punkty przechowywane są w strukturze, która reprezentuje klaster. Zastosowane więc zostaje założenie iż każdy punkt jest początkowo klasterem.

--Pobiera punkty i wpisuje je do clusterów, zwraca je

CREATE OR REPLACE FUNCTION  _get_clusters(lpoint geometry, rpoint geometry, limit_ int) RETURNS _cluster[] AS $$
DECLARE
   my_arr _cluster[]; -- tablica poczatkowa
   tmp_clust _cluster;
   counter integer:= 0;
   i integer := 0;
BEGIN
FOR tmp_clust IN SELECT 1 as count,coords  FROM opimap_point WHERE coords &< rpoint AND coords |&> rpoint AND coords &> lpoint AND coords &<| lpoint LIMIT limit_  LOOP -- get points from area
       my_arr[counter] := tmp_clust; -- initalize array
       counter :=counter+1;
        ---RAISE NOTICE 'Counter is %', counter;  -- Prints 80
        -- RAISE NOTICE 'Geomenty is %', tmp_clust.point_;  -- Prints 80
   END LOOP;
 WHILE i < counter LOOP -- loop, calculate 
 -- RAISE NOTICE 'LOOP POINT is %', ST_AsText(my_arr[i].point_);  -- Prints 80
 i := i +1;
END LOOP;

        -- RAISE NOTICE 'wielkosc tablicy to %', array_length(my_arr, 1);  -- Prints 80


   RETURN my_arr;
END
$$
LANGUAGE plpgsql;

1. O tworzymy z klastrów macierz trójkątną, na której wartościami jest odległość między środkami ciężkości klastrów ( czyli początkowo między punktami)
2. znajdujemy parę dla której wartość ( czyli odległość jest najmniejsza ), dla pierwszej iteracji będą to para punktów znajdujących się najbliżej
3. Sprawdzany jest warunek stopu, jeżeli dla odległość dla pary znalezionej w punkcie 2, jest mniejsza od zadanej maksymalnej (warunek stopu) wykonujemy pozostałe punkty. Jeżeli odległość minimalna jest większa od warunku stopu, uznajemy iż klasteryzacja została zakończona i zwracamy wynik
4. para która posiada minimalną odległóść dla aktualnej iteracji tworzy nowy klaster. Współrzędnymi klastra jest środek ciężkości pomiędzy dwoma klastrami które ją tworzą. Ilość punków jest sumą punktów dwóch tworzących ją klastrów
5. Następnie rekurencyjnie, dla nowego zbioru klastrów (n-1), powtarzany jest cały algorytm aż do punktu stopu

Poniżej przedstawione jest działanie algorytmu na rysunku

-- znajduje najmniejszy dystans między klustrami, max stop_dist is 999999

CREATE OR REPLACE FUNCTION  calculate_distance(_clusters _cluster[], stop_dist int) RETURNS _cluster[] AS $$
DECLARE
   counter integer := 0;
   i integer := 0;
   j integer := 0;
   i_min integer := 0;
   j_min integer := 0;
 new_count integer := 0;
   min_distance float := 999999;
   tmp_distance float := 0;
   middle_geometry geometry;
   new_arr _cluster[];
BEGIN
   counter := array_length(_clusters, 1);
   WHILE i < counter LOOP -- loop, calculate 
     j := i + 1;
     WHILE j < counter LOOP
       
	         -- RAISE NOTICE 'computed points are  % and %', ST_AsText(_clusters[i].point_),ST_AsText(_clusters[j].point_);  -- Prints 80 
 -- RAISE NOTICE 'distance is  %',  ST_Distance(    ST_Transform(  _clusters[i].point_,900913 ), ST_Transform( _clusters[j].point_,900913)  )  ;  

  tmp_distance := ST_Distance(    ST_Transform(  _clusters[i].point_,900913 ), ST_Transform( _clusters[j].point_,900913)  )  ;

IF tmp_distance < min_distance THEN
min_distance = tmp_distance;
i_min = i;
j_min = j;
END IF;
          
     
       j := j + 1;
     END LOOP;	
      i := i + 1;
		
END LOOP;
-- RAISE NOTICE 'COUNTER is  %', counter  ; 
-- RAISE NOTICE 'MINIMAL distance is  %',  min_distance  ;  
IF min_distance < stop_dist THEN

 -- RAISE NOTICE 'MINIMAL points are  % and %', ST_AsText(_clusters[i_min].point_),ST_AsText(_clusters[j_min].point_);  -- Prints 80 


middle_geometry := ST_Line_Interpolate_Point(geometryFromText(

       'LINESTRING('|| xmin(_clusters[i_min].point_)||' ' || ymin(_clusters[i_min].point_) ||','||
xmin(_clusters[j_min].point_) || ' '||  ymin(_clusters[j_min].point_) || ')', 4326
    ), 0.5);

-- RAISE NOTICE 'MIDDLE POINT  is  %', ST_AsText( middle_geometry );
-- RAISE NOTICE 'Distance From Point A  %', ST_Distance(     ST_Transform(middle_geometry,900913), ST_Transform(_clusters[i_min].point_,900913)  )  ;
-- RAISE NOTICE 'Distance From Point B  %', ST_Distance(      ST_Transform(middle_geometry,900913), ST_Transform(_clusters[j_min].point_,900913)  )  ;
-- RAISE NOTICE 'I  is  % AND J is %', i_min, j_min;


i := 0;
j:= 0;
    WHILE i < counter LOOP
	IF (i != i_min) AND (i != j_min) 
        THEN
	  new_arr[ j] := _clusters[i];
	  j :=  j + 1;
          
        END IF;	
	i := i + 1;
    END LOOP;
new_count := _clusters[i_min].count_ + _clusters[j_min].count_;
--RAISE NOTICE '@!!! new count = % hey', new_count;

    new_arr[j] =  ROW(new_count,    middle_geometry );

--PERFORM show_clusters(new_arr);
    RETURN calculate_distance(new_arr ,stop_dist);

ELSE
RETURN _clusters ;

END IF;


  
END
$$
LANGUAGE plpgsql;

Na etapie tworzenia aplikacje, ważnym aspektem była kontrola poprawności. Stworzona została pomocnicza funkcja show, której zadaniem jest wypisanie przesłanej na wejściu tablicy.

-- wypisuje klustry

CREATE OR REPLACE FUNCTION  show_clusters(_clusters _cluster[]) RETURNS integer AS $$
DECLARE
counter integer;
i integer := 0;
BEGIN
   counter := array_length(_clusters, 1);
   WHILE i < counter LOOP -- loop, calculate 
     
RAISE NOTICE '!!  % cluster geometry =  % and count = %',i, ST_AsText(_clusters[i].point_),_clusters[i].count_;  -- Prints 80 
i := i +1;
END LOOP;

RETURN 1;
  
END
$$
LANGUAGE plpgsql;


-- CREATE TYPE _cluster AS (count_ int, point_ geometry);
-- CREATE TYPE _xy_cluster AS ( count_ int, x float, y float)

Początkowo stworzony kontroler, zwraca strukturę _cluster. Zawiera ona ilość punktów i współrzędne geograficzne przechowywane jako geometry. Przechowywanie współrzędnych w tej postaci jest odpowiednie dla pozostałych metod używanych w aplikacji. Format binarny jest jednak nie do zaakceptowania dla aplikacji klienckiej, oczekuje ona prostej prezentacji - współrzędnych X i Y. Dodatkowym utrudnieniem jest zwracanie przez kontroler danych w postaci tablicy (array []), format ten nie jest obsługiwany w aplikacji klienckiej - wymaga napisania specjalnego parsera.

Stworzony została więc pomocnicza struktura _xy_cluster i kontroler który zwraca ją w postaci wierszy

CREATE TYPE _xy_cluster AS ( count_ int, x float, y float)

Znosi on następujące wady wcześniejszego kontrolera:

• zwraca strukturę przekonywującą współrzędne geograficzne jako x i y
• zwraca kolekcje w postaci wierszy ( nie tablicy )

-- kontroller, zarządza wszystkim, xy

CREATE OR REPLACE FUNCTION _xy_cluster_controller(lpoint geometry, rpoint geometry, stop_dist int, limit_ int) RETURNS SETOF _xy_cluster AS $$
DECLARE
 _clusters _cluster[];
   counter integer := 0;
   i integer := 0;
   tmp _xy_cluster;
BEGIN

   _clusters = _cluster_controller(lpoint, rpoint,stop_dist, limit_ );
   counter := array_length(_clusters, 1);
   WHILE i < counter LOOP -- loop, calculate 
	tmp = ROW(_clusters[i].count_,   xmin(_clusters[i].point_),  ymin(_clusters[i].point_) );
	RETURN NEXT tmp;
    i := i + 1;
    
     END LOOP;	
 RETURN;

END
$$
LANGUAGE plpgsql;

Aby zagwarantować działanie aplikacji musimy dysponować serwerem z zainstalowanym PostgreSQL i Postgis. Struktura tabeli przechowującej punkty które będą klasteryzowane może być dowolna. Wymagane jest jedynie by kolumna przechowująca dane geograficzne była typu geometry. Następnie należy dostosować pobieranie punktów poprzez edycje procedury get_cluster(); Instalacje polega na wgraniu procedur do silnika bazodanowego.

Klasteryzacja 400 punktów dla warunku stopu 1000 metrów

SELECT _xy_cluster_controller(ST_GeometryFromText('srid=4326;POINT(49 20)'), ST_GeometryFromText('srid=4326;POINT(51 18)'), 1000,400);

Klasteryzacja dla 20 punktów przy warunku stopu 500 metrów

SELECT _xy_cluster_controller(ST_GeometryFromText('srid=4326;POINT(49 20)'), ST_GeometryFromText('srid=4326;POINT(51 18)'), 500,20);

Testowane były czasy wywołania. Testy podzielono na dwie grupy: dla warunku stopu równego 1000 metrów i warunku stopu równego 1 metr. Dla każdej z grup wykonano testy odpowiednio dla 10, 50, 100, 200, 300 i 400 punktów. Dla wykonanych testów, w procedurze get_cluster, wyświetlane były 2 wiadomości typu notice - o aktualnej minimalnej odległości i ich współrzędnych ich klastrów

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=4.745..4.761 rows=9 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1000, 10)"
"Total runtime: 4.797 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=191.030..191.078 rows=34 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1000, 50)"
"Total runtime: 191.145 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=1862.090..1862.173 rows=55 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1000, 100)"
"Total runtime: 1862.274 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=16864.409..16864.552 rows=104 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1000, 200)"
"Total runtime: 16864.714 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=67506.663..67506.865 rows=145 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1000, 300)"
"Total runtime: 67507.061 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=205355.356..205355.586 rows=170 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1000, 400)"
"Total runtime: 205355.813 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=3.015..3.031 rows=10 loops=1)"
"Total runtime: 3.104 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=18.639..18.712 rows=50 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1, 50)"
"Total runtime: 18.806 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=64.056..64.202 rows=100 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1, 100)"
"Total runtime: 64.349 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=249.291..249.560 rows=200 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1, 200)"
"Total runtime: 249.797 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=590.032..590.433 rows=300 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1, 300)"
"Total runtime: 590.807 ms"

"Result  (cost=0.00..0.26 rows=1 width=0) (actual time=1127.768..1128.336 rows=400 loops=1)"
"  Output: _xy_cluster_controller('0101000020E610000000000000008048400000000000003440'::geometry, '0101000020E610000000000000008049400000000000003240'::geometry, 1, 400)"
"Total runtime: 1128.840 ms"

Analiza otrzymanych wyników wskazuje na wykładniczy wzrost czasu wykonania klasteryzacji wraz ze wzrostem ilości punktów. Można optymalizować algorytm, usunąć wszelkie miejsca gdzie dokonywane jest wypisanie na ekran - spowoduje to jednak jedynie zmniejszenie współczuwając, jednak samo wykonanie przebiegać będzie wykładniczo. Podsumowując, klasteryzacja wykonywana jest przy akceptowalnym progu czasowym dla małej ilości punktów, dla większej ich ilości należy dobrać inny algorytm klasteryzacji.