• Manuale PostGIS

Installando il pacchetto Debian postgresql-8.2-postgis (versione PostGIS 1.2.1 per PostgreSQL 8.2) il template_gis non viene creato. Ecco un riassunto delle istruzioni per crearlo. Notare che ci si rivolge alla porta TCP 5433 invece che alla standard 5432, questo perché sul server è stato installato Postgres 8.2 insieme a Postgres 8.1; il pacchetto Debian ha provveduto a configurare i due server su due porte differenti (vedere alcune note in PostgreSQL cluster):

su - postgres
PATH=/usr/lib/postgresql/8.2/bin:$PATH
createdb --port 5433 --template template0 template_gis
createlang --port 5433 plpgsql template_gis
psql --port 5433 -f /usr/share/postgresql-8.2-postgis/lwpostgis.sql -d template_gis
psql --port 5433 -f /usr/share/postgresql-8.2-postgis/spatial_ref_sys.sql -d template_gis

Per trasformare il database in un template ci si collega all'SQL (/usr/lib/postgresql/8.2/bin/psql –port 5433) e si esegue:

UPDATE pg_database SET datistemplate='t' WHERE datname = 'template_gis';

ATTENZIONE: per motivi ignoti a volte il restore eseguito con psql fallisce, sia che per leggere il file dump si usi l'opzione -f, oppure la ridirezione da shell <, oppure ancora il comando interno \i.

Sembra che il problema stia nel comando COPY su record molto grossi o che contengono caratteri particolari. Alcuni errori possibili sono:

ERROR:  syntax error at or near "5" at character 1
Invalid command \..
Invalid command \N.

La soluzione è utilizzare pg_dump -Fc (dump in formato custom invece che plain-text) e quindi effettuare il restore con pg_restore invece di psql. Anche l'opzione --inserts di pg_dump risolve il problema, ma è sconsigliata per l'estrema lentezza del restore.

Quando si aggiorna la versione di PostGIS il database geografico deve essere ricostruito. Mentre i dati utente possono essere salvati con un dump e poi recuperati con un restore, gli oggetti PostGIS (funzioni, tipi, ecc.) devono essere aggiornati alla nuova versione e quindi non vanno recuperati dal dump. Anche la tabella spatial_ref_sys dovrebbe essere ricostruita con l'ultima versione del file sql fornito con PostGIS.

Quella che segue è una ricetta per effettuare un upgrade, è stata sperimentata nel passaggio:

PostgreSQL 8.1, PostGIS 1.1.6 ⇒ PostgreSQL 8.2, PostGIS 1.2.1

Nel nostro caso il vecchio database sta in ascolto sulla porta TCP 5432, mentre quello nuovo è sulla porta 5433. Si consiglia di eseguire i vari comandi PostgreSQL (psql, pg_dump, pg_restore) richiamandoli con il percorso completo (/usr/lib/postgresql/<ver>/bin/) rispettando la versione di database con il quale di volta in volta si interagisce.

Anzitutto si esegue il dump del vecchio database (notare l'opzione -Fc):

su - postgres
pg_dump --port 5432 -Fc dbname > dbname.dump

Poi da utente privilegiato si crea il nuovo database a partire dal nuovo template_gis, in modo da avere l'ultima versione dell'estensione PostGIS:

su - postgres
psql --port 5433

CREATE DATABASE dbname OWNER dbuser TEMPLATE template_gis;
\connect dbname
GRANT ALL ON geometry_columns TO dbuser;
GRANT SELECT ON spatial_ref_sys TO dbuser;

Poi si deve fare il restore solo dei dati utente. Un trucco potrebbe essere quello di effettuare il restore con un utente che non ha i privilegi per sovrascrivere gli oggetti PostGIS; vengono visualizzati molti errori durante il restore, ma il risultato finale dovrebbe essere corretto.

Una soluzione più pulita è quella di effettuare un restore selettivo con pg_restore.

Si estrae la lista degli oggetti contenuti nel dump:

pg_restore --list dbname.dump > dbname_list

Si edita la lista degli oggetti (dbname_list) togliendo quelli generati da PostGIS. In pratica si tolgono tutti i riferimenti alla tabella spatial_ref_sys e le righe che contengono:

ACL
AGGREGATE
CAST
COMMENT
FUNCTION
OPERATOR
PROCEDURAL LANGUAGE
SCHEMA
TYPE

Quindi si effettua il restore selettivo da utente non privilegiato:

pg_restore --port 5433 -U dbuser -W -h localhost -L dbname_list -d dbname dbname.dump

Se si eseguono query su oggetti geografici conviene creare un indice sulla colonna geografica. PostGIS utilizza in questo caso indici di tipo GiST. Ecco come creare l'indice (tracksegments è il nome della tabella e trkseg è il campo con la geometria):

CREATE INDEX tracksegments_trkseg_idx ON tracksegments USING GIST (trkseg GIST_GEOMETRY_OPS);
VACUUM ANALYZE tracksegments (trkseg);

Probabilmente l'operatore GIST_GEOMETRY_OPS può essere omesso, essendo l'unico disponibile.

Importante il comando VACUUM ANALYZE, altrimenti i successivi query plans non saranno ottimizzati. Su un campione di 1061 record di tipo LINESTRING (per un totale di circa 315000 punti) una query senza indice impiega circa 29 ms:

EXPLAIN ANALYZE SELECT count(trkseg) FROM tracksegments
WHERE trkseg && SetSRID(MakeBox2D(MakePoint(10.588, 43.026), MakePoint(11.661, 43.771)), 4326);
 Aggregate  (cost=52.00..52.01 rows=1 width=4548) (actual time=29.292..29.293 rows=1 loops=1)
   ->  Seq Scan on tracksegments
   ...
   ...
 Total runtime: 29.343 ms

dopo la creazione dell'indice impiega un terzo del tempo:

 Aggregate  (cost=51.26..51.27 rows=1 width=4794) (actual time=9.720..9.722 rows=1 loops=1)
   ->  Bitmap Heap Scan on tracksegments  (cost=9.96..50.71 rows=220 width=4794)
   ...
   ...
 Total runtime: 9.783 ms

Il template_gis (template_gis#SONAME#) non viene creato in automatico.

Esiste la documentazione on-line per PostGIS all'indirizzo http://postgis.refractions.net/docs/, utile soprattutto per la reference delle funzioni. Sono disponibili anche alcune note sulla versione 0.9. Interessanti anche queste dispense del Mapserver User Meeting.

Sono richiesti anzitutto i pacchetti postgresql-contrib e libpg-perl. La libreria postgis ancora non fa parte della distribuzione ufficiale Debian, prelevato il pacchetto libpostgis1-pg74_1.0.0RELEASE-1_i386.deb dal sito pkg-grass e installato.

Durante l'installazione del pacchetto viene creato un template_gis che conviene utilizzare come base per creare tutti i database gemoetrici. Su questo template sono state fatte automaticamente alcune operazioni preliminari:

• Creazione del linguaggio plpgsql
• Creazione delle funzioni PostGIS
• Creazione della tabella geometry_columns
• Creazione e popolazione della tabella spatial_ref_sys

NOTA: Alcune operazioni su PostgreSQL (creare un database, attivare il linguaggio PL/pgSQL, creare le estensioni PostGIS) bisogna eseguirle come utente privilegiato postgres del database. Con l'installazione predefinita di Debian è sufficiente operare come utente di sistema postgres (da root effetturare su - postgres), in quanto il database effettua l'autenticazione ident sameuser.

Per gli esempi che seguono si presuppone che PostgreSQL accetti anche l'autenticazione con password (md5) per chi si connette via TCP/IP (host), aggiungere al file /etc/postgresql/pg_hba.conf una riga del tipo:

host all all 127.0.0.1 255.255.255.255 md5

Ora si può creare un database geometrico. La procedura completa può essere riassunta in questo modo:

$ su -
# su - postgres
$ psql template1

template1=# CREATE USER "strade" PASSWORD 'supersecret';
template1=# CREATE DATABASE strade OWNER strade TEMPLATE template_gis ENCODING 'UTF8';
template1=# \connect strade
strade=# ALTER TABLE geometry_columns OWNER TO strade;
strade=# ALTER TABLE spatial_ref_sys OWNER TO strade;
strade=# \dt
             List of relations
 Schema |       Name       | Type  | Owner
--------+------------------+-------+--------
 public | geometry_columns | table | strade
 public | spatial_ref_sys  | table | strade
(2 rows)

strade=# \q

Per chi interessa, le procedure per attivare le estensioni geometriche senza usare il template_gis sarebbero queste:

$ su -
# su - postgres
$ createlang plpgsql strade
$ psql -f /usr/share/postgresql/contrib/lwpostgis.sql -d strade
$ psql -f /usr/share/postgresql/contrib/spatial_ref_sys.sql -d strade

Creare un linguaggio significa che sarà possibile aggiungere a quel database funzioni o trigger scritti in quel linguaggio. Questa operazione è necessaria perché PostGIS comprende funzioni scritte in PL/pgSQL. In generale aggiungere un linguaggio ad un database richiede privilegi di superutente (postgres) e va fatto per ogni database che lo necessita.

Per default PostgreSQL aggiunge il linguaggio con la clausola TRUSTED, ciò significa che il linguaggio è considerato sicuro (non consente di aggirare i meccanismi di protezione) e quindi potrà essere utilizzato anche dagli utenti non privilegiati. Controllare lo stato di trust del linguaggio con il il comando createlang -l strade. Vedere il manuale SQL a proposito di CREATE LANGUAGE.

Anche lo script lwpostgis.sql deve essere eseguito con i privilegi di superutente (postgres). Il motivo è che vengono definite diverse funzioni di tipo LANGUAGE 'C' fornite dalla shared library /usr/lib/postgresql/lib/liblwgeom.so.1.0; agganciare un eseguibile esterno ad una funzione PostgreSQL è consentito al solo superutente.

Oltre alle funzioni vengono create anche le due tabelle: geometry_columns e spatial_ref_sys (e due relativi indici geometry_columns_pk e spatial_ref_sys_pkey). Questi oggetti appena creati appartengono all'utente postgres e quindi è un problema far usare il database o farne fare un dump all'utente non privilegiato. Si cambia proprietario alle tabelle (e contestualmente agli indici) con le istruzioni ALTER TABLE viste sopra.

Per consentire all'utente non privilegiato di lavorare, sarebbe sufficiente garantire la SELECT sulla tabella spatial_ref_sys e garantire tutti i privilegi (ALL) sulla geometry_columns. Vedere il manuale SQL a proposito di GRANT.

Nella tabella spatial_ref_sys sono memorizzate le caratteristiche dei vari sistemi di riferimento spaziale (usata ad esempio per fare le conversioni) e in generale verrà usata in sola lettura. I dati vengono caricati dallo script spatial_ref_sys.sql fornito con PostGIS.

La tabella geometry_columns contiene informazioni sui campi geometrici creati nel database. Le funzioni PostGIS, ad esempio la AddGeometryColumn(), modificano i record in questa tabella. Le funzioni rese disponibili da PostGIS si elencano dalla shell psql con il comando \df, ma sono mescolate a quelle predefinite. L'unica differenza è che appartengono allo schema public invece che allo schema pg_catalog.

$ pg_dump -h 127.0.0.1 --create -U strade -W strade > strade_plpgsql_postgis.dump

Dentro il dump ci saranno ovviamente le tabelle con i dati dell'utente, ma anche oggetti che devono essere ricreati con i privilegi di superutente (il linguaggio PL/pgSQL, le funzioni di PostGIS). Inoltre se si vuole effettuare l'upgrade ad una versione più recente di PostGIS il passaggio potrebbe non essere così semplice. Vedere in proposito la pagina dump e restore di un database PostGIS.

Vedere un po' come fare il restore, perché alcuni oggetti vanno restorati come superutente, altri possono essere restorati come utente non privilegiato.

Creazione di una tabella per memorizzare waypoint e inserimento di un record:

psql -h 127.0.0.1 -U strade -W strade
strade=> CREATE TABLE waypoints (
           id          SERIAL,
           elevation   float,
           name        varchar(6),
           comment     varchar(20),
           description varchar(50),
           idsymbol    integer
         );
strade=> SELECT AddGeometryColumn('waypoints','wpt',4326,'POINT',2);
strade=> INSERT INTO waypoints (elevation, name, idsymbol, wpt)
         VALUES (46.749512, 'NICCO', 2, GeometryFromText('POINT(11.097390857 43.858187562)',4326));

Nella funzione AddGeometryColumn() il parametro 4326 è lo SRID (spatial referencing system identifier), in questo caso corrisponde al datum WGS84) utilizzato dal GPS. Deve corrispondere allo srid di un record in spatial_ref_sys. L'ultimo parametro è la dimensione, in questo caso si lavora sul piano in due dimensioni. Una colonna di tipo geometrico deve essere aggiunta con la funzione AddGeometryColumn(), perché essa provvede a memorizzare informazioni supplementari in geometry_columns.

Per rimuovere da una tabella una colonna di tipo geometry:

strade=> SELECT DropGeometryColumn('', 'waypoint','wpt');

Creare una tabella con il vettoriale delle strade, ogni record di questa tabella contiene un track segment in un campo di tipo LINESTRING:

CREATE TABLE tracksegments (id SERIAL);
SELECT AddGeometryColumn('tracksegments','trkseg',4326,'LINESTRING',3);
INSERT INTO tracksegments (trkseg) VALUES (
    GeometryFromText('LINESTRING(
        11.132305 43.820643,
        11.131639 43.820643,
        11.132927 43.819184,
        11.134493 43.819592,
        11.136060 43.819249
    )', 4326)
);

Problema: come fare INSERT di una LINESTRING con magari qualche migliaio di punti? Sicuramente con una sola INSERT (specialmente da PHP) si potrebbero avere dei problemi sulla eccessiva lunghezza della istruzione SQL.

Per risolvere il problema di cui sopra si può creare una tabella in cui inserire i singoli trackpoint e poi usare la funzione di aggregazione MakeLine() (assente in PostGIS 0.9) per generare un record nella tabella tracksegments.

CREATE TABLE trackpoints (
    id         serial,
    idtrkseg   integer,
    trkpt_time timestamp without time zone,
    elevation  double precision
);
SELECT AddGeometryColumn('trackpoints','trkpt',4326,'POINT',2);
 
INSERT INTO tracksegments (trkseg) VALUES (NULL);
SELECT CURRVAL('tracksegments_id_seq'); -- Supponiamo che il valore corrente sia 1
 
INSERT INTO trackpoints (idtrkseg, trkpt_time, elevation, trkpt) VALUES
    (1, '2005-04-25 16:23:09', 418.791504,
    GeometryFromText('POINT(10.987358 43.786654)', 4326));
INSERT INTO trackpoints (idtrkseg, trkpt_time, elevation, trkpt) VALUES
    (1, '2005-04-25 16:23:28', 424.559448,
    GeometryFromText('POINT(10.987401 43.786612)', 4326));
INSERT INTO trackpoints (idtrkseg, trkpt_time, elevation, trkpt) VALUES
    (1, '2005-04-25 16:23:32', 426.001343,
    GeometryFromText('POINT(10.987508 43.786805)', 4326));
...
...
 
UPDATE tracksegments
    SET trkseg = (
        SELECT MakeLine(trkpt) FROM (
            SELECT trkpt FROM trackpoints WHERE idtrkseg = 1 ORDER BY trkpt_time
        ) AS trkseg
    ) WHERE id = 1;

L'istruzione UPDATE è quella che provvede ad aggregare i trackpoint e ad assegnare il valore al campo geometrico di tipo LINESTRING nella tabella tracksegments. Contiene due SELECT annidate, questo è il significato di ciascuna parte:

-- La seguente istruzione restituisce un geometry set:
SELECT trkpt FROM trackpoints WHERE idtrkseg = 1 ORDER BY trkpt_time;
 
-- La seguente istruzione trasforma un geometry set in una LINESTRING:
SELECT MakeLine(trkpt) FROM <geometry set> AS trkseg;
 
-- Con questa istruzione si aggiorna il campo trkseg del record:
UPDATE tracksegments SET trkseg = <geometry> WHERE id = 1;

Per il momento non si utilizza la terza e la quarta dimensione (altezza e tempo) offerta dalle entità geometriche di PostGIS 1.0. Sarebbe ovviamente più sensato rispetto all'uso dei campi elevation e time. Però sembra che QGIS - almeno fino alla versione 0.6 - possa lavorare solo con layer PostGIS a due dimensioni. ARGH!!!! QGIS non funziona con PostGIS 1.0!!! No, il pacchetto Debian qgis 0.6.0-2 ha applicata una patch che lo rende funzionante!

Creare una tabella con i contorni di un paese, il campo che memorizza il dato sarà di tipo POLYGON, ovviamente ogni paese può essere costituito da più poligoni (isole, ecc.). Ogni poligono viene creato da una o più LINESTRING, dove la prima definisce il contorno del poligono e le successive definiscono gli eventuali buchi del poligono stesso (ad esempio un lago, oppure lo stato di San Marino all'interno dell'Italia). Per unire le LINESTRING in un POLYGON si usa la funzione di aggregazione MakePolygon(), nell'esempio che segue i poligoni non hanno buchi, ma solo il contorno. E' indispensabile che ogni LINESTRING sia chiusa, cioè il primo e l'ultimo punto devono coincidere.

CREATE TABLE countries (
  idcountry integer NOT NULL,
  idpolygon integer NOT NULL
);
SELECT AddGeometryColumn('countries', 'poly', 4008, 'POLYGON', 2);
 
INSERT INTO countries (idcountry, idpolygon) VALUES (39, 133);
 
UPDATE countries SET poly = (
  SELECT MakePolygon(trkseg) FROM (
    SELECT trkseg FROM tracksegments WHERE idtrkseg = 133
  ) AS trkseg
) WHERE idcountry = 39 AND idpolygon = 133;

A volte è più comodo creare una VIEW e prendere il dato geometrico da questa invece che lavorare direttamente con le tabelle, ad esempio quando si visualizzano i dati in QGIS. Ecco un esempio di come creare la VIEW e aggiornare la tabella geometry_columns in modo appropriato. Ovviamente il nome del campo contenente il dato geometrico, la dimensione, lo SRID e il tipo devono essere conosciuti.

CREATE VIEW my_view AS SELECT oid, nam, na2, the_geom FROM vmap0 WHERE nam LIKE 'SICILIA';
INSERT INTO geometry_columns VALUES('', 'public', 'my_view', 'the_geom', 2, 4326, 'MULTIPOLYGON');

Nota E' necessario che compaia il campo oid altrimenti QGIS non riesce a fare la query. Per poter fare il DROP oppure il CREATE OR REPLACE della query è necessario chiudere QGIS, altrimenti la VIEW risulta in uso e non può essere modificata.

Utilissima questa mail: Using Views as a source of geometric data sull'uso di una VIEW come sorgente dei dati per MapServer.

Può capitare che una query effettuata su PostGIS da QGIS oppure da MapServer risulti estremamente lenta rispetto ad una analoga SELECT effettuata da riga di comando. Tenere presente che QGIS e MapServer attingono dati dichiarando dei BINARY CURSOR invece di utilizzare delle semplici SELECT e quindi il planner di Postgres si comporta in modo diverso. Vedere in proposito questa mail. Purtroppo il comando EXPLAIN ANALYZE per valutare le performance della query non funziona sui BINARY CURSOR.

Le entità geometriche supportate sono quegli oggetti descritti dalle specifiche Simple Features for SQL del consorzio OpenGIS Consortium (OGC). Le specifiche OGC prevedono solo due dimensioni, mentre PostGIS supporta anche la terza (altezza). La quarta dimensione (tempo) pare ancora non supportata (versione PostGIS versione 0.9.1). Attenzione comunque anche ad usare la terza dimensione, pare che QGIS 0.6 non riesca ad usare layer PostGIS con tre dimensioni.

POINT(X Y)
A Point is a 0-dimensional geometry and represents a single location in coordinate space. A Point has a x-coordinate value and a y-coordinate value.

MULTIPOINT(X Y, X Y, ...)
A MultiPoint is a 0 dimensional geometric collection. The elements of a MultiPoint are restricted to Points. The points are not connected or ordered.

LINESTRING(X Y, X Y, ...)
A LineString is a Curve with linear interpolation between points. Each consecutive pair of points defines a line segment.

MULTILINESTRING((X Y, X Y, ...), ...)
A MultiLineString is a collection whose elements are LineStrings.

POLYGON((X Y, X Y, X Y, ...), (X Y, X Y, X Y, ...), ...))
A Polygon is a planar Surface, defined by 1 exterior boundary and 0 or more interior boundaries. Each interior boundary defines a hole in the Polygon. Polygons are topologically closed. The boundary of a Polygon consists of a set of LinearRings (simple and closed Linestrings) that make up its exterior and interior boundaries. No two rings in the boundary cross, the rings in the boundary of a Polygon may intersect at a Point but only as a tangent…

MULTIPOLYGON(...)
A MultiPoligon is a collection whose elements are Polygon.

GEOMETRYCOLLECTION(POINT(X Y),LINESTRING(X Y, X Y, ...), ...)
A GeometryCollection is a collection whose elements others Geometry objects.