Sem PostGIS Roger: Unterschied zwischen den Versionen

Version vom 1. November 2007, 23:19 Uhr

Wiki-Seite zum Selbststudium im GISpunkt-Seminar PostGIS.

Daten: Amtliche Vermessung Kanton Basel-Stadt

Verwendete Daten

Als Übungsdaten werden zwei Datensätze der amtlichen Vermessung des Kantons Basel-Stadt verwendet. Der eine Datensatz enthält die Rohrleitungen (INTERLIS-Topic), der zweite die Bodenbedeckung (INTERLIS-Topic), aus der im Laufe der Übung die Tankanlagen (Bodenbedeckungsart 17) extrahiert werden.

Die Datensätze liegen als Shapefiles Bodenbedeckung_merged.shp und Rohrleitungen.shp vor.

Metadaten-Record

Bemerkung: Die Angabe beim DC-Element modified bezieht sich auf das Datum des Datenbezugs.

Rohrleitung

 dc:title           = Rohrleitung
 dct:abstract       = Bestandteil des Datenmodells des Bundes. Beschreibt die geometrische Form und 
                      Lage der Rohrleitungen, die von übergeordneter Bedeutung für die 
                      Eigentumsverhältnisse sind.
 dc:format          = INTERLIS, Shapefile
 dct:spatial        = name=Kanton Basel-Stadt; northlimit=47.60200; southlimit==47.51858; 
                      eastlimit=7.69527; westlimit=7.55281
 dct:modified       = 2007-10-19
 dc:publisher       = Grundbuch- und Vermessungsamt Basel-Stadt
 dc:language        = de
 dc:rights          = Abgabebedingungen: Berechtigter Interessennachweis. Die Daten können aufgrund 
                      einer Benützungsbewilligung bezogen werden. Entsprechende Bezugsformulare sind 
                      auf www.gva.bs.ch erhältlich. Für die Verwendung der Daten gelten die allgemeinen 
                      Bedingungen für die Benützung von Geodaten des Grundbuch- und Vermessungsamtes 
                      Basel-Stadt.

Tankanlagen (Art der Bodenbedeckung)

 dc:title           = Bodenbedeckung
 dct:abstract       = Bestandteil des Datenmodells des Bundes. Beschreibt die geometrische Form und Lage 
                      der Flächenarten mit einer Mindestfläche von in der Regel 100 m2 (Gebäude, befestigt, 
                      humusiert, bestockt, Gewässer, vegetationslos, etc.).
 dc:format          = INTERLIS, Shapefile
 dct:spatial        = name=Kanton Basel-Stadt; northlimit=47.60200; southlimit==47.51858; 
                      eastlimit=7.69527; westlimit=7.55281
 dct:modified       = 2007-10-19
 dc:publisher       = Grundbuch- und Vermessungsamt Basel-Stadt
 dc:language        = de
 dc:rights          = Abgabebedingungen: Die Daten können aufgrund einer Benützungsbewilligung bezogen 
                      werden. Entsprechende Bezugsformulare sind auf www.gva.bs.ch erhältlich. Für die 
                      Verwendung der Daten gelten die allgemeinen Bedingungen für die Benützung von 
                      Geodaten des Grundbuch- und Vermessungsamtes Basel-Stadt.

Laden des Datensatzes

Die vorliegenden Shapefiles Bodenbedeckung_merged.shp und Rohrleitungen.shp wurden mit shp2pgsql.exe in .sql-Files konvertiert, diese wiederum mit psql.exe in die PostGIS-Datenbank geladen.

Die Shapefiles wurden mit zwei Batchfiles in die Datenbank importiert: Das erste Batchfile übergibt die zu importierenden Files als Parameter an das zweite Batchfile:

 echo ==============================
 echo  Shapefile in PostGIS-DB laden
 echo ==============================
 echo.
 echo  Übergibt Shapefilenamen (ohne
 echo  Endung) an Hauptbatchfile.
 echo.
  
 set shp2pg_main="_shp2pg_main.bat"
  
 call %shp2pg_main% Bodenbedeckung_merged
 call %shp2pg_main% Rohrleitungen
  
 pause

Das zweite Batchfile lädt das so übergebene Shapefile in die Datenbank:

 echo ==============================
 echo  Shapefile in PostGIS-DB laden
 echo ==============================
 echo.
 echo  Erwartet Shapefilename (ohne
 echo  Endung) als Parameter (%1).
 echo.
  
 set shp2pgsql_exe="C:\Programme\PostgreSQL\8.2\bin\shp2pgsql.exe"
 set psql_exe="C:\Programme\PostgreSQL\8.2\bin\psql.exe"
  
 set srid="21781"
 set schema="public"
 set server="localhost"
 set dbname="postgis-seminar"
 set dbuser="postgres"
  
 echo --------------------
 echo  SQL-File erstellen
 echo --------------------
 %shp2pgsql_exe% -s %srid% -I %1 %schema%.%1 > %1.sql
  
 echo -------------------
 echo  Daten in DB laden
 echo -------------------
 %psql_exe% -h %server% -d %dbname% -U %dbuser% -f %1.sql

Unter Verwendung der UNIX-Pipes können Konvertierung und Laden in die Datenbank auch in einem Schritt durchgeführt werden (s. S. 62 in [1]).

Queries

1. Aufbereiten des Übungsdatensatzes Tankanlagen

Tankanlagen sind als Teil der Bodenbedeckung in der Tabelle bodenbedeckung_merged vorhanden. Für eine einfachere Handhabung werden Sie in eine separate Tabelle tankanlagen ausgelagert. (Die Rohrleitungen sind bereits in Tabelle rohrleitungen aus dem Ursprungsdatensatz rohrleitungen.shp vorhanden.)

 -- SQL-Script: Tabelle_tankanlagen_erstellen_und_fuellen.sql 

 -- Neue Tabelle tankanlagen mit Primary Key erstellen
 CREATE TABLE tankanlagen
 (
   gid integer PRIMARY KEY,
   art_txt text,
   art text
 );
 
 -- Geometriespalte the_geom hinzufügen
 SELECT AddGeometryColumn('public', 'tankanlagen', 'the_geom', 21781, 'MULTIPOLYGON', 2);
 
 -- Daten von Tabelle bodenbedeckung_merged in Tabelle tankanlagen importieren (Tankanlagen sind art = 17)
 INSERT INTO tankanlagen
 (
   SELECT gid, art_txt, art, the_geom 
   FROM bodenbedeckung_merged
   WHERE art = '17'
 );
 
 -- Räumlichen Index für Tabelle tankanlagen erstellen
 CREATE INDEX tankanlagen_the_geom_gist
   ON tankanlagen 
   USING gist (the_geom GIST_GEOMETRY_OPS);
 
 -- Statistiken für optimierte Abfragen anlegen
 VACUUM ANALYZE tankanlagen;

2. Umprojizieren (CH03->WGS84)

Die Daten sind ursprünglich im Koordinatensystem CH03/LV03 (EPSG-Code oder SRID 21781) vorhanden. Für den anschliessenden KML-Export müssen sie jedoch in geografische Koordinaten (WGS84, EPSG-Code 4326) transformiert werden. Zu Übungszwecken wird auf zwei Arten projiziert, einmal anhand einer zweiten Geometriespalte (Rohrleitungen), und einmal on-the-fly (Tankanlagen; s. Abschn. "Als KML speichern").

2.1 Rohrleitungen (mit zusätzlicher Geometriespalte)

Neben der schon vorhandenen Geometriespalte the_geom mit SRID 21781, wird für den an-schliessenden KLM-Export eine zweite Geometriespalte mit dem SRID 4326 hinzugefügt. Doppeltes Ablegen von Geometrien kann bei häufigen Abfragen Sinn machen.

 -- Zweite Geometriespalte mit EPSG-Code 4326 (WGS84) hinzufügen
 SELECT AddGeometryColumn('public', 'rohrleitungen', 'the_geom_4326', 4326, 'MULTILINESTRING', 2)
 
 -- Neue Geometriespalte the_geom_4326 füllen
 UPDATE rohrleitungen
 SET the_geom 4326 = transform(the_geom, 4326);
 
 -- Räumlichen Index für Tabelle rohrleitungen (the_geom_4326) erstellen
 CREATE INDEX rohrleitungen_the_geom_4326_gist
   ON rohrleitungen 
   USING gist (the_geom_4326 GIST_GEOMETRY_OPS);
 
 -- Statistiken für optimierte Abfragen anlegen
 VACUUM ANALYZE rohrleitungen;

2.2 Tankanlagen (on-the-fly)

Die Daten in der Tabelle tankanlagen werden bei der KML-Generierung on-the-fly in WGS84 umprojiziert (s. Abschn. "Als KML speichern").

3. Als KML speichern

Die KML-Dateien wurden mit der Funktion AsKML(text, text, geometrie) gemäss Anleitung auf [2] erstellt. Mit dieser Funktion wird die Ausgabe KML-konform formatiert.

Wichtig ist, dass beim SELECT AsKML(...) die Geometrien zusammengefasst werden, z.b. mit geomunion(), sonst wird der Inhalt mehrerer KML-Dateien (eine für jeweils eine Geometre) in eine einzige geschrieben (d.h. mehrere <?xml>-Tags).

3.1 Tankanlagen

Die Tankanlagen wurden mit folgendem SELECT on-the-fly in WGS84-Koordinaten transformiert und in KLM-Form zurückgegeben:

 SELECT AsKML(
   'Tankanlagen', 
   'Amtliche Vermessung Basel-Stadt, Bodenbedeckung Art 17', 
   geomunion( Transform( the_geom, 4326 ) ) )
FROM tankanlagen;

Ergebnis: tankanlagen.kml

3.2 Rohrleitungen

Die Geometrien sind nicht nur in der Spalte the_geom in Schweizer Landeskoordinaten ab-gelegt (EPSG-Code 21781), sondern auch in einer zweiten Spalte the_geom_4326 in geografischen Koordinaten, die für den KML-Export benötigt werden. Hier ist somit keine on-the-fly-Transformation nötig.

 SELECT AsKML( 
   'Rohrleitungen', 
   'Amtliche Vermessung Basel-Stadt', 
   geomunion( the_geom_4326 ) )
 FROM rohrleitungen;

Ergebnis: rohrleitungen.kml

4. Ausschnitt (Perimeter, definiert mittels BBox)

4.1 Bounding Box generieren

Im Folgenden wird aus den Eckkoordinaten 610750/269000 und 612600/271100 eine Bounding Box generiert.

 SELECT ST_MakeBox3D(
   ST_GeomFromText( 'POINT( 610750 269000 )', 21781 ),
   ST_GeomFromText( 'POINT( 612600 271100 )', 21781 ) 
 );

Ausgabe der Query:

 BOX3D( 610750 269000, 612600 271100 )

4.2 Geometrien mit Bounding Box ausschneiden

Die in der oben erstellten Bounding Box enthaltenen Objekte aus der Tabelle Tankanlagen sollen in einer separaten Tabelle tankanlagen_aoi abgelegt werden. &&-Operator....

 -- Tabelle tankanlagen_aoi aus BBox-Abfrage erstellen und füllen.
 CREATE TABLE "public"."tankanlagen_aoi" AS
 (
   SELECT *
   FROM tankanlagen 
   WHERE the_geom
         && SetSRID( 'BOX3D( 610750 269000, 612600 271100 )'::box3d, 21781 )
  );
 
 -- Primärschlüssel hinzufügen
 ALTER TABLE tankanlagen_aoi ADD PRIMARY KEY( gid );
 
 -- Räumlichen Index kreieren
 CREATE INDEX "tankanlagen_aoi_the_geom_gist" 
   ON "public"."tankanlagen_aoi" using gist ("the_geom" gist_geometry_ops);
 
 - Statistiken für optimierte Abfragen anlegen
 VACUUM analyze tankanlagen_aoi;

4.3 Bounding Box erstellen und darin enthaltene Objekte wählen (in einer einzigen Query)

 CREATE TABLE "public"."tankanlagen_aoi" AS
 (
   SELECT *
   FROM tankanlagen 
   WHERE the_geom
     && SetSRID( 
     (
       SELECT ST_MakeBox3D(
         ST_GeomFromText( 'POINT( 610750 269000 )', 21781 ),
         ST_GeomFromText( 'POINT( 612600 271100 )', 21781 ) 
       )
      )::box3d, 21781 )
 );

Die Funktion SetSRID() ordnet der BBox das entsprechende SRID zu; andernfalls würde es defaultmässig -1 betragen.

Nun kann, wie oben beschrieben, noch ein Primärschlüssel und ein räumlicher Index hinzugefügt werden.

5. Buffer rund um einen Ausschnitt

Um die weiter oben erstellte Bounding Box soll ein Buffer berechnet werden. Die auf diese Weise ausgegebene Bounding Box (BOX3D(…)) soll gleich in die Buffer-Query integriert werden. Um die Bounding Box vom Type box3d in eine Geometrie umzuwandeln, stand früher die Funktion GeometryFromText() zur Verfügung; in neueren PostGIS-Versionen ist diese nicht mehr implementiert, stattdessen kann SetSRID() verwendet werden[3]. Als Buffer wurden 1000m gewählt.

 SELECT AsText(
   ST_Buffer( 
     SetSRID( 'BOX3D( 610750 269000, 612600 271100 )'::box3d, 
       21781 ),
     1000
   )
 );

Ergebnis: Ein Polygon, das ein Rechteck mit "abgerundeten Ecken" repräsentiert.

6. Spatial Join (Overlay, Flächenverschnitt)

Es soll ermittelt werden, welche Teilstücke der Rohrleitungen wenigstens teilweise innerhalb eines Buffers von 50 Metern um die Tankanlagen liegen, diese Buffer also überlappen:

 SELECT r.gid, r.the_geom
   FROM rohrleitungen r, tankanlagen t
   WHERE r.the_geom && ( SELECT ST_Buffer( t.the_geom, 50 ) )
   GROUP BY r.gid, r.the_geom;

Ergebnis: rohrleitungen_nahe.kml (Erstellung: zuerst wurde aus obenstehendem SELECT analog dem Abschnitt "Geometrien mit Bounding Box ausschneiden" die Tabelle rohrleitungen_nahe kreiert und gefüllt, die dann analog dem Abschnitt "Als KML speichern" als .kml ausgegeben wurde).