PostGIS - Tipps und Tricks: Unterschied zwischen den Versionen

Aus Geoinformation HSR
Wechseln zu: Navigation, Suche
(→‎Weitere Beispiele)
(→‎Select nearest point on polyline)
Zeile 247: Zeile 247:
  
 
Ab PostGIS 1.3.5 (wegen generate_series() )?
 
Ab PostGIS 1.3.5 (wegen generate_series() )?
−
 
+
  SELECT P.gid, P.geom As pclosest
−
SELECT P.gid, P.geom As pclosest
+
  FROM (SELECT  
−
FROM (SELECT  
+
    gid,  
−
  gid,  
+
    ST_PointN(the_geom, generate_series(1, ST_NPoints(the_geom))) AS geom
−
  ST_PointN(the_geom, generate_series(1, ST_NPoints(the_geom))) As geom
+
  FROM point_table) AS P
−
FROM point_table) AS P
+
  WHERE ST_Dwithin(P.geom, ST_SetSRID(ST_MakePoint(lon,lat),somesrid), 0.01)
−
WHERE ST_Dwithin(P.geom, ST_SetSRID(ST_MakePoint(lon,lat),somesrid), 0.01)
 
−
 
 
  
 
=== PostGIS views ===
 
=== PostGIS views ===

Version vom 18. Oktober 2009, 02:56 Uhr

Siehe auch:


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
 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?

Shapefiles in PostGIS 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-Daten darstellen

Übersicht

PostGIS und Google Earth

PostGIS über http-Tool mit Google Earth verknüpfen:

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

PostGIS und Koordinatenreferenzsystem-Angaben

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/

Beim Aufruf von GDAL/OGR wird die GDAL_DATA environment Variable benötigt.

Benutzerdefinierte PostgreSQL/PostGIS-Funktionen

Implementiert als sog. 'Stored Procedures'.

Siehe auch PostGIS-Beispiele.

Eliminate 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)?

Remarks:

  • Similar like the ELIMINATE command in ArcGIS.
  • See also CLEAN topology
 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)

Clean topology

Remarks:

  • Similar to ArcGIS' CLEAN.
  • See also Eliminate sliver polygons.
 tbd.

PostGIS optimieren

Siehe auch:

Version:

  • 1.4 and GEOS 3.1 will bring prepared geometries which will make things faster.

Reorder on disk:

  • import the data (into a temporary table)
  • and then sort the data using some sort of spatial key. A centroid should be OK but consider using a larger bucket such as a gridded area.
 > CREATE TABLE gis_roads 
     AS SELECT * FROM tmp_gis_roads 
     ORDER BY CENTROID(the_geom);

Cluster on disk: 

 > CREATE INDEX idx_mytable_the_ON mytable 
     USING GIST(the_geom);
 > CLUSTER mytable 
     USING idx_mytable_the_geom;

Don'ts:

  • Don't use just one big diagonal lines. Break it down because one diagonal line has a huge bounding box (which is what the index works on) while 20 equivalent shorter lines have a much smaller combined box area.
  • Don't use this construction: ST_Intersects(Buffer(GeomA, D), GeomB) Use this one: ST_DWithin(GeomA, GeomB, D)

Weitere Beispiele

Siehe auch:

Extent of database

If you want to find out the EXTENT parameters of your database, run this command: 

 # su postgres
 # psql -d <yourdatabase> -c "SELECT extent(geometry) FROM <yourtable>"

Replace <your...> with the relevant database and table names.

Select nearest point on polyline

A user clicks somewhere on a polyline and we want to get the closest point that is part of the representation of the line. You can replace 0.01 with some other tolerance number to compensate for the fact that a person probably won't click right on the line.

Ab PostGIS 1.3.5 (wegen generate_series() )? 

 SELECT P.gid, P.geom As pclosest
 FROM (SELECT 
   gid, 
   ST_PointN(the_geom, generate_series(1, ST_NPoints(the_geom))) AS geom
 FROM point_table) AS P
 WHERE ST_Dwithin(P.geom, ST_SetSRID(ST_MakePoint(lon,lat),somesrid), 0.01)

PostGIS views

Q: I like to polygonize linestrings in a view to see the result in qgis. But qgis need gid. How can I create the gid column after an aggregate function in a view?

A1. In postgresql 8.3 you can force oid creation if WITH OIDS is specified when the table is created: 

 CREATE VIEW myview AS SELECT oid AS gid, the_geom FROM mytable;

A2. Create a sequence and use then nextval(sequence_name) to get a primary key.

A3. Grab the key from the underlying geometry table, & use it for a key column in your aggregate view, eg. "create view v1 as select min(key) as key" ...

Standard direction for polylines

Q: Is there a standard direction for polylines?

A: There isn't a standard, but you can force your polygons into an orientation with ST_ForceRHR, which forces a clockwise orientation for exterior rings and anti-clockwise for interior rings.

Portion of river geometry along administrative boundary

Q: I have a rivers linestring table and a political boundaries (counties) table polygon. I am wishing to extract the portion of the river geometry that lies along a county border. Of course, these two layers don't have the same individual point precision, so I suspect some buffering is necessary.

A: See following SQL query: 

 intersection(
      buffer(exteriorring(geometryn(n.geom,1)),0.02),
      r.geom
 )

So turn the county polygon into an exterior ring and then buffer it out a bit. Then see what portions of the river intersect that buffered polygon.

Von «https://giswiki.hsr.ch/index.php?title=PostGIS_-_Tipps_und_Tricks&oldid=24369»