PostGIS - Tipps und Tricks

2012-02-15T09:15:14Z

Zpang: /* PostGIS und Rasterdaten */

Siehe auch:
* [[PostGIS]], [[PostGIS-Beispiele]] und [[PostGIS Snippets]]
* [[GISpunkt-Seminar_PostGIS]]
* [[PostgreSQL]], [[PostgreSQL - Tipps und Tricks]]
* Weitere Tipps und Tricks findet man im [http://spatialdbadvisor.com/postgis_tips_tricks SpatialDBAdvisor Blog für PostGIS]

== Tutorial zum Erstellen einer PostGIS-Datenbank ==

Quelle: [http://www.bostongis.com/?content_name=postgis_tut01 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
oder
shell> createuser -s -D -U postgres gisdb

* Datenbank 'gisdb' mit Template 'template_postgis' und User 'gisdb' erzeugen:
shell> createdb -O gisdb -T template_postgis -U postgres gisdb description

* Session und Grants festlegen:
shell> psql -U postgres template1
template1=# \c gisdb
You are now connected to database "gisdb".
gisdb=# GRANT ALL ON TABLE geometry_columns TO postgres,public;
GRANT
gisdb=# GRANT ALL ON TABLE spatial_ref_sys TO postgres,public;
GRANT

* 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
oder
shell> shp2pgsql -s 4326 -W latin1 Kantone_WGS84 kantone >kantone.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
oder
shell> psql -U postgres -d gisdb2 -f kantone.sql

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 ([http://thesteve0.wordpress.com/2007/10/08/updating-projections-and-postgis/ 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 laden ==

* Mittels (mit PostGIS ausgeliefertem) Kommandozeilen-Tool 'shp2pgsql'.
* Mittels INSERT und WKT oder Konstruktoren (vgl. Snippet "Creating Geometry (...) Types" unten).

== PostGIS-Daten darstellen ==

=== Übersicht ===
* Lokal, mittels Desktop GIS:
** [[QGIS]]
** [[OpenJUMP]]
** GDV Spatial Commander [http://www.gdv.com/down/scommander.php]
** gvSIG [http://www.gvsig.gva.es/index.php?id=gvsig&L=2] gvSIG Wiki [http://webmapping.info/mwgvsig/index.php?title=Hauptseite]
** etc.
* über Webapplikationen ([[WMS]]):
** [[GeoServer]]
** [[UMN MapServer]]
** [[ArcGIS]] Server

=== PostGIS und KML ===

PostGIS kennt die Funktion ST_AsKML(), die eine Geometrie KML-konform ausgibt. Allerdings werden nur diejenigen KML-Elemente der Geometrie selbst ausgegeben. Kopf- und Fusszeile des KML-Dokuments muss man selber hinzufügen. Dazu erzeugen wir nun eine Funktion askmldoc, die ein gültiges KML-Dokument ausgibt:

CREATE OR REPLACE FUNCTION askmldoc(name text, description text, the_geom geometry) RETURNS text AS $$
DECLARE
result text;
BEGIN
result := '<?xml version="1.0" encoding="UTF-8"?>' || E'\n' ||
'<kml xmlns="http://www.opengis.net/kml/2.2">' || E'\n' ||
'<Document>' || E'\n' ||
'<name>' || name || '</name>' || E'\n' ||
'<description>' || description || '</description>' || E'\n\n' ||
'<Style id="defaultStyle">' || E'\n' ||
' <LineStyle>' || E'\n' ||
' <color>ff00ff00</color>' || E'\n' ||
' <width>1</width>' || E'\n' ||
' </LineStyle>' || E'\n' ||
' <PolyStyle>' || E'\n' ||
' <color>5f00ff00</color>' || E'\n' ||
' </PolyStyle>' || E'\n' ||
'</Style>' || E'\n\n' ||
'<Placemark>' || E'\n' ||
'<styleUrl>#defaultStyle</styleUrl>' || E'\n';
result := result || ST_AsKML(the_geom) || E'\n';
result := result ||
'</Placemark>' || E'\n\n' ||
'</Document>' || E'\n' ||
'</kml>';
RETURN result;
END;
$$ LANGUAGE plpgsql;

Exportieren Sie damit die Geodaten in der Tabelle, z.B. mit folgendem Befehl auf der System-Kommandozeile (cmd):

% psql -A -t -d gisdb -c "SELECT askmldoc('MyTable', 'A Comment...', the_geom) FROM mytable WHERE gid=1;" -o mykmlfile.kml

=== PostGIS und Google Earth ===

PostGIS über http-Tool mit Google Earth verknüpfen:
* [http://cse-mjmcl.cse.bris.ac.uk/blog/2006/05/27/1148737394599.html mit Java, GML und Google Maps] - Rendering roads on Google Maps using Java and PostGIS.
* [http://www.thetimoneygroup.net/foss/ thetimoneygroup.net] - On-the-Fly Spatial Analysis With PostGIS and Google Earth.
* ([http://www.oracle.com/technology/pub/articles/schalk-googlemaps.html 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';

== ST_Relate ==

Documentation:
* Strobl paper
* PostGIS documentation ch04
* ...

Syntax:
ST_Relate(geom_a, geom_b, '<<param>>')
...

Example:
* ST_Contains(polygon_a, polygon_b) => ST_Relate(.., .., '****1****')
* ST_Crosses
* ST_Intersects
* ST_Touches

== PostGIS und Rasterdaten ==

Ab PostGIS Version 2.0 (geplant 1. Quartal 2011) gibt es neu einen Raster-Datentyp.

Siehe [[PostGIS Raster]].

== 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.

See:
* Function cleanGeometry(): http://trac.osgeo.org/postgis/wiki/UsersWikiCleanPolygons
* Topology function: http://trac.osgeo.org/postgis/wiki/UsersWikiPostgisTopology

== PostGIS optimieren ==

Siehe auch:
* [[PostgreSQL - Tipps und Tricks#PostgreSQL_optimieren]] (tuning).
* [http://www.foss4g2007.org/presentations/view.php?abstract_id=117 Tips for the PostGIS Power User ]

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:
* [[PostGIS Snippets]]
* [[PostGIS-Beispiele]] mit British-Columbia-Daten
* [http://trac.osgeo.org/postgis/wiki/UsersWikiMain PostGIS-Wiki]

== Fragen und Antworten ==

;Was ist der Unterschied von geometry und geography type?: Geography type besteht aus geografischen Koordinaten (z.B. lat/lon) und wird für grossflächige Daten verwendet, die z.B. "länderübergreifend" sind.

;Wie wird ein Geometrie-Attribut erzeugt?: als separates Statement, gleich nach CREATE TABLE... Mit AddGeometryColumn. Ab PostGIS 1.5 ist es endlich möglich, den konkreten Geometrie-Typ direkt im CREATE STATEMENT anzugeben.

;Was bedeutet SRID und wie kann man das ändern?: SRID bedeutet Spatial Reference ID und ist ein Fremdschlüssel zum [[Koordinatenreferenzsystem]]. SRID information is stored *both* in the geometry_columns table (srid column), and inside each geometry itself. Ist opbligatorischer Parameter u.a. von AddGeometryColumn oder beim Geometrie-Erzeugen. Kann mit UpdateGeometrySRID geändert werden.

;Was bedeutet der Error "addgeometrycolumn(...) does not exist"?: Die PostGIS-Funktionen sind nicht geladen. Man installiere die PostGIS-Erweiterung (falls nicht schon geschehen) und erzeuge eine neue PostGIS-Datenbank (siehe "Wie erzeugt man eine neue PostGIS-Datenbank?").

;Wie erzeugt man eine neue PostGIS-Datenbank?: 1. Falls noch keine 'normale' (=Nicht-GIS-) Datenbank vorhanden ist, "create a new database and choose the template_postgis as your template". 2. Falls eine 'normale' Datenbank mit Tabellen etc. bereits existiert, "spatially enable your database by running the following two SQL files: a) Postgresql\8.2\share\contrib\lwpostgis.sql and b) Postgresql\8.2\share\contrib\spatial_ref_sys.sql. lwpostgis.sql benötigt plpgsql language support installed (Syntax: createlang plpgsql yourexistingdb).

;Was ist das für eine Zahl, die ST_distance() zurückgibt?: Siehe [http://postgis.refractions.net/support/wiki/ PostGIS Wiki].

;Was ist der Unterschied zwischen ST_Overlaps, ST_Crosses, ST_Intersects und ST_Contains?: Siehe [http://www.postgis.org/documentation/manual-1.5SVN/ST_Overlaps.html PostGIS-Doku. zu ST_Overlaps].

;Was ist der Unterschied zwischen ST_Extent und ST_Expand?: Siehe [http://www.bostongis.com/postgis_extent_expand_buffer_distance.snippet].

;Was für Mass-Einheiten (units of measurement) benützen die PostGIS-Funktionen?: Siehe [http://postgis.refractions.net/support/wiki/ PostGIS Wiki].

;Wie erhält man Distanz-Masse in Meter?: Siehe [http://postgis.refractions.net/support/wiki/ PostGIS Wiki].

;Was sind die Grenzen der 3D-Typen und -Funktionen in PostGIS?: (1) 3D volumetric objects are not supported. (2) 3D non-volume are supported partly -- e.g. a 2D polygon in 3 d space, a line, point in 3d space. (3) Spatial relationships however only consider the spatial component of the x, y plane.

;MySQL hat auch eine Spatial Extension und ist erst noch bei Internet-Service-Providern verbreiteter; was ist der Unterschied zu PostGIS?: Drei Gründe: (1) Operationen werden nur mittels Bounding-Box gemacht ([http://dev.mysql.com/doc/refman/5.1/en/functions-that-test-spatial-relations hips-between-geometries.html], (2) es gibt keine Unterstützung für Koordinatensysteme, (3) die Spatial Extension ist praktisch eine One-Man-Show.

== Snippets ==

=== Get area of a polygon which has SRID EPSG:4326 (WGS84/GPS) ===

All area functions are squared by definition.
But a sq degree varies in its area depending on latitude.
Transform your data using a SRID referring to an equal-area (or near enough to
equal area for your purposes) projection in the query for area:

-- Get area from a geometry given the geometry is in 4326.
% SELECT ST_Area2d(ST_transform(geometry, <SRID>));

=== Creating Geometry and Geography Types ===
Allowed and non-functioning constructors for creating geometry types (Varianten für geometry types):
/*0*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)')) -- '''Preferred simplest text form without SRID'''
/*1*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)', 4326)) -- '''Preferred for text form with SRID'''
/*2*/ SELECT ST_AsEWKT(ST_GeomFromText('SRID=4326;POINT(-71.06 42.28)')) -- Note: srid 'inline'
/*3*/ SELECT ST_AsEWKT(ST_GeomFromEWKT('SRID=4326;POINT(-71.06 42.28)')) -- Alternative to ST_GeomFromText()
/*4*/ SELECT ST_AsEWKT('SRID=4326;POINT(-71.06 42.28)'::geometry) -- With cast; prefer *1*/ST_GeomFromText
/*5*/ SELECT ST_AsEWKT(ST_SetSRID('POINT(-71.06 42.28)'::geometry,4326)) -- With cast; prefer *1*/ST_GeomFromText
/*6*/ SELECT ST_AsEWKT(ST_MakePoint(-71.06, 42.28, 4326)) -- '''Preferred symbolic form''' (Hint: returns WKT, not EWKT)

Allowed and non-functioning constructors for creating geography types (Varianten für geography types):
/*1*/ -- SELECT ST_AsEWKT(ST_GeogFromText('POINT(-71.06 42.28)', 4326)) -- ERROR: unknown ST_GeogFromText()
/*2a*/ -- SELECT ST_AsEWKT(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)')) -- ERROR: no ST_AsEWKT: why not?
/*2b*/ SELECT ST_AsText(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)')) -- Preferred with srid 'inline'; (Hint: ST_AsEWKT doesn't work)
/*3*/ -- SELECT ST_AsText(ST_GeogFromEWKT('SRID=4326;POINT(-71.06 42.28)')) -- ERROR: no ST_GeogFromEWKT()
/*4*/ SELECT ST_AsText('SRID=4326;POINT(-71.06 42.28)'::geography) -- With cast; prefer *1*; (Hint: ST_AsEWKT doesn't work)
/*5*/ --- SELECT ST_AsText(ST_SetSRID('POINT(-71.06 42.28)'::geography,4326)) -- ERROR: ok SetSRID w. geography unnecessary
/*6*/ -- Equivalent to ST_MakePoint() for geography types missing. -- .

=== 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.

=== Nearest Neighbors ===

Find n (e.g. 5) nearest neighbors for given point coordinate (a classic query).
From '[postgis-users] Find n Nearest Neighbors for given Point using PostGIS?'.

CREATE TABLE geoname ( -- (from geonames.org)
geoname varchar(255),
geonameid integer,
geom geometry
);
-- PRIMARY KEY UNIQUE BTREE index for geonameid

UPDATE geoname SET geom = ST_SetSRID(ST_Point(longitude,latitude) 4326);
...
CREATE INDEX geom_index ON geoname USING GIST (geom);)
CLUSTER geom_index ON geoname;

Possible solution (inspired by [http://www.bostongis.com/PrinterFriendly.aspx?content_name=postgis_nearest_neighbor]):

SELECT
start.asciiname,
ende.asciiname,
distance_sphere(start.geom,ende.geom) as distance
FROM geoname AS start, geoname AS ende
WHERE start.geonameid = 2950159
AND start.geonameid <> ende.geonameid
AND ST_DWithin(start.geom, ende.geom, 30)
ORDER BY distance LIMIT 5;

There will be a KNNGIS index probably in PostgreSQL 9.1 which makes such queries even faster. It returns the found tupels in decreasing order until LIMIT without computing and sorting whole set before returning it.

The way currently is to write a stored procedure that expands the search if it fails to get the number of results. So make a set returning function with the body something like this (pseudo code) and you should get good performance:

radius := 0.01; -- assuming degrees
loop
select into cnt count(*) from mytable
where expand(mypnt, radius) && the_geom limit 5;
if found and cnt = 5 or radius > maxradius then
for rr in select * from mytable
where expand(mypnt, radius) && the_geom limit 5
loop
return rr;
end loop;
return;
else
radius := radius * 2;
end if;
end loop;

=== 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.

=== Convert multipolygon to polygon ===
Q: How to convert multipolygon to polygon geometries

A: If multipolygon only have one member, then st_geometryn(geom, 1) will do it. Otherwise, look at function st_dump().

=== Getting all polygon's coordinates ===

Q. Given a geometry attribute with polygons, how can one get all the polygon's coordinates (as a tupel set)? E.g. either

1 POINT(1 2)
2 POINT(2 2)
3 POINT(3 2)
...
... or even better:
X Y
1 1 2
2 2 2
3 3 2
...

1. Maybe a "SELECT ST_AsText(polygon_column) FROM Table" is enough

Output: POLYGON( (1 1, 2 2, 3 3) )

2.
SELECT
DISTINCT ST_X(ST_POINTN(foo.lines, foo.index)) AS x,
ST_Y(ST_POINTN(foo.lines, foo.index)) AS y
FROM (SELECT thelines.lines, generate_series(1,thelines.n) AS index
FROM (SELECT polys.geo AS lines, ST_NUMPOINTS(polys.geo) AS n
FROM (SELECT ST_BOUNDARY((ST_DUMP( the_geom )).geom) AS geo
FROM "YOURTABLE" WHERE YOURQUERY
) AS polys
) AS thelines
) AS foo;

SELECT DISTINCT
-- foo.index,
ST_X(ST_POINTN(foo.lines, foo.index)) AS x,
ST_Y(ST_POINTN(foo.lines, foo.index)) AS y
FROM (SELECT thelines.lines, generate_series(1,thelines.n) AS index
FROM (SELECT polys.geo AS lines, ST_NUMPOINTS(polys.geo) AS n
FROM (SELECT ST_BOUNDARY((ST_DUMP( the_geom )).geom) AS geo
FROM "gemeinden" WHERE gid=4
) AS polys
) AS thelines
) AS foo

Should work for multipolygons and polygons.

3. You can import a plpgsql implementation of this until it's implemented in C (http://trac.osgeo.org/postgis/attachment/ticket/76/my_st_dump_points.2.sql).
[[Media:]]

Usage:
SELECT ST_AsText( (My_ST_DumpPoints( poly )).geom )
FROM (
SELECT 'POLYGON((0 0, 1 1, 1 0, 0 0))'::geometry AS poly
) AS foo;

ST_AsText ------------
POINT(0 0)
POINT(1 1)
POINT(1 0)
POINT(0 0)
(4 rows)

=== Split a polygon to N linestrings ===

Source: http://postgis.refractions.net/pipermail/postgis-users/2010-January/025818.html

CREATE TEMP TABLE mypolygontable AS
SELECT
'MULTIPOLYGON(((0 0, 0 4, 4 4, 4 0, 0 0),
(1 1, 1 2, 2 2, 2 1, 1 1)),
((5 5, 5 6, 6 6, 6 5, 5 5))
)'::geometry geom;

-- make line segments from every startpoint and endpoint
SELECT ST_AsText( ST_MakeLine(sp,ep) )
FROM
-- extract the endpoints for every 2-point line segment for each
(SELECT
ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
ST_PointN(geom, generate_series(2, ST_NPoints(geom) )) as ep
FROM
-- extract the individual linestrings
(SELECT (ST_Dump(ST_Boundary(geom))).geom
FROM polygons
) AS linestrings
) AS segments;

st_astext
---------------------
LINESTRING(0 0,0 4)
LINESTRING(0 4,4 4)
LINESTRING(4 4,4 0)
LINESTRING(4 0,0 0)
LINESTRING(1 1,1 2)
LINESTRING(1 2,2 2)
LINESTRING(2 2,2 1)
LINESTRING(2 1,1 1)
LINESTRING(5 5,5 6)
LINESTRING(5 6,6 6)
LINESTRING(6 6,6 5)
LINESTRING(6 5,5 5)
(12 rows)

More general solution:
* http://www.spatialdbadvisor.com/postgis_tips_tricks/128/exploding-a-linestring-or-polygon-into-individual-vectors-in-postgis/

=== Aggregate linestrings ===

Q. Trying to aggregate linestrings together based on their attributes and the fact that they are touching each other.

A. There is a ST_Linemerge function that takes a collection of linestrings as an argument and merges them together into a multilinestring.

With a table "thetable" having an "attrib" field and a "geom" field, one could write the query like this :

SELECT attrib, (st_dump(merged_geom)).geom
FROM (
SELECT attrib, ST_Linemerge(ST_Collect(geom)) AS merged_geom
FROM thetable
GROUP BY attrib
) AS subq;

=== Polygons that don't overlap ===

Q. I have a table with a MULTIPOLYGON field. Some of these records have self intersections and other problems that make ST_IsValid return false. The polygons overlap a lot and I wanted to generate another table that will be the union of all polygons. The table looks like: id, the_geom. What would be the best way to union all of the geometries into a new table where there is only POLYGONS that do not overlap?

A. See SQL:

SELECT id, ST_Union(ST_Accum(CASE WHEN ST_IsValid(a.the_geom) = 't' THEN a.the_geom ELSE st_buffer(a.the_geom,0) end) AS the_geom
FROM (SELECT id, the_geom
FROM mytable
WHERE the_geom IS NOT NULL
AND geometrytype(the_geom) != 'GEOMETRYCOLLECTION'
UNION ALL
SELECT gid, (ST_Dump(the_geom)).geom AS the_geom
FROM mytable
WHERE the_geom IS NOT NULL
AND geometrytype(the_geom) = 'GEOMETRYCOLLECTION'
) AS a;

The split based on GEOMETRYCOLLECTION is that when I first tested the STR indexed ST_Union I found that performance was affected if the union set includes an mpoly.

[[Kategorie:PostGIS]] [[Kategorie:HowTo]]

Geometa Lab OST - Benutzerbeiträge [de-ch]

PostGIS - Tipps und Tricks