Stefan: hat „GISpunkt-Seminar PostGIS Uebungen“ nach „GISpunkt-Seminar PostGIS-Uebungen“ verschoben

2011-09-13T11:26:21Z

hat „GISpunkt-Seminar PostGIS Uebungen“ nach „GISpunkt-Seminar PostGIS-Uebungen“ verschoben

← Nächstältere Version

Version vom 13. September 2011, 11:26 Uhr

Stefan: Die Seite wurde neu angelegt: „Zurück zu GISpunkt-Seminar PostGIS. Hier einige Abfragen die mit dem Übungsdatensatz gemacht werden können. Die Übungsdaten können Sie von [http://www.c…“

2011-09-12T13:19:29Z

Die Seite wurde neu angelegt: „Zurück zu GISpunkt-Seminar PostGIS. Hier einige Abfragen die mit dem Übungsdatensatz gemacht werden können. Die Übungsdaten können Sie von [http://www.c…“

Neue Seite

Zurück zu [[GISpunkt-Seminar PostGIS]].

Hier einige Abfragen die mit dem Übungsdatensatz gemacht werden können. Die Übungsdaten können Sie von [http://www.carto.net/postgis/pgch_uebungsdaten.zip http://www.carto.net/papers/pgch_uebungsdaten.zip] herunterladen. Sie benötigen die shapefiles "gemeinden" und "fluesse" und das SQL File staedte_ch.sql

=== Neue räumliche Tabelle erstellen und bevölkern ===
<pre>
-- Eine neue räumliche Tabelle "seen" soll erstellt und alle Seen von der Tabelle
-- "gemeinden" eingespielt werden. Danach soll ein räumlicher Index auf der neuen Tabelle
-- "seen" erstellt werden und die Seen aus der Tabelle "gemeinden" gelöscht werden

-- neue Tabelle kantone mit Primary Key erstellen
CREATE TABLE seen
(
gid integer PRIMARY KEY,
name text
);
-- räumliche Geometriespalte dazufügen
SELECT AddGeometryColumn('public', 'seen', 'the_geom', 21781, 'MULTIPOLYGON', 2);
-- Daten von Gemeinde hereinspielen wo kt = 0 (sind Seen), sollten 22 Datensätze sein
INSERT INTO seen
SELECT gid, name, the_geom FROM gemeinden WHERE kt = 0 AND gmde >= 9000;
-- Seen und liechtensteinische Gemeinden aus gemeindetabelle löschen
DELETE FROM gemeinden WHERE kt = 0;
-- neuen räumlichen index erstellen bei der tabelle seen
CREATE INDEX seen_the_geom_gist
ON seen USING gist (the_geom GIST_GEOMETRY_OPS);
</pre>

Resultat: neue Tabelle Seen mit Index und 36 Datensätzen

=== Flächenberechnungen ===
<pre>
-- Selektiere alle Gemeinden, die eine Fläche grösser als 100 km2 haben
SELECT name, kt, Round(ST_Area(the_geom) / 1000000) As flaeche FROM gemeinden
WHERE (ST_Area(the_geom) / 1000000) > 100
ORDER BY flaeche DESC;
</pre>

Resultat: Bagnes, 23, 283; Davos, 18, 255; Zermatt, 23, 243; ... insgesamt 48 Datensätze (wenn Seen in vorherigem Schritt gelöscht)

=== Südlichste Geometrie ermitteln ===
<pre>
-- Selektiere die südlichste Gemeinde des Kantons Graubünden
-- kt id 18 ist Graubünden

SELECT name FROM gemeinden WHERE kt = 18 ORDER BY ST_YMIN(ST_ENVELOPE(the_geom)) ASC LIMIT 1;
</pre>

Resultat: Roveredo (GR)

=== Koordinatentransformation ===
<pre>
-- Die Koordinate von Bern (Tabelle staedte_schweiz), urspr. Landeskoordinate,
-- als Lat/Lon, WGS84 ausgeben
SELECT AsText(Transform(the_geom,4326)) FROM staedte_schweiz WHERE "FULL_NAME" = 'Bern';
</pre>
Resultat: POINT(7.46666670742301 46.9166667102846)

<pre>
-- Eine Koordinate von einem fremden Koordinatensystem in WGS 84 transformieren:
SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(484344.9842453949 1480036.13293874)',32648),4326))
</pre>
Resultat: POINT(104.855408409113 13.3878890725056) -- WGS 84 Grad (EPSG:4326)

=== Distanzberechnung ===
<pre>
-- Sphärische Distanz (spheroid) zwischen Bern und Zürich, basierend auf Bessel 1841 Ellipsoid
SELECT ST_distance_spheroid(
(SELECT Transform(the_geom,4326) FROM staedte_schweiz WHERE "FULL_NAME" = 'Bern'),
(SELECT Transform(the_geom,4326) FROM staedte_schweiz WHERE "FULL_NAME" = 'Zürich'),
'SPHEROID["Bessel 1841",6377397.155,299.1528128]');
</pre>

Resultat: 96139.7469334825

=== Kombinierte Distanzabfragen ===
<pre>
-- Selektiere Orte (Tabelle staedte_schweiz) und zugehörige Seen (Tabelle Seen) die innerhalb eines Radius
-- von 75 km von Bern sind und ein Ort nicht mehr als 5km von einem See entfernt ist (Tabelle Seen)
-- Achtung: kartesische Distanz, nicht ellipsoidische!
SELECT sta."FULL_NAME", se.name FROM staedte_schweiz sta, seen se
WHERE ST_DWithin(sta.the_geom,(SELECT the_geom FROM staedte_schweiz WHERE "FULL_NAME" = 'Bern'),75000) AND
ST_DWithin(sta.the_geom,se.the_geom,5000) ORDER BY sta."FULL_NAME" ASC, se.name ASC;
</pre>

Resultat: Adligenswil, Vierwaldstättersee; Aecherli, Lungerer See; Aecherlis, Sarnersee; ...

=== Selektiere angrenzende Flächen mit ST_Touches() ===
<pre>
-- Selektiere alle Gemeinden die an den Zürichsee grenzen

SELECT name FROM gemeinden
WHERE ST_TOUCHES(the_geom,(SELECT the_geom FROM seen WHERE name = 'Zürichsee')) ORDER BY name;
</pre>

Resultat: 25 Records; Altendorf, Erlenbach (ZH, Freienbach, ...

=== Crosses(), Selektion von Flächen die von Linien geschnitten werden ===
<pre>
-- Es sollen sämtliche Gemeinden abgefragt durch die der Fluss Emme fliesst

SELECT g.name FROM gemeinden g, fluesse f
WHERE f.name = 'Emme' AND f.the_geom && g.the_geom AND ST_Crosses(f.the_geom,g.the_geom)
ORDER by name ASC;
</pre>

=== Geomunion() und Crosses() Abfrage (MULTIPOLYGON und MULTILINESTRING) ===
<pre>
-- Prüfe ob die kombinierten Einzel-Linienelemente des Rheinflusses die Gemeindegeometrie von Basel durchlaufen
-- Resultat ergibt "t" (true) oder "f" (false)
SELECT ST_Crosses(
(SELECT the_geom FROM gemeinden WHERE name = 'Basel'),
(SELECT geomunion(the_geom) FROM fluesse WHERE name = 'Rhein' GROUP BY name));
</pre>

Resultat: t

=== Flächen zusammenfassen mit Geomunion() ===
<pre>
-- Eine neue räumliche Tabelle "kantone" soll erstellt und Gemeinden aus der Tabelle
-- "gemeinden" zu Kantonen zusammengefasst eingespielt werden. Danach soll ein räumlicher Index auf der neuen Tabelle
-- "kantone" erstellt werden

-- neue Tabelle kantone mit Primary Key erstellen
CREATE TABLE kantone
(
gid integer PRIMARY KEY,
name text
);
-- räumliche Geometriespalte dazufügen
SELECT AddGeometryColumn('public', 'kantone', 'the_geom', 21781, 'MULTIPOLYGON', 2);
-- Kantonsgeometrien von Gemeindegeometrien zusammenfassen mit Geomunion, gruppiert nach Kantons-ID "kt"
-- ST_Multi() erzwingt Multigeometrien wegen dem constraint, es sollten 26 Datensätze resultieren
INSERT INTO kantone
SELECT kt, '', ST_Multi(Geomunion(the_geom)) FROM gemeinden GROUP BY kt;
-- neuen räumlichen index erstellen bei der tabelle seen
CREATE INDEX kantone_the_geom_gist
ON kantone USING gist (the_geom GIST_GEOMETRY_OPS);
-- danach manuell die Kantonsnamen ergänzen, etwa in QGIS, OpenJump oder GDV Spatial Commander
</pre>

Resultat: neue Tabelle mit 26 Kantonsgeometrien und räumlicher Index

=== Point in Polygon mit ST_Within() ===
<pre>
-- Es soll abgefragt werden wieviele Städte/Orte in jedem Kanton liegen
-- Voraussetzung ist, dass Sie Ihren Kantonen Namen vergeben haben

SELECT k.name, COUNT(*) FROM
staedte_schweiz s, kantone k
WHERE ST_WITHIN(s.the_geom,k.the_geom) GROUP BY k.name ORDER BY k.name;
</pre>

Resultat: Aargau, 146; Appenzell Ausserrhoden, 20; Appenzell Innerrhoden, 8; ...

=== Buffer() und Point in Polygon ST_Within() ===
<pre>
-- Selektierte Orte die innerhalb eines Buffers von 10 km des Flusses Emme liegen

SELECT DISTINCT "FULL_NAME" FROM staedte_schweiz
WHERE ST_Within(
the_geom,
(SELECT ST_Buffer(the_geom,10000) FROM fluesse WHERE name = 'Emme'))
ORDER BY "FULL_NAME" ASC;
</pre>

Resultat: 144 Orte; Aedermannsdor, Aelgäu, Aeschau, Affoltern, ...

Dazu gäbe es eine schnellere Alternative:
<pre>
SELECT DISTINCT staedte_schweiz.full_name FROM staedte_schweiz, fluesse
WHERE ST_DWithin(staedte_schweiz.the_geom, fluesse.the_geom, 10000)
AND fluesse.name = 'Emme'
ORDER BY full_name ASC;
</pre>

=== Anzahl Teilgeometrien ermitteln ===
<pre>
-- Selektiere alle Kantone die aus 4 oder mehr Teilgeometrien bestehen

SELECT name, ST_NumGeometries(the_geom) AS NrGeom
FROM kantone
WHERE ST_NumGeometries(the_geom) >= 4
ORDER BY name ASC;
</pre>

Resultat: Bern, 4; Fribourg, 4; Solothurn, 4

=== ST_DWithin() ===

<pre>
-- Given a set of Point objects, filter out the Points from this set
-- that lie on a straight line within a certain amount of tolerance

SELECT p.* FROM yourPoints AS p, yourLine AS line
WHERE ST_DWithin(p.the_geom, line.the_geom, yourTolerance)
</pre>

=== Räumliche Beziehungen mit ST_Relate() ===
<pre>
-- Es sollen sämtliche Seen abgefragt werden die im Kanton Bern sind oder an den Kanton Bern grenzen
-- Voraussetzung ist, dass Sie Ihren Kantonen Namen vergeben haben

-- Selektiere alle Seen die innerhalb des Kantons Bern liegen
SELECT s.name, 'innerhalb' AS Relation FROM seen s, kantone k
WHERE k.name = 'Bern' AND ST_RELATE(s.the_geom,k.the_geom) = 'FF2F1F212' ORDER BY k.name;

-- Selektiere alle Seen die an der Grenze zum Kanton Bern liegen
SELECT s.name, 'angrenzend' AS Relation FROM seen s, kantone k
WHERE k.name = 'Bern' AND ST_RELATE(s.the_geom,k.the_geom) = 'FF2F11212' ORDER BY k.name;
</pre>

Resultat 1: Brienzersee, innerhalb; Thunersee, innerhalb

Resultat 2: Lac de Neuchâtel, angrenzend; Bielersee / Lac de Bienne, angrenzend

=== Suche Punkt entlang einer Linie mit Linear Referencing ===
<pre>
-- Selektiere die Gemeinde die sich 20 km vom Ursprung der Emme befindet
-- Achtung, die Emme wurde offenbar von der Mündung bis zum Ursprung digitalisiert
-- die Geometrie muss daher umgedreht werden
-- ST_GeometryN() wird verwendet, da es sich um einen Multilinestring handelt

SELECT name FROM gemeinden WHERE
ST_WITHIN((SELECT ST_line_interpolate_point(ST_Reverse(ST_GeometryN(the_geom,1)),(20000 / ST_LENGTH(the_geom)))
FROM fluesse WHERE name = 'Emme'),the_geom);
</pre>

Resultat: Eggiwil

=== Cut Shapefiles At A Border ===

Use ST_Intersection in conjunction with ST_Intersects. ST_Intersection will
clip the geometries so what is left is what falls inside both geometries.

So something like:

<pre>
SELECT i.gid, ST_Intersection(country.the_geom, i.the_geom)
FROM country
INNER JOIN someothertable i ON ST_Intersects(country.the_geom, i.the_geom)
</pre>

For your countries that are linestrings, you will need to convert them to
polygons with something like ST_BuildArea, see
http://postgis.refractions.net/documentation/manual-svn/ST_BuildArea.html

=== Set of closest Points to a Location (ST_Distance) ===

Say you have geometry A and geometry B and you want to know what set of points on Geometry A are closest to Geometry B: Which points in ST_Distance are selected? Return distance, startpoint, endpoint.

SELECT ST_Intersection(A, ST_Buffer(B, ST_Distance(A,B) + 0.001) );

Solution: 1) First find the minimum distance. 2) Figure out which point/pointset on the geometry satisfies the minimum distance using buffer and intersection (with minimum distance being the buffer distance).

Then if you really wanted the points -- you would use the imaginary not yet invented ST_DumpPoints should the result not be point or be a geometry collection of some sort.

GISpunkt-Seminar PostGIS-Uebungen - Versionsgeschichte

Stefan: hat „GISpunkt-Seminar PostGIS Uebungen“ nach „GISpunkt-Seminar PostGIS-Uebungen“ verschoben

Stefan: Die Seite wurde neu angelegt: „Zurück zu GISpunkt-Seminar PostGIS. Hier einige Abfragen die mit dem Übungsdatensatz gemacht werden können. Die Übungsdaten können Sie von [http://www.c…“