Skip to content
This repository has been archived by the owner on Jan 16, 2019. It is now read-only.

PostGIS

Jump to bottom
Falk Schuetzenmeister edited this page May 8, 2015 · 7 revisions

Only if you did not attend yesterday's workshop

CREATE ROLE dbuser LOGIN SUPERUSER CREATEDB;
ALTER ROLE dbuser PASSWORD 'unsicher';
CREATE DATABASE addresses;
ALTER DATABASE addresses OWNER to dbuser;

Or do the equivalent steps in pgadmin.

Then connect to Database addresses using user dbuser.

Schemas

Schemas address three problems:

  • different user permissions for core data and user views and experiments
  • SQL queries do not work across Databases but across schemas in the same database allows for data organization
  • sort of a folder system to organize data

Create Schema

DROP SCHEMA IF EXISTS geostuff;
CREATE SCHEMA geostuff AUTHORIZATION dbuser;

Getting Geodata

Preparing the database for geospatial data - Installing the PostGIS extension

Postgis needs to be installed on the system (above we choose packages that were already bundled with Postgis).

CREATE EXTENSION postgis;

Import some geo data (ESRI shapefile)

FROM COMMAND LINE:

shp2pgsql -s 4269 cb_2013_us_state_500k.shp geostuff.states | psql -U dbuser -d addresses

Or use

shp2pgsql-gui (Ubuntu, Mac?)

pgShapeloader (Windows)

Use same parameters as above.

See also http://suite.opengeo.org/4.1/dataadmin/pgGettingStarted/pgshapeloader.html

  • A good way to look at and evaluate the result is [QGIS] (http://www.qgis.org/en/site/). It allows direct connections to your Postgresql installation.

Creating geometries see http://postgis.net/docs/using_postgis_dbmanagement.html

BEGIN;
   -- CREATE SOME POINT LOCATIONS
   DROP TABLE IF EXISTS geostuff.locations;
   DROP SEQUENCE IF EXISTS geostuff.locations_id;
   CREATE SEQUENCE geostuff.locations_id START 1; 
   CREATE TABLE geostuff.locations (
    id integer PRIMARY KEY DEFAULT nextval('geostuff.locations_id'),
    name varchar(30),
    geom geometry(POINT, 4326));
   ALTER TABLE geostuff.locations OWNER TO dbuser;
   INSERT INTO geostuff.locations(name, geom) VALUES ('Sacramento',  ST_PointFromText('POINT(-121.468889 38.555555)', 4326));
   -- it is x y which is lon lat (and not the other way around)
   INSERT INTO geostuff.locations(name, geom) VALUES ('Berlin',  ST_PointFromText('POINT(-81.794444 40.561111)', 4326)); 
   INSERT INTO geostuff.locations(name, geom) VALUES ('Dresden',  ST_PointFromText('POINT(-82.0131 40.1214)', 4326)); 
   INSERT INTO geostuff.locations(name, geom) VALUES ('Berlin',  ST_PointFromText('POINT(-72.7725 41.613889)', 4326)); 
COMMIT;
SELECT name FROM geostuff.states WHERE ST_INTERSECTS(geom, ST_TRANSFORM((SELECT geom FROM geostuff.locations WHERE name like 'Sacramento'), 4269));

Create geo-spatial Indexes on geometry field to support geographic operations.

CREATE INDEX geom_gix ON geostuff.states USING GIST(geom);

Spatial Joins

SELECT a.name, b.name FROM geostuff.locations AS a 
    JOIN geostuff.states AS b 
	ON ST_Intersects(b.geom, ST_Transform(a.geom, 4269));

Now we are adding a city outside the US (which does not have a state associated with it). Plus, demonstrating a few things on the side (GeoJSON import and SRID assignment which is different from transform).

DELETE FROM geostuff.locations WHERE name like 'Canberr%';
INSERT INTO geostuff.locations (name, geom) 
    VALUES ('Canberra', 
            ST_SetSRID(ST_GeomFromGeoJSON(
		'{"type": "Point", "coordinates": [149.124417, -35.3075]}'),
		4326));
SELECT name, ST_AsGeoJSON(geom) as geojson FROM geostuff.locations;

Now we run the example above again. Compare the effect of the following queries:

-- this one is identical with just JOIN
SELECT a.name, b.name FROM geostuff.locations AS a 
    INNER JOIN geostuff.states AS b 
	ON ST_Intersects(b.geom, ST_Transform(a.geom, 4269));
-- this one is identical with LEFT OUTER JOIN
SELECT a.name, b.name FROM geostuff.locations AS a 
    LEFT JOIN geostuff.states AS b 
	ON ST_Intersects(b.geom, ST_Transform(a.geom, 4269));
-- this one is identical with RIGHT OUTER JOIN
SELECT a.name, b.name FROM geostuff.locations AS a 
    RIGHT JOIN geostuff.states AS b 
	ON ST_Intersects(b.geom, ST_Transform(a.geom, 4269));

Taking it up a notch: Triple Join

SELECT a.first_name, a.last_name, a.city, c.state FROM people AS a LEFT JOIN (SELECT aa.name AS city, bb.name AS state 
	FROM geostuff.locations AS aa
	JOIN geostuff.states AS bb
	ON ST_INTERSECTS(aa.geom, ST_TRANSFORM(bb.geom, 4326))) as c
	    ON a.city = c.city;
  • Be aware of the problem that occurs from ambiguous joins (by name) here: records with all combinations get created.

Some geospatial analysis

  • The number of geospatial analysis and type conversion functions is extreme. PostGIS has also raster support. Well Create extension created about 1000 functions. You have to consult the documentation. But we start with some simple tasks.

1. Distance between Sacramento and Berlin, OH:

  • Reproject in a equidistant coordinate system, well the one I picked actually isn't.
SELECT ST_Distance(
   ST_Transform((SELECT geom FROM geostuff.locations WHERE name like 'Sacramento'), 2163),
   ST_Transform((SELECT geom FROM geostuff.locations WHERE name like 'Berlin' LIMIT 1), 2163))/1609 
   AS distance_in_miles;

2. Uniting California and Nevada

BEGIN;
    DROP TABLE IF EXISTS geostuff.regions;
    DROP SEQUENCE IF EXISTS geostuff.regions_id;
    CREATE SEQUENCE geostuff.regions_id START 1;
    CREATE TABLE geostuff.regions (
	id integer PRIMARY KEY DEFAULT nextval('geostuff.regions_id'),
	name varchar(30),
	geom geometry('Multipolygon', 4326)
    );
    ALTER TABLE geostuff.regions OWNER TO dbuser;
    INSERT INTO geostuff.regions (name, geom) VALUES (
	'CALIVADA',
	(SELECT ST_UNION(ST_TRANSFORM(geom, 4326)) FROM geostuff.states WHERE name IN ('Nevada', 'California'))
    );
COMMIT;

3. States bigger than Ohio

SELECT * FROM geostuff.states 
    WHERE ST_AREA(geom) > ST_AREA((SELECT geom FROM geostuff.states WHERE name='Ohio'))

4. States neighboring Nevada

SELECT * FROM geostuff.states 
    WHERE ST_Intersects(geom, (SELECT geom FROM geostuff.states WHERE name = 'Nevada')) 
    AND NOT name='Nevada'

5. Ouput transformations

Bounding boxes

SELECT gid, ST_Envelope(geom) FROM geostuff.states

Centroids

SELECT gid, ST_Centroid(geom) FROM geostuff.states

Buffers Simplified

Clone this wiki locally