PostGIS - Tipps und Tricks
Siehe auch:
Inhaltsverzeichnis
- 1 Tutorial zum Erstellen einer PostGIS-Datenbank
- 2 Koordinatenreferenzsystem-Angaben (SRID)
- 3 Shapefiles importieren (shp2pgsql)
- 4 PostGIS über http-Tool mit Google Earth verknüpfen
- 5 Von PostGIS direkt nach Google Earth
- 6 Benutzerdefinierte Funktionen (Stored Procedures)
- 7 PostGIS optimieren
- 8 Weitere Beispiele
Tutorial zum Erstellen einer PostGIS-Datenbank
Quelle: BostonGIS.
Erläuterungen: "shell>" ist der Prompt einer DOS- oder Linux-Shell. "gisdb=#" ist der Prompt von psql, eingeloggt als gisdb-User.
- PostGIS-Versionen kontrollieren
- PostgreSQL-Version: SELECT version();
- PostGIS-Version:
- SELECT postgis_full_version(); -- (bei aktiver Datenbank)
- sonst template kontrollieren?
- PostGIS-Datenbank erzeugen (Name gisdb)
- Zur Vorbereitung: User 'gisdb' erzeugen:
shell> createuser -a -D gisdb
- Datenbank 'gisdb' mit Template 'template_postgis' und User 'gisdb' erzeugen:
shell> createdb gisdb -O gisdb -T template_postgis -U postgres
- Session und Grants festlegen:
shell> psql -U postgres template1 shell> \c gisdb You are now connected to database "gisdb". gisdb=# GRANT ALL ON TABLE geometry_columns TO postgres,public; gisdb=# GRANT ALL ON TABLE spatial_ref_sys TO postgres,public;
- DB-Projektion (default SRID=-1) vormerken:
gisdb=# select srid, srtext, proj4text from spatial_ref_sys where srtext ILIKE '%Massachusetts%';
- Download vornehmen
- Downloaden: ftp://data.massgis.state.ma.us/pub/shape/state/towns.exe
- Auspacken und überprüfen:
shell> unzip towns.exe shell> ogrinfo -so -al towns_poly.shp
- Shapefile in PostGIS laden
Shapefile nach SQL konvertieren: (-S erzeugt 'single' Polygone (was shapefiles sind) und keine Multipolygone):
shell> shp2pgsql -s 26986 -S TOWNS_POLY towns >towns_poly.sql
- SQL-Daten in PostGIS laden
- ('2>towns_psql_error.txt' leitet Error-Meldungen in eine Datei um):
shell> psql -d gisdb -h localhost -U gisdb -f towns_poly.sql >towns_psql_log.txt 2>towns_psql_error.txt
Beim Import können u.a. Topologie-Fehler zum Vorschein kommen ("Error: geometry contains non-closed rings"). Solche Geometrien werden nicht importiert.
- PostGIS-Index erzeugen und Installation überprüfen
PostGIS-Index erzeugen:
shell> psql -U gisdb gisdb=# CREATE INDEX idx_towns_the_geom ON towns USING gist(the_geom); gisdb=# CREATE INDEX idx_towns_town ON towns USING btree(town);
PostGIS-DB-Installation überprüfen:
gisdb=# SELECT extent(the_geom) FROM towns WHERE town = 'BOSTON'; gisdb=# SELECT area(geomunion(the_geom)) FROM towns WHERE town = 'BOSTON'; gisdb=# SELECT town,astext(transform(the_geom,4326)) FROM towns WHERE town = 'BOSTON'; gisdb=# SELECT count(town) FROM towns WHERE town = 'BOSTON'; gisdb=# SELECT st_askml(geomunion(transform(the_geom,4326))) AS the_geom FROM towns WHERE town = 'BOSTON';
- PostGIS-DB/-Tabelle löschen
- PostGIS-DB:
shell> dropdb gisdb
- Tabelle löschen inkl. aller Einträge in den "Systemtabellen":
abhängigkeiten in den 'Ref-System'-Tabellen?
Koordinatenreferenzsystem-Angaben (SRID)
Der SRID (projection identifier) wird an drei Orten verwaltet:
- In der geometry column: select SRID(the_geom) from mytable limit 5;
- As a constraint on the table for the geometry column: \d mytable
- In the geometry_columns system table: select * from geometry_columns;
Mit der Methode UpdateGeometrySRID([<schema_name>], <table_name>, <column_name>, <srid>) kann man die Kolonne aktualisieren.
Das wohl bekannteste Koordinatenreferensystem ist wohl WGS 84 (long/lat), das von GPS und KML geprägt ist und den Identifier 'EPSG:4326' hat.
Liste von SRID/SRS/CRS: http://www.spatialreference.org/
Shapefiles importieren (shp2pgsql)
Man beachte beim Import von Shapefiles das .prj File. Ohne weitere Angaben wird der SRID auf -1 gesetzt. Siehe SRID (Quelle).
Shapefile-Daten in Datenbank geo1 laden:
% shp2pgsql -s 21781 -I -W ISO-8859-1 gemeinden public.gemeinden > gemeinden.sql % psql -d geo1 -U geo1 -f gemeinden.sql
Eine Alternative zu shp2pgsql ist ogr2ogr von OGR.
PostGIS über http-Tool mit Google Earth verknüpfen
- mit Java, GML und Google Maps - Rendering roads on Google Maps using Java and PostGIS.
- thetimoneygroup.net - On-the-Fly Spatial Analysis With PostGIS and Google Earth.
- (mit Oracle XML-DB).
Von PostGIS direkt nach Google Earth
Man starte psql (Beispiel mit Box um Victoria, BC, Kanada):
-- Set output to unaligned \a -- Show only tuples \t -- Dump query to a file \o my_kml_file.kml -- Run your query SELECT askml('my_line', 'description', 'SRID=3005;LINESTRING( 1190000 390000, 1200000 390000, 1200000 380000, 1190000 380000, 1190000 390000 )'::geometry); -- Flush output your file and close the output stream \o
Dann Doppelklick auf die KML-Datei und Google Earth startet!
Hinweis: Das KML muss ev. mit KML-Kopf und -Fuss ergänzt werden.
Beispiel für einen Ausschnitt der Schweiz:
SELECT AsKML('the_geom', 'SRID=21781;LINESTRING(480000 300000, 835000 300000, 835000 70000, 480000 70000, 480000 300000)'::geometry );
Beispiel für gisdb:
SELECT ST_AsKML(geomunion(transform(the_geom,4326))) as the_geom from towns where town='BOSTON';
Benutzerdefinierte Funktionen (Stored Procedures)
Clean sliver polygons
Given a polygon table that has many small areas and holes. How to remove "small" areas and holes (smaller than a given area in m2)? Remark: Sounds similar like the CLEAN sliver polygons command in ArcGIS.
CREATE OR REPLACE FUNCTION Filter_Rings(geometry,float) RETURNS geometry AS $$ SELECT ST_MakePolygon(c.outer_ring, d.inner_rings) as final_geom FROM (/* Get outer ring of polygon */ SELECT ST_ExteriorRing(b.the_geom) as outer_ring FROM (SELECT (ST_DumpRings($1)).geom As the_geom, path(ST_DumpRings($1)) as path) b WHERE b.path[1] = 0 /* ie the outer ring */ ) c, (/* Get all inner rings > a particular area */ SELECT ST_Accum(ST_ExteriorRing(b.the_geom)) as inner_rings FROM (SELECT (ST_DumpRings($1)).geom As the_geom, path(ST_DumpRings($1)) as path) b WHERE b.path[1] > 0 /* ie not the outer ring */ AND ST_Area(b.the_geom) > $2 ) d $$ LANGUAGE 'sql' IMMUTABLE;
Usage example:
% SELECT ST_AsText( Filter_Rings( ST_PolyFromText( 'POLYGON((10 10,10 20,20 20,20 10,10 10),(0 0,0 1,1 1,1 0,0 0),(5 5,5 7,7 7,7 5,5 5))' ) ,1::float ) ); % "POLYGON((10 10,10 20,20 20,20 10,10 10),(5 5,5 7,7 7,7 5,5 5))"
Shorter alternative to Filter_Rings:
CREATE OR REPLACE FUNCTION Filter_Rings2(geometry,float) RETURNS geometry AS $$ SELECT ST_BuildArea(ST_Collect(a.geom)) as final_geom FROM ST_DumpRings($1) AS a WHERE a.path[1] = 0 OR (a.path[1] > 0 AND ST_Area(a.geom) > $2) $$ LANGUAGE 'sql' IMMUTABLE;
With following restrictions: Possibly slower and squashes 3D geometries to 2D. (Source: http://postgis.refractions.net/pipermail/postgis-users/2009-January/022325.html)
PostGIS optimieren
Siehe PostgreSQL - Tipps und Tricks.