[postgis-users] Example queries - a compendium

[Paul Ramsey]

Shaun Kolomeitz wrote:
> Dear PostGIS'ers,
> 
> Has anyone got a reference for, or created for themselves a compendium
> of useful examples of SQL GIS type queries ?
> I'm looking for a good smattering of examples covering a range of
> feature types (points, lines, poly's) and the sorts of spatial
> operations one is likely to require (from basic/simple to intermediate
> and advanced involving complex joins). It would be nice if this could
> cover simple reading of features, creation, through to advanced updating
> techniques for features.
> I've seen a few examples included with the PostGIS Manual as well as
> some HTML type help that comes with PostGresQL's install of PostGIS.
> These help, but don't seem to go far enough.
> 
> Any pointers most appreciated.
> As a quid-pro-quo I'd re-distribute back to the group a summary.
> 
> Kind regards,
> Shaun
> 
> ___________________________
> Disclaimer
> 
> WARNING: This e-mail (including any attachments) has originated from a Queensland Government department and may contain information that is confidential, private, or covered by legal professional privilege, and may be protected by copyright.  
> 
> You may use this e-mail only if you are the person(s) it was intended to be sent to and if you use it in an authorised way.  No one is allowed to use, review, alter, transmit, disclose, distribute, print or copy this e-mail without appropriate authority.  If you have received this e-mail in error, please inform the sender immediately by phone or e-mail and delete this e-mail, including any copies, from your computer system network and destroy any hardcopies.
> 
> Unless otherwise stated, this e-mail represents the views of the sender and not the views of the Environmental Protection Agency.
> 
> Although this e-mail has been checked for the presence of computer viruses, the Environmental Protection Agency provides no warranty that all viruses have been detected and cleaned. Any use of this e-mail could harm your computer system.  It is your responsibility to ensure that this e-mail does not contain and is not affected by computer viruses, defects or interference by third parties or replication problems (including incompatibility with your computer system).
> 
> E-mails sent to and from the Environmental Protection Agency will be electronically stored, managed and may be audited, in accordance with the law and Queensland Government Information Standards (IS31, IS38, IS40, IS41 and IS42) to the extent they are consistent with the law.
> 
> ___________________________
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users


-- 

   Paul Ramsey
   Refractions Research
   http://www.refractions.net
   pramsey at refractions.net
   Phone: 250-383-3022
   Cell: 250-885-0632
-------------- next part --------------
SQL Commands for the Demonstration

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);