PostGIS-Beispiele

Für weitere PostGIS-Beispiele siehe auch GISpunkt-Seminar PostGIS.

British Columbia-Daten

Quelle: Paul Ramsey, Refractions Research, 11. Jul. 2007, [1]

Es folgt eine Reihe von SQL-Commandos:

HINWEIS: Die Daten zu dieser Übung (bc_pubs.shp, bc_data.sql) findet man hier (FOSS4G 2007 Workshop 4).

Simple Spatial SQL

 create table points ( pt geometry, name varchar );
 insert into points values ( 'POINT(0 0)', 'Origin' );
 insert into points values ( 'POINT(5 0)', 'X Axis' );
 insert into points values ( 'POINT(0 5)', 'Y Axis' );
 select name, AsText(pt), ST_Distance(pt, 'POINT(5 5)') from points;
 drop table points;

Loading Shape Files

 cd \postgis-workshop\data
 pg_setenv
 shp2pgsql -i -s 3005 bc_pubs.shp bc_pubs > bc_pubs.sql
 notepad bc_pubs.sql
 pg_shpsql
 psql -U postgres -f bc_data.sql -d postgis
   

Creating Spatial Indexes

 create index bc_roads_gidx on bc_roads using gist ( the_geom );
 create index bc_pubs_gidx on bc_pubs using gist ( the_geom );
 create index bc_voting_areas_gidx on bc_voting_areas using gist ( the_geom );
 create index bc_municipality_gidx on bc_municipality using gist ( the_geom );
 create index bc_hospitals_gidx on bc_hospitals using gist ( the_geom );

Using Spatial Indexes

 -- Unindexed query
 select gid, name from bc_roads where _ST_Crosses(the_geom, ST_GeomFromText('LINESTRING(1220446 477473,1220417 477559)', 3005));
 -- Indexed query
 select gid, name from bc_roads where ST_Crosses(the_geom, ST_GeomFromText('LINESTRING(1220446 477473,1220417 477559)', 3005));

Indexes and Query Plans

 -- Unindexed query plan
 explain select gid, name from bc_roads where _ST_Crosses(the_geom, ST_GeomFromText('LINESTRING(1220446 477473,1220417 477559)', 3005));
 -- Indexed query plan
 explain select gid, name from bc_roads where ST_Crosses(the_geom, ST_GeomFromText('LINESTRING(1220446 477473,1220417 477559)', 3005));

PostgreSQL Optimization

  open postgresql.conf

Spatial Analysis in SQL

 -- What is the total length of all roads in the province, in kilometers?
 select sum(ST_Length(the_geom))/1000 as km_roads from bc_roads;

 -- How large is the city of Prince George, in hectares?
 select ST_Area(the_geom)/10000 as hectares from bc_municipality where name = 'PRINCE GEORGE';

 -- What is the largest municipality in the province, by area?
 select name, ST_Area(the_geom)/10000 as hectares from bc_municipality order by hectares desc limit 1;

Data Integrity

 -- How many invalid voting area polygons?
 select count(*) from bc_voting_areas where not ST_IsValid(the_geom);

Distance Queries

 -- Where is the Tabor Arms pub?
 select ST_AsText(the_geom) from bc_pubs where name ilike 'Tabor Arms%';

 -- How many Unity Party supporters live within 2km of Tabor Arms pub?
 select sum(upbc) as unity_voters from bc_voting_areas where ST_DWithin(the_geom, ST_GeomFromText('POINT(1209385 996204)', 3005), 2000);

Spatial Joins

 -- Find all pubs located within 250 meters of a hospital.
 select h.name, p.name from bc_hospitals h, bc_pubs p where ST_DWithin(h.the_geom, p.the_geom, 250);

 -- Summarize the 2000 provincial election results by municipality.
 select m.name, sum(v.ndp) as ndp, sum(v.lib) as liberal, sum(v.gp) as green, sum(v.upbc) as unity, sum(v.vtotal) as total from bc_voting_areas v, bc_municipality m where ST_Intersects(v.the_geom, m.the_geom) group by m.name order by m.name;

Overlays

 -- Clip the voting areas by the Prince George boundary.
 create table pg_voting_areas as select ST_Intersection(v.the_geom, m.the_geom) as intersection_geom, ST_Area(v.the_geom) as va_area, v.*, m.name from bc_voting_areas v, bc_municipality m where ST_Intersects(v.the_geom, m.the_geom) and m.name = 'PRINCE GEORGE';

 -- What is the area of the clipped features?
 select sum(ST_Area(intersection_geom)) from pg_voting_areas;

 -- How does that compare to the clipping polygon? (Should be the same.)
 select ST_Area(the_geom) from bc_municipality where name = 'PRINCE GEORGE';

Coordinate Projection

 -- What is the proj4 definition of srid 3005?
 select proj4text from spatial_ref_sys where srid=3005;

 -- What are the native coordinate of the first road in the roads table?
 select ST_AsText(the_geom) from bc_roads limit 1;

 -- What do those coordinates look like in lon/lat?
 select ST_AsText(ST_Transform(the_geom,4326)) from bc_roads limit 1;

Exercises

 -- What is the perimeter of the municipality of Vancouver?
 select ST_Perimeter(the_geom) from bc_municipality where name = 'VANCOUVER';
 
 -- What is the total area of all voting areas in hectares?
 select sum(ST_Area(the_geom))/10000 as hectares from bc_voting_areas;
 
 -- What is the total area of all voting areas with more than 100 voters in them?
 select sum(ST_Area(the_geom))/10000 as hectares from bc_voting_areas where vtotal > 100;
 
 -- What is the length in kilometers of all roads named Douglas St?
 select sum(ST_Length(the_geom))/1000 as kilometers from bc_roads where name = 'Douglas St';

Advanced Exercises

 -- What is the length in kilometers of Douglas St in Victoria?
 select sum(ST_Length(r.the_geom))/1000 as kilometers from bc_roads r, bc_municipality m where ST_Contains(m.the_geom, r.the_geom) and r.name = 'Douglas St' and m.name = 'VICTORIA';
 
 -- What two pubs have the most Green Party supporters within 500 meters of them?
 select p.name, p.city, sum(v.gp) as greens from bc_pubs p, bc_voting_areas v where ST_DWithin(p.the_geom, v.the_geom, 500) group by p.name, p.city order by greens desc limit 2;

 -- What is the latitude of the most southerly hospital in the province?
 select ST_Y(ST_Transform(the_geom,4326)) as latitude from bc_hospitals order by latitude asc limit 1;

 -- What were the percentage NDP and Liberal vote within the city limits of Prince George in the 2000 provincial election?
 select 100*sum(v.ndp)/sum(v.vtotal) as ndp, 100*sum(v.lib)/sum(v.vtotal) as liberal from bc_voting_areas v, bc_municipality m where ST_Contains(m.the_geom, v.the_geom) and m.name = 'PRINCE GEORGE';

 -- What is the largest voting area polygon that has a hole?
 select gid, id, ST_Area(the_geom) as area from bc_voting_areas where ST_NRings(the_geom) > 1 order by area desc limit 1;

 -- How many NDP voters live within 50 meters of Simcoe St in Victoria?
 select sum(v.ndp) as ndp from bc_voting_areas v, bc_municipality m, bc_roads r where ST_Contains(m.the_geom, r.the_geom) and r.name = 'Simcoe St' and m.name = 'VICTORIA' and ST_DWithin(r.the_geom, v.the_geom, 50);

Weitere

List of districts with given watershed

There are two tables:

1. ) One table consist of 150 watershed polygons
2. ) Second table consista of 500 district polygons.

and a watershed name called "Sussex".

Some of the districts are completely whithin one certain watershed (100%), but some has intersection with more than one watershed. How I can write SQL to know percentage of district intersection area with certain watershed (i.e. fraction of intersected district area divided by area of one watershed). As an output I will have 500 rows consiting name of district and percentage of its area in certain watershed.

1) If you want to list all 500 districts for a given watershed regardless of if any of it lies in watershed, then you would do

SELECT d.district_name, w.ws_name,

 CASE WHEN ST_Intersects(w.the_geom, d.the_geom)
   THEN
     ST_Area(ST_Intersection(w.the_geom,d.the_geom))/ST_Area(d.the_geom)
   ELSE 0 END AS per_distinws
     FROM districts AS d CROSS JOIN water_sheds AS w
 WHERE w.ws_name = 'Sussex'

2) If you only care about the districts that fully or partly full in shed then

 SELECT d.district_name, w.ws_name,
   ST_Area(ST_Intersection(w.the_geom,d.the_geom))/ST_Area(d.the_geom)  AS per_distinws
   FROM districts AS d
   INNER JOIN water_sheds AS w ON ST_Intersects(w.the_geom,d.the_geom)
   WHERE w.ws_name = 'Sussex'