Sem PostGIS Juliane: Unterschied zwischen den Versionen

Aktuelle Version vom 10. Februar 2009, 19:32 Uhr

Wiki-Seite zum Selbststudium im GISpunkt-Seminar PostGIS. Zurück zu Sem PostGIS Selbststudium.

Daten: OpenStreetMap von Zürich und Umgebung

Verwendete Daten

Als Übungsdaten werden OpenStreetMap-Daten der Schweiz verwendet. Diese stehen kostenfrei unter http://download.geofabrik.de/osm/europe/ zum Download zur Verfügung. Der Datensatz switzerland.shp.zip enthält Gebäude (buildings.shp), Bodenbedeckung (natural.shp), Bahnen (railways.shp), Strassen (roads.shp), Gewässernetz (waterways.shp) und Points of Interest (points.shp). Aus dem Datensatz points.shp werden zur beispielhaften Bearbeitung der Aufgabe die Bushaltestellen extrahiert (vgl. Query 1).

Zudem werden die Schweizer Landesgrenze (ponet.e00) und die Grenzen der Siedlungsgebiete (pppoly.e00) des "Digital Chart of the World" (DCW)-Datensatzes bezogen. Diese Daten können ebenfalls kostenfrei heruntergeladen werden.

Metadaten-Records

OpenStreetMap

 dc:title           = OpenStreetMap
                      Buildings, Natural, Railways, Roads, Waterways, Points
                      
 dct:abstract       = Die Dateien werden regelmässig mit dem Tool Osmosis aus dem so genannten "Planet File" ausgeschnitten, das 
                      einen Gesamt-Abzug der "OpenStreetMap"-Datenbank darstellt. Solche Auszüge werden von der Geofabrik GmbH 
                      erstellt, aktualisiert und als Shapefiles bereitgestellt. Folgende Elemente sind in den einzelnen Shapefiles 
                      gespeichert und teilweise bereits vollständig attributiert:
                       
                      buildings.shp
                      Gebäude, Öffentliche Gebäude, etc.
                      
                      natural.shp
                      Wälder, Parkanlagen, Seen, etc.
                      
                      railways.shp
                      Eisenbahnlinien, Tramlinien, Bahnhöfe, U-Bahnen, Schmalspurbahnen, etc. 
                      
                      roads.shp 
                      Fusswege, Strassen (1. Klasse, 2. Klasse usw.), Radwege, etc.
                      
                      waterways.shp
                      Kanäle, Bäche, Flüsse, Wehre, etc.
                       
                      points.shp
                      Busstationen, Almhütten, Museen, Cafés, Hotels, Kindergärten, Universitäten, Tierparks, etc.
                      
 dc:format          = ESRI Shapefile
                        
 dct:spatial        = buildings.shp
                      northlimit=47.734204; southlimit=46.006147; eastlimit=10.409632; westlimit=6.037668
                      
                      natural.shp
                      northlimit=47.808295; southlimit=45.905356; eastlimit=10.481553; westlimit=5.957567 
                   
                      railways.shp
                      northlimit=47.745569; southlimit=45.832066; eastlimit=10.287146; westlimit=5.993640
  
                      roads.shp
                      northlimit=47.808603; southlimit=45.821394; eastlimit=10.488137; westlimit=5.965783
                  
                      waterways.shp
                      northlimit=47.783465; southlimit=45.904293; eastlimit=10.487938; westlimit=5.978676
                   
                      points.shp
                      northlimit=47.792457; southlimit=45.832330; eastlimit=10.453054; westlimit=6.000755
                  
 dct:modified       = 02-Feb-2009 07:41
 dc:publisher       = Geofabrik GmbH Karlsruhe 
 dc:language        = en
 dc:rights          = Alle Daten aus dem OSM-Projekt unterliegen der Creative Commons Attribution Share-Alike-Lizenz.

Digital Chart of the World

Bemerkung: Die Angabe beim DC-Element modified bezieht sich auf das Datum des Datenbezugs. Die Original DCW-Daten stammen vermutlich aus dem Jahr 1993.

 dc:title           = Digital Chart of the World
                      Ponet, Pppoly
            
 dct:abstract       = Verschiedene Kartenebenen zahlreicher Länder können im Arc/INFO Exportformat von ESRI's "Digital Chart of the
                      World"-Datensatz heruntergeladen werden:
              
                      ponet.e00
                      Landesgrenze der Schweiz
        
                      pppoly.e00
                      Siedlungsgebiete der Schweiz
     
 dc:format          = Arc/INFO coverage 
          
 dct:spatial        = ponet.e00
                      northlimit=47.810918; southlimit=45.796475; eastlimit=10.552661; westlimit=5.912654
            
                      pppoly.e00
                      northlimit=47.721973; southlimit=45.982600; eastlimit=9.668165; westlimit=6.040542
                 
 dct:modified       = 04-Feb-2009 11:28 
 dc:publisher       = maproom, Pennsylvania State University 
 dc:language        = en
 dc:rights          = Die ursprünglichen Daten der DCW wurden von der US Defence Mapping Agency erstellt und sind über die
                      National Imagery and Mapping Agency erhältlich.

Laden der Datensätze

Die Coverages ponet.e00 und pppoly.e00 werden zuerst in QuantumGIS geladen und als Shapefiles exportiert (rechte Maustaste auf entsprechenden Layer > Save as shapefile). Zu Übungszwecken für eine anschliessende Umprojektion (vgl. Query 2) wird als Coordinate reference System CH1903/LV03 (EPSG: 21781) ausgewählt.

Alle Shapefiles werden mit shp2pgsql.exe in sql-Files konvertiert.

-- Statement 
shp2pgsql -I -W codierung shapefilename tablename > outputsql.sql

Konvertieren der "OpenStreetMap"-Daten:

shp2pgsql –s 4326 –I –W ISO-8859-15 buildings gebaeude > buildings.sql
shp2pgsql –s 4326 –I –W ISO-8859-15 natural flaechen > natural.sql
shp2pgsql –s 4326 –I –W ISO-8859-15 points points > points.sql
shp2pgsql –s 4326 –I –W ISO-8859-15 railways bahnen > railways.sql
shp2pgsql –s 4326 –I –W ISO-8859-15 roads strassen > roads.sql
shp2pgsql –s 4326 –I –W ISO-8859-15 waterways gewaesser > waterways.sql

Konvertieren der Daten der "Digital Chart of the World":

shp2pgsql –s 21781 –I –W ISO-8859-15 ponet grenze > ponet.sql
shp2pgsql –s 21781 –I –W ISO-8859-15 pppoly siedlungen > pppoly.sql

Anschliessend werden die sql-Files über Terminal mit psql.exe in die PostGIS-Datenbank auf dem Server geladen.

-- Statement
psql –U rolename –d datenbankname –h host –p port –f filename.sql

Queries

Benutzte Versionen

• PostGIS Version 1.1.4 (PC) mit PostgreSQL Version 8.1.4
• PostGIS Version 1.3.3-1 (Mac) mit PostgreSQL Version 8.3.5 für kml-Generierung

Query 1: Datensatz mit mindestens einem Thema

Der Datensatz points.shp beinhaltet zahlreiche Points of Interest (siehe Metadaten-Records). Für die weitere Bearbeitung werden alle Bushaltestellen selektiert und in eine neue Tabelle gespeichert.

-- Neue Tabelle busstationen anlegen
CREATE TABLE busstationen
(
  gid serial NOT NULL,
  osm_id bigint,
  "name" character varying(32),
  "type" character varying(16),
  the_geom geometry,
  CONSTRAINT busstationen_pkey PRIMARY KEY (gid),
  CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2),
  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'POINT'::text OR the_geom IS NULL),
  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 4326)
)
WITHOUT OIDS;
ALTER TABLE busstationen OWNER TO swai;
  
-- Index kreieren
CREATE INDEX busstationen_the_geom_gist
  ON busstationen
  USING gist
  (the_geom);
 
-- Busstationen aus Tabelle points selektieren und in die zuvor angelegte Tabelle busstationen einfügen
INSERT INTO busstationen
 (
  SELECT gid, osm_id, "name", "type", the_geom 
  FROM points
  WHERE "type" = 'bus_station'
 );

Query 2: Umprojizieren (CH03->WGS84)

Die "Digital Chart of the World"-Daten liegen im Koordinatensystem CH03/LV03 (EPSG: 21781) vor. Für den anschliessenden KML-Export müssen sie jedoch in das geografische Koordinatensystem WGS84 (EPSG: 4326) transformiert werden. Die Umprojektion wird mit Hilfe einer zweiten Geometriespalte (wgs84) durchgeführt.

-- Geometriespalte hinzufügen
SELECT AddGeometryColumn('public','ponet','wgs84',4326,'MULTILINESTRING', 2);
SELECT AddGeometryColumn('public','pppoly','wgs84',4326,'MULTILINESTRING', 2);

-- Daten umprojizieren
UPDATE ponet SET wgs84 = transform(the_geom, 4326);
UPDATE pppoly SET wgs84 = transform(the_geom, 4326);

Die "OpenStreetMap"-Daten müssen nicht umprojiziert werden. Sie liegen bereits im WGS84 Koordinatensystem vor.

Query 3: KML speichern

Die KML-Dateien werden mit der Funktion AsKML(text, text, geometrie) [1] erstellt. Mit dieser Funktion wird die Ausgabe KML-konform formatiert.

-- KML ponet
SELECT askml(
  'schweiz',
  'landesgrenze',
  geomunion(transform(wgs84,4326)))
FROM ponet;

Die Konvertierung alle anderen Datensätze dauert sehr lange. Aus diesem Grund werden alle weiteren KML’s aus den selektierten BBox Dateien (vgl. Query 4) erzeugt. Im Folgenden beispielhaft die Speicherung als KML für die Tabellen pppoly_zh und buildings_zh:

-- KML aus BBox pppoly_zh
SELECT askml(
  'zuerich',
  'siedlungsflaeche',
  geomunion(transform(wgs84,4326)))
FROM pppoly_zh;

 
-- KML aus BBox buildings_zh
SELECT askml(
  'zuerich',
  'gebaeude',
  geomunion(transform(the_geom,4326)))
FROM buildings_zh;

Query 4: Ausschnitt anfragen (Perimeter, definiert mittels BBox)

Als Bounding Box-Umgebung wird die Region Zürich mit folgenden Eckkoordinaten definiert: 8.45 47.45, 8.6 47.35. Die in der Bounding Box enthaltenen Objekte der verschiedenen Tabellen werden jeweils in eine separate Tabelle tabellenname_zh abgelegt. Anschliessend wird noch ein Primärschlüssel und ein räumlicher Index hinzugefügt (hier beispielhaft für die Tabellen pppoly und buildings).

-- Tabelle pppoly_zh aus BBox-Abfrage erstellen und füllen
CREATE TABLE pppoly_zh AS
(
  SELECT *
  FROM pppoly 
  WHERE wgs84 && SetSRID( 'BOX3D( 8.45 47.45, 8.6 47.35 )'::box3d, 4326 )
);

-- Primärschlüssel hinzufügen
ALTER TABLE pppoly_zh ADD PRIMARY KEY( gid );

-- Räumlichen Index kreieren
CREATE INDEX pppoly_zh_wgs84_gist 
  ON pppoly_zh 
  using gist (wgs84);
 
 
-- Tabelle buildings_zh aus BBox-Abfrage erstellen und füllen
CREATE TABLE buildings_zh AS
(
  SELECT *
  FROM buildings 
  WHERE the_geom && SetSRID( 'BOX3D( 8.45 47.45, 8.6 47.35 )'::box3d, 4326 )
);

-- Primärschlüssel hinzufügen
ALTER TABLE buildings_zh ADD PRIMARY KEY( gid );

-- Räumlichen Index kreieren
CREATE INDEX buildings_zh_the_geom_gist 
  ON buildings_zh 
  using gist (the_geom);

Query 5: Buffer rund um ein Linien-Element

Es wird ein Buffer von ca. 400 m um die Limmat gelegt und anschliessend ebenfalls als KML-Datei gespeichert.

SELECT Buffer(the_geom,0.005) as buffer, 
* FROM waterways_zh 
WHERE "name"='Limmat';

-- Buffer als KML speichern
SELECT askml(
  'zuerich',
  'buffer limmat',
  geomunion(transform(Buffer(the_geom,0.005),4326))) 
FROM waterways_zh 
WHERE "name"='Limmat';

Bemerkung: Umrechnung ° in km

Längengradabstand a in km = cos(φ) × 2 π r / 360

Für Zürich φ = 47°40’ und r = 6370 km -> 0.1° ca. 75 km

Query 6: Spatial Join (Overlay, Flächenverschnitt) zweier Feature-Klassen vom Typ Fläche

Es soll ein Spatial Join aller Gebäude von buildings_zh, die innerhalb der Polygone von natural_zh liegen, erstellt werden.

SELECT b."name", b.the_geom, b.gid, n."type"
FROM natural_zh n, buildings_zh b 
WHERE Contains(n.the_geom, b.the_geom);

Query 7: Eigene Query

Es soll ermittelt werden von welchem Gebäude innerhalb natural_zh (vgl. Query 6) die Distanz zur nächsten Bushaltstelle am kürzesten ist.

Bemerkung: Umrechnung ° in km

Längengradabstand a in km = cos(φ) × 2 π r / 360

Für Zürich φ = 47°40’ und r = 6370 km -> 0.1° ca. 75 km

SELECT s.gid, s."name", b."name", b.osm_id,
  MIN (distance(s.the_geom, b.the_geom)*750) AS min_dist
FROM spatial_join s, busstationen_zh b
GROUP BY s.gid, s."name", b."name", b.osm_id
ORDER BY min_dist ASC
LIMIT 1;