Essential PostGIS

by: Leo Hsu and Regina Obe

PostGIS is a free, open source spatial extender for the PostgreSQL open source database. Spatial extenders leverage existing database capabilities to handle geographical and image data. Common applications include proximity analysis for addressing questions of how far something is, geocoding for finding latitude and longitude of addresses, and web mapping for supplying data to various map servers.

PostGIS has features equal to, and often exceeding, popular commercial offerings such as Oracle Spatial and SQL Server. Like Oracle, it offers 2D and 3D support for geometry. PostGIS geography type mimics what you’ll find in SQL Server.

This Refcard provides a quick reference for common tasks you perform when using PostGIS. We based our examples on PostgreSQL 9.1+ and PostGIS 2.1. For a thorough treatment of PostGIS, we encourage you to buy PostGIS In Action, 2nd Edition and to chomp the PostGIS manuals at http://postgis.net/documentation.

This Refcard assumes that you are reasonably familiar with PostgreSQL and have more or less heard of basic concepts required to work with geographical and image data. No experience with PostGIS is assumed. For more on PostgreSQL, check out our Essential PostgreSQL Refcard (http:// refcardz.dzone.com/refcardz/essential-postgresql).

To find installer packages for most operating systems, visit http://postgis. net/install. To compile from source, visit http://postgis.net/source. Make sure you already have a working PostgreSQL server.

To create a new database in PostgreSQL, execute the following command in psql or use pgAdmin:

CREATE DATABASE mygisdb;

​

You then need to enable PostGIS for your new database. To do so, you first connect to your database with psql and then run one of the following:

PostGIS introduces four families of spatial data types:

shp2pgsql is a command-line tool for importing ESRI shape files. Use it with psql. Shp2pgsql generates a SQL script that you then execute to perform the import itself.

Common Options

Example

Generate script using (loads into geography format):

shp2pgsql -s 4326 -G -I shapefile.shp a_table > import.sql

​

Execute script using psql:

psql -d my_db -f import.sql

​

Load wgs84 long lat (4326) and transform to web Mercator (3857):

shp2pgsql -s 4326:3857 shapefile.shp a_table | psql -d my_db

​

raster2pgsql is a command-line tool for importing raster data. Supports many formats of raster data depending on GDAL library.

Common Options

Example

Generate sql load script for new project_pictures table (loads in all pictures in a folder):

raster2pgsql -e -F *.jpg project_pictures > pics.sql

​

Execute script using psql:

psql -d my_db -f pics.sql

​

Load aerial coverage tiling 512x512, index, constraints in one step:

raster2pgsql -I -e -C *.jpg aerials -t 512x512 | psql -d my_db

​

Geometry, geography, and raster data types often share the same function, but not always. You can usually intuit which functions apply or don’t apply to a specific data type. For instance, you’d never try to count the number of bands against a geometry type.

Not all functions require parameters of the same data type. For instance, to see if a point lies within a raster, you’d pass a geometry type and a raster type to the ST_Intersects function.

PostGIS offers close to a thousand functions, each with multiple overloads, but mastering only about a dozen or so for each data type will suffice.

We cover the most commonly used functions below along with the most typical arguments. To explore the ever growing set of functions, visit http:// postgis.net.

Functions and Operators

Geometry and Geography Functions

Raster Functions

Topology Functions

These are examples that showcase some PostGIS functions.

Find All Places Within 100 Meters of a Road

SELECT p.place_name, ST_Distance(p.geog, r.geog) As distFROM places AS p INNER JOIN roads As r ON ST_DWithin(p.geog, r.geog, 100)WHERE r.road_name = ‘Some Road’;

​

Find Ten Closest Places to a Given Point

Using indexed nearest neighbor bounding box centroid operator. Only approximate because the check is approximate because we’re treating the earth as planar instead of spherical.

SELECT postal_code 

FROM postal 

ORDER BY ST_SetSRID(ST_Point(-71.06, 42.36), 4326) <-> geom LIMIT 10;

​

For Each Place Find the Closest Hospital

Using indexed distance operator - approximate if data is long lat or if place and hospital are not both points:

SELECT p.place_name  , (SELECT h.place_name        FROM places As h        WHERE h.type = 'hospital'         ORDER BY  h.geom <-> p.geom LIMIT 1) AS closest_hospitalFROM places As p;

​

For Each Postal Find the Closest Road

This example assumes our geometry is stored in WGS 84 long lat. Use indexed box operator <#> to find 100 closest roads by bounding box distance and then narrow down to true geometry proximity by doing a spherioidal distance check.

SELECT p.postal_code  , ( SELECT road_name        FROM (SELECT r.road_name , ST_Distance(r.geom::geography      , p.geom::geography) As dist_m 

       FROM roads As r        ORDER BY  r.geom <#> p.geom LIMIT 100) As span_1ORDER BY dist_m LIMIT 1 ) AS closest_roadFROM postal As p WHERE place_name = 'Boston';

​

Number and Limit Records by Distance Order

For each road segment, find the 2 closest postals within 1000 meters of the road. If no postal within 1000 meters, then NULL is returned for postal_code.

SELECT rnum,road_name, postal_code FROM (SELECT r.gid , ROW_NUMBER() OVER(PARTITION BY r.gid         ORDER BY ST_Distance(p.geog,r.geog) ) As rnum , p.postal_code, r.road_nameFROM  roads As r LEFT JOIN  postal As p    ON ST_Dwithin(p.geog, r.geog,1000) ) As xWHERE rnum < 3ORDER BY x.road_name, x.gid;

​

Rectangle Box

Returns geometries within the WGS 84 planar rectangle defined by ST_ MakeEnvelope(min_lon,min_lat,max_lon,max_lat,4326):

SELECT geomFROM roads As rWHERE r.geom && ST_MakeEnvelope(-71.06,42.36,-70.10,41.40,4326);

​

Returns geographies within the WGS 84 spheroid rectangle defined by ST_MakeEnvelope(min_lon,min_lat,max_lon,max_lat,4326):

SELECT geogFROM roads As rWHERE r.geog  && ST_MakeEnvelope(-71.06,42.36,-70.10,41.40,4326)::geography;

​

Insert Geometry Statement

INSERT INTO roads(road_name,geom)VALUES (‘Main St’, ST_GeomFromText(‘LINESTRING(-71.06 42.36

,-71.50 42.50)’, 4326));

​

Insert Geography Statement

INSERT INTO roads(road_name, geog)VALUES (‘Main St’, ST_GeogFromText(‘LINESTRING(-71.06 42.36, 

-71.50 42.50)’));

​

Update Geometry Point

Updates a geometry column with data from lon,lat column

 UPDATE postal SET geom = ST_SetSRID(ST_Point(lon,lat),4326);

​

Update Geography Point

Updates a geometry column with data from lon,lat column:

 UPDATE postal SET geog = ST_Point(lon,lat)::geography;

​

Add Missing Roads

Insert statement, selects only those records from rnew table with road geometry not already present in roads:

INSERT INTO roads(road_name,geom)SELECT rnew.road_name, rnew.geomFROM rnew   LEFT JOIN roads ON (ST_OrderingEquals(rnew.geom, roads.geom) )WHERE roads.gid IS NULL;

​

Spatial Update Count

For each county populate a field with count of restaurants in county:

UPDATE counties 

 SET num_restaurants = (SELECT COUNT(r.id) 

     FROM restaurants As r

       WHERE ST_Intersects(counties.geom,r.geom) );

​

Linestring path from GPS points

This example uses an ordered aggregate clause to ensure the linestring formed is in order of time travel. For each user, for each day will create linestring representing travel for that day:

SELECT user_name,travel_ts::date As travel_date , ST_MakeLine(geom ORDER BY travel_ts) AS geom FROM gps_pointsGROUP BY user_name, travel_ts::date;

​

Return Elevation at Each Point in a Table of Points

SELECT p.place_name, ST_Value(r.rast,1,p.geom) as elev

FROM places As p INNER JOIN dems As r

  ON ST_Intersects(p.geom, r.rast);

​

Bounding Box Raster Filter

SELECT rast

FROM dem As rWHERE r.rast && ST_MakeEnvelope(-71.06,42.36,-70.10,41.40,4326);

​

Return a Raster Clipped to Within 100 Unit Box of a Parcel

SELECT ST_Union(ST_Clip(a.rast,ST_Expand(p.geom,100)))FROM aerials As a    INNER JOIN      parcels As p ON ST_DWithin(a.rast::geometry, p.geom, 100)WHERE p.pid = ‘57-169E’;

​

Resize and Output to PNG

Creates a rast of 25% of original size using Lanczos resampling algorithm and then outputs to PNG format

SELECT ST_AsPNG(     ST_Resize(rast, 0.25,0.25, ‘Lanczos’) ) As pngFROM project_pictures

​

Convert a Geometry to a Raster Using Reference Raster

Rasters formed will be 1 band 8BUI value 1 (defaults when not specified) and pixel size will be the same and clipped to the extent of the dem rasters and return first.

SELECT ST_AsRaster(r.geom, d.rast) As rastFROM roads r INNER JOIN dems d ON ST_Intersects(r.geom,d.rast)WHERE r.road_name = ‘Some Road’ LIMIT 2;

​

Convert to Geometry Elevation Range for Particular Bounding Box

Raster tiles that intersect bounding box and of a particular elevation range, convert to polygon/multipolygon

SELECT ST_Union( (gval).geom)FROM (SELECT ST_DumpAsPolygons(rast,1, true) As gvalFROM usgs_srtmWHERE ST_Intersects(rast  , ST_MakeEnvelope(85.99958, 49.0004, 85.99959, 49.0005, 4326) )   ) AS elevWHERE (gval).val BETWEEN 2403 and 2405

​

Convert to PostGIS Raster, Resize, Convert to PNG

Create 100px width png thumbnails of a JPEG or other supported raster formats:

WITH x AS(SELECT ST_FromGDALRaster(pic_image) As rast   FROM project_pictures    WHERE project = ‘Gulliver Redevelop’)SELECT ST_AsPNG(ST_Resize(rast, 100, (ST_Height(rast)/ST_Width(rast)*1.0*100)::integer) ) As thumbFROM x;

​

The spatial catalog tables provide an inventory of spatial columns you have in your database as well as additional properties about them often used by third-party tools and viewers. They are similar in concept to the standard information_schema.columns table of PostgreSQL proper.

These are just a few of the most popular PostGIS extensions.

pgrouting

http://pgrouting.org

FOSS PostgreSQL extension that extends PostGIS to provide geospatial routing functionality. pgRouting 2.0+ can now be installed with

CREATE EXTENSION pgrouting;

​

PL/R

http://www.joeconway.com/plr

PostgreSQL language extension that allows writing stored functions in R.

CREATE EXTENSION plr;

​

hstore

Key/value column type that allows storing and quering multiple ad-hoc key- values in a single field with very fast indexing support using PostgreSQL gist index. Comes packaged with all binary distributions of PostgreSQL.

CREATE EXTENSION hstore;

​

All spatial indexes use the PostgreSQL gist index type for indexing, but have a different operator class that controls what operations can be used on them and how those operations behave.

PostGIS comes bundled with several tools for loading and exporting spatial data. There are many other FOSS tools that are commonly used with PostGIS as well. Listed below is a common subset.

Usage:

pgsql2shp [<options>] <database> [<schema>.]<table>pgsql2shp [<options>] <database> <query>

​

Some common options:

Get help:

pgsql2shp

​

Dump out a table:

pgsql2shp -f ca_roads -g geog gisdb ca.roads 

​

Dump out a query and reproject:

pgsql2shp -f roads_merc gisdb “SELECT road_name, ST_Transform(geog::geometry, 3857) AS geom FROM ca.roads”

​

Backup and Restore

Below are common backup and restore statements. Note: backing up a spatial database is not different from backing up any other database.

Create a compressed backup:

pg_dump -h localhost -p 5432 -U postgres -F c -b -v -f “/

somepath/gisdb.backup” gisdb

​

Restore compressed backup without upgrade:

psql -h localhost -d postgres -U postgres -c “CREATE DATABASE gisdb” 

​

Upgrading

You’ll need to do a hard upgrade if upgrading from major to new major: e.g. 1.5 -> 2.0 requires backup and restore, but 2.0->2.1 just requires standard upgrade with ALTER EXTENSION. You need to install the legacy script if you have applications that only support 1.5

psql —h localhost -d gisdb -U someuser -c “CREATE EXTENSION postgis;”psql -h localhost -d gisdb -U someuser -f legacy.sql

​

Optimization

After a large load or delete, to ensure good planner stats:

vacuum analyze verbose sometable;

​

Uninstalling PostGIS

Should you ever need to uninstall, run:

DROP EXTENSION postgis;

​

To both uninstall and drop all spatial columns, stored data, and dependent extensions:

DROP EXTENSION postgis CASCADE;

​