Introduction to PostGIS

Geometry, Geography, and Raster

Regina Obe and Leo Hsu

http://www.postgis.us

http://www.paragoncorporation.com

Buy our books

PostGIS In Action 2nd Edition better be out in print Q1 2015. Can buy draft now.

PostgresQL: Up and Running 2nd Edition Just went to print.

Essential PostGIS RefCardz

Buy other people's book

PostGIS Cookbook (covers up to 2.1)

Agenda

Prerequisites

Spatially enable a database

Primer on kinds of PostGIS spatial types

Load Spatial Data

Geometry Queries

Geography Queries

Tiger Geocoder

Raster

Prerequisites

Slides posted at http://bit.ly/1zfXMhK

Refer to our quick start guide if missing any PostGIS 2.1 install guide

You need PostgreSQL 9.3 or higher.

You need PostGIS 2.1 installed in your PostgreSQL instance.

Viewing tools: OpenJump or QGIS

pgAdmin

Spatial data types

• geometry Vector. Cartesian coordinate. Lots of functions
• geography Vector. Spherical coordinate. Fewer functions.
• raster Matrix of pixels. Multiple bands possible. Rapidly growing list of functions.

Geometry data type

Create table with point column

                        
DROP TABLE IF EXISTS worlds;

CREATE TABLE worlds (
id serial, 
name varchar(20), 
geom geometry(POINT), 
CONSTRAINT pk_worlds 
PRIMARY KEY(id)
);

CREATE INDEX idx_words_geom ON worlds USING gist (geom);
                        
                    

Add more spatial geometry columns

                        
ALTER TABLE worlds ADD COLUMN points geometry(point);
ALTER TABLE worlds ADD COLUMN linestrings geometry(linestring);
ALTER TABLE worlds ADD COLUMN polygons geometry(polygon);
                        
                    

See list of geometry columns

                        
SELECT f_table_name as tbl, f_geometry_column As geoc, srid, type 
FROM geometry_columns;
                        
                    

tbl            | geoc        | srid | type
----------------+-------------+------+------------
raster_columns | extent      |    0 | GEOMETRY
worlds         | geom        |    0 | POINT
worlds         | points      |    0 | POINT
worlds         | linestrings |    0 | LINESTRING
worlds         | polygons    |    0 | POLYGON
                    

Add features to our world

                        
TRUNCATE TABLE worlds;
                    
INSERT INTO worlds (name, points) 
VALUES 
    ('City Hall', ST_GeomFromText('POINT(0 0)')),
    ('Fire Station', ST_GeomFromText('POINT(2 4)')),
    ('Tower', ST_GeomFromText('POINT(10 0)'));

INSERT INTO worlds (name, linestrings) 
VALUES ('Main St', ST_GeomFromText('LINESTRING(-12 0,0 0,12 12)'));
    
INSERT INTO worlds (name, polygons) 
VALUES (
    'Centerville', 
    ST_GeomFromText('
    POLYGON((-10 -10,-10 10,10 10, 10 -10,-10 -10))'
    )
);
                        
                    

View in OpenJump or QGIS

Alter data type

                        
ALTER TABLE worlds 
ALTER COLUMN geom TYPE geometry 
USING(COALESCE(polygons,linestrings,points));
                        
                    

Geography data type

Exactly the same as geometry except use ST_GeogFromText and geography type. Spatial reference system defaults to 4326 (WGS 84 Lon Lat) if not specified. Coordinates always in decimal degrees.

                        
ALTER TABLE worlds ADD COLUMN geog geography;
UPDATE worlds 
SET geog = COALESCE(polygons,linestrings,points)::geography;
                        
                    

List geography columns

                        
SELECT f_table_name as tbl, f_geography_column As geoc, srid, type 
FROM geography_columns;
                        
                    

tbl     | geoc | srid | type
--------+------+------+----------
worlds  | geog |    0 | Geometry
                    

Geography units

Coordinates x=Lon, y=Lat

Measurements always in meters

                        
SELECT name, ST_Area(geom) As deg_aread, ST_Area(Geog) as m_area
FROM worlds
WHERE name ='Centerville';
                        
                    

name        | deg_aread | m_area
------------+-----------+------------------
Centerville |       400 | 4969694980894.16
                    

Spatial reference systems

Cataloged in spatial_ref_sys table.

Geoid

Ellipsoid

Common spatial projections

Three most popular ones

Plate Carrée

AKA as Lon-Lat in geometry.

• geometry in 4326 is Cartesian, NOT spheroidal
• geometry in 4326 measurements are pseudo degrees
• geometry in 4326 is NOT geography in 4326

Plate Carrée: 4326 (WGS 84 Lon Lat) in geometry

Geodetic (Geography) 4326 (WGS 84 Lon Lat) in geography

Web Mercator

Goes by many SRIDs : 900913 (defunct), EPSG:3857, EPSG:3785 (defunct)

Functions for working with spatial reference systems

ST_Transform

Converts from one spatial reference system to another. (Supports: geometry and raster)

ST_SetSRID

Embeds the spatial reference system in the geometry and raster.

Working with real data

We'll talk about getting data, loading data, and querying data.

Vector data sources

Office of Geographic Information (MassGIS Datalayers)

Datasets we'll use in Massachusetts State Plane Meters (SRID: 26986)

• MBTA rapid transit
• Bus routes and stops
• Cities and towns polygons

Tools packaged with PostGIS for importing and exporting vector data

• shp2pgsql — Command-line tool for loading ESRI shape files and DBase file
• pgsql2shp — Command-line tool for exporting ESRI shape file and DBase Files
• shp2pgsql-gui — GUI for importing and exporting. Not packaged with all distributions of PostGIS. You'll find this in PostGIS windows Stackbuilder distribution and Boundless OpenGeo distributions.
• raster2pgsql — Command-line tool for loading rasters

Other popular tools for loading data

• GDAL command-line tools — ogr2ogr for vector data, gdal_translate, gdalwarp for loading raster data
• osm2pgsql — Command-line tool for loading OpenStreetMap data to PostGIS (Linux, Mac, Windows).
• imposm — Tool for loading OSM to PostGIS. No binaries available for Windows. (Only tested on Linux and Mac)
• PostgreSQL FDWs — For transferring between data sources.
• oracle_fdw — Has SDO_GEOMETRY to PostGIS geometry translation. Binaries available for windows and source for easy compile on Linux/Unix. Requires Oracle OCI.dll.

Basic query to try to find matching SRID

                        
SELECT srid,  srtext, proj4text
FROM spatial_ref_sys 
WHERE 
    srtext ILIKE '%NAD83%' AND 
    srtext ILIKE '%Massachusetts%' AND 
    srtext ILIKE '%UNIT["met%",1%';
                        
                    

 srid  |                            srtext                            |      proj4text
-------+--------------------------------------------------------------+----------------------
  2805 | PROJCS["NAD83(HARN) / Massachusetts Mainland",GEOGCS["NAD83( | +proj=lcc +lat_1=42.
  2806 | PROJCS["NAD83(HARN) / Massachusetts Island",GEOGCS["NAD83(HA | +proj=lcc +lat_1=41.
 26986 | PROJCS["NAD83 / Massachusetts Mainland",GEOGCS["NAD83",DATUM | +proj=lcc +lat_1=42.
 26987 | PROJCS["NAD83 / Massachusetts Island",GEOGCS["NAD83",DATUM[" | +proj=lcc +lat_1=41.
  3583 | PROJCS["NAD83(NSRS2007) / Massachusetts Island",GEOGCS["NAD8 | +proj=lcc +lat_1=41.
  3585 | PROJCS["NAD83(NSRS2007) / Massachusetts Mainland",GEOGCS["NA | +proj=lcc +lat_1=42.
                    

Setting path variable

This is generally more of an issue on windows especially since you might have multiple versions of PostgreSQL.

SET PATH=%PATH%;"C:\Program Files\PostgreSQL\9.3\bin"

Setting PostgreSQL specific variables

Useful for loading data using shp2pgsql or raster2pgsql

Windows users replace export with SET

export PATH=...

                        
export PGPORT=5432
export PGHOST=localhost
export PGUSER=postgres
export PGPASSWORD=whatever
export PGDATABASE=avidgeo
                        
                    

Basic Geometry Load

-g column_name, -I index, -s spatial_ref sys, -S single, -D copy mode

                        
shp2pgsql -s 26986 -D -g geom -I -W "latin1" "MBTA_ARC" mbta_arc | psql
shp2pgsql -s 26986 -D -I -W "latin1" "MBTA_NODE" mbta_node | psql
shp2pgsql -s 26986 -D -I "MBTABUSSTOPS_PT" mbta_bus_stops | psql
shp2pgsql -s 26986 -D -I -S "MBTABUSROUTES_ARC" mbta_bus_routes | psql
                        
                    

Basic Geography Load with transform

-g column_name, -G force geography, -I index, -s spatial_ref sys. You can't use -D if you will transform -s from:to. So it's a bit slower.

shp2pgsql -s 26986:4326 -G -I -W "latin1" "TOWNS_POLY" towns | psql

Adding corresponding geometry /geography columns

                        
ALTER TABLE towns ADD COLUMN geom geometry (MULTIPOLYGON, 26986);
UPDATE towns SET geom = ST_Transform(geog::geometry, 26986);
CREATE INDEX idx_towns_geom ON towns USING gist(geom);                   
                        
                    

Inspecting spatial catalogs

SELECT * FROM geometry_columns;

SELECT * FROM geography_columns;

Geometry queries

Units are always in the units of our spatial reference system. For our Massachusetts data that is state-plane meters, our units are in meters.

Stations within 1000 meters

Use ST_DWithin

                        
WITH ref As (
    SELECT 
        ST_Transform(
            ST_GeomFromText(
                'POINT(-71.0824 42.3659)', 4326
            ),
            26986
        ) As loc
)
SELECT n.station, n.route, ST_Distance(n.geom, ref.loc) As dist_m 
FROM 
    mbta_node As n 
    INNER JOIN 
    ref 
    ON ST_DWithin(n.geom, ref.loc, 1000)
ORDER BY dist_m;
                        
                    

   station   |                  route                |      dist_m
-------------+---------------------------------------+------------------
 Kendall/MIT | A - Ashmont B - Braintree C - Alewife | 497.356410111099
 Lechmere    | E - Lechmere                          | 705.062362967058
                    

Three closest stops

Use KNN distance operators

                        
WITH ref As (
    SELECT 
        ST_Transform(
            ST_GeomFromText(
                'POINT(-71.0824 42.3659)', 4326
            ),
            26986
        ) As loc
)
SELECT 
    n.station, 
    n.route, 
    ST_Distance(
        n.geom, 
        (SELECT loc FROM ref)
    )::numeric(10,2) As dist_m,
    (n.geom <-> (SELECT loc FROM ref))::numeric(10,2) As dist_knn 
FROM mbta_node As n 
ORDER BY n.geom <-> (SELECT loc FROM ref) 
LIMIT 3;
                        
                    

station     |                  route                  | dist_m  | dist_knn
------------+-----------------------------------------+---------+---------
Kendall/MIT | A - Ashmont  B - Braintree  C - Alewife |  497.36 |   497.39
Lechmere    | E - Lechmere                            |  705.06 |   705.08
Charles/MGH | A - Ashmont  B - Braintree  C - Alewife | 1070.40 |  1070.39
                    

Longest bus route

Use ST_Length

                        SELECT mbta_route, SUM(ST_Length(geom)) As len_m 
FROM mbta_bus_routes 
GROUP BY mbta_route 
ORDER BY len_m DESC LIMIT 5;
                        
                    

 mbta_route |      len_m
------------+------------------
 34E        | 440190.649593145
 134        | 213819.732897472
 354        | 196239.007361261
 240        | 184848.966892802
 435        | 177963.889656812

TIGER geocoder

PostGIS 2.2 dev (trunk) been upgraded to use newly released Tiger 2014 dataset. You can use the postgis_tiger_geocoder extension files. Copy it into your PostgreSQL from Windows Builds, even if on Linux.

Install if you haven't already

                        
CREATE EXTENSION fuzzystrmatch SCHEMA contrib; 
CREATE EXTENSION address_standardizer SCHEMA contrib; -- optional
CREATE EXTENSION postgis_tiger_geocoder;
                        
                    

Standardizing addresses

Usually the first step before geocoding. TIGER without street data can normalize in most cases, but not all.

Verify your install of regular norm_addy

                        
SELECT na.address, na.streetname, na.streettypeabbrev, na.zip
FROM normalize_address('275 Third St, Cambridge, MA, 02142') AS na;
                        
                    

 address | streetname | streettypeabbrev |  zip
---------+------------+------------------+-------
     275 | Third      | St               | 02142
                    

Geocoding

                        
SELECT pprint_addy(addy), ST_AsText(ST_SnapToGrid(geomout,0.0001)), rating 
FROM geocode('275 3rd St, Cambridge, MA, 02142',1);
                        
                    

           pprint_addy           |        st_astext        | rating
---------------------------------+-------------------------+--------
275 3rd St, Cambridge, MA 02142 | POINT(-71.0824 42.3659) |      0
                    

Puts us at the door. Google geocodes to 42.365946, -71.082573 which puts us inside the building.

Geocoding intersections

Intersection name on 3rd St

                        
SELECT 
    pprint_addy(addy), 
    ST_AsText(ST_SnapToGrid(geomout,0.0001)) As loc, 
    rating 
FROM geocode_intersection('3rd St','Binney St','MA','Cambridge');
                        
                    

          pprint_addy           |           loc           | rating
--------------------------------+-------------------------+--------
282 3rd St, Cambridge, MA 02141 | POINT(-71.0824 42.3657) |      0
281 3rd St, Cambridge, MA 02142 | POINT(-71.0824 42.3657) |      0
280 3rd St, Cambridge, MA 02141 | POINT(-71.0824 42.3657) |      0
279 3rd St, Cambridge, MA 02142 | POINT(-71.0824 42.3657) |      0
                    

Geocoding Intersections

Intersection name on Binney St

                        
SELECT 
    pprint_addy(addy), 
    ST_AsText(ST_SnapToGrid(geomout,0.0001)) As loc, 
    rating 
FROM geocode_intersection('Binney St','3rd St','MA','Cambridge');
                        
                    

           pprint_addy             |           loc           | rating
-----------------------------------+-------------------------+--------
148 Binney St, Cambridge, MA 02142 | POINT(-71.0824 42.3657) |      0
149 Binney St, Cambridge, MA 02143 | POINT(-71.0824 42.3657) |      0
146 Binney St, Cambridge, MA 02142 | POINT(-71.0824 42.3657) |      0
147 Binney St, Cambridge, MA 02142 | POINT(-71.0824 42.3657) |      0
                    

Geocoding intersections, include zip and number matches

You get better performance if you provide zip and max returns

                        
SELECT 
    pprint_addy(addy), 
    ST_AsText(ST_SnapToGrid(geomout,0.00001)) As loc, 
    rating 
FROM geocode_intersection('Binney St','3rd St','MA','Cambridge','02142',1);
                        
                    

           pprint_addy             |            loc            | rating
-----------------------------------+---------------------------+--------
148 Binney St, Cambridge, MA 02142 | POINT(-71.08239 42.36575) |      0
                    

Reverse Geocoding

addy is an array of possible normalized addresses of a point location. You can expand aliases with unnest.

                        
SELECT pprint_addy(unnest(rc.addy)) 
FROM reverse_geocode(ST_Point(-71.0824, 42.3659)) As rc;
                        
                    

pprint_addy
----------------------------------
277 3rd St, Cambridge, MA 02142
151 Binney St, Cambridge, MA 02143
                    

Working with Rasters

Download some raster data

• Wind Speed data (ESRI GRID format -- driver built in to gdal by default)
• Color Orthos (JP2 ECW need your raster2pgsql built with MrSID or ECW (Jasper seems unable to handle these).

Unzip and extract into folder

Load with raster2pgsql

raster2pgsql -s 26986 -F -I -C -M wind/ne_spd_30m -t 200x200 wind_spd_30m | psql -d avidgeo

Aerials are in UTM Zone 10 NAD 83 meters SRID:26919

Load and generate -I add spatial index, -l overviews (really speeds up display in QGIS), -Y use copy, -e no transaction, -F include filename, -C add constraints, -t tile size, -M analyse (update planner stats)

PostGIS 2.1 raster2pgsql can guess at the SRID if enough info in meta data. So you can skip the SRID.

raster2pgsql -Y -e -F -I -C -M aerials/*.jpg -t 500x500 -l 2,4 aerials | psql -d avidgeo

Get wind speed at 30 meters above ground at point of interest

                        
SELECT ST_Value(rast, 1, geom) As wdval 
FROM 
    wind_spd_30m As w 
    INNER JOIN 
    ST_Transform(
            ST_GeomFromText('POINT(-71.082573 42.365946)',4326),26986
    ) As geom
    ON ST_Intersects(w.rast, geom);
                        
                    

wdval
----------------
5.16392469406128
                    

Get histogram of wind speed

500 meter box around point of interest, 3 bins

                        
SELECT (h).*						 
FROM (
    SELECT ST_Histogram(ST_Union(ST_Clip(rast,geom)),1,3) As h
    FROM 
        wind_spd_30m As w 
        INNER JOIN 
        ST_Expand(
            ST_Transform(
                ST_GeomFromText('POINT(-71.082573 42.365946)',4326),
                26986
            ),
            300
        ) As geom
        ON ST_Intersects(w.rast,geom)
) AS c;
                        
                    

      min        |       max        | count |      percent
-----------------+------------------+-------+------------------
4.92105531692505 | 5.04267676671346 |     3 | 0.333333333333333
5.04267676671346 | 5.16429821650187 |     2 | 0.222222222222222
5.16429821650187 | 5.28591966629028 |     4 | 0.444444444444444
                    

Output as aerial PNG area of interest

This took 2781 ms

                        
SELECT 
    ST_AsPNG(
        ST_Resize(ST_Union(ST_Clip(w.rast,geom)),0.5,0.5)
    ) As png 
FROM 
    aerials AS w 
    INNER JOIN 
    ST_Buffer(
        ST_Transform(
            ST_GeomFromText('POINT(-71.0824 42.3659)',4326),
            26919
        ),
        150
    ) As geom
    ON ST_Intersects(w.rast,geom);
                        
                    

Use overview table

Using overview is much faster (o2 = ~921 ms, o4 = ~270ms). o2 achieves same result.

                        
SELECT ST_AsPNG(ST_Union(ST_Clip(w.rast, geom))) As png 
FROM 
    o_2_aerials AS w 
    INNER JOIN 
    ST_Buffer(
        ST_Transform(
            ST_GeomFromText('POINT(-71.0824 42.3659)',4326),
            26919
        ),
        150
    ) As geom
    ON ST_Intersects(w.rast,geom);
                        
                    

View in QGIS

                        
CREATE OR REPLACE VIEW aerial_test AS 
    SELECT 1 As rid, ST_Union(ST_Clip(w.rast, geom)) As rast
    FROM 
        o_2_aerials AS w INNER JOIN 
        ST_Buffer(
            ST_Transform(
                ST_GeomFromText('POINT(-71.0824 42.3659)', 4326),
                26919
            ), 
            150
        ) As geom
        ON ST_Intersects(w.rast,geom);
                        
                    

Then in DBManager drag the view onto map.

Links of Interest

• PostGIS
• Planet PostGIS

THE END

Buy our books here http://www.postgis.us