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

Create user

CREATE ROLE avidgeo LOGIN PASSWORD 'g3o' INHERIT CREATEDB;

Create a new database

                            
CREATE DATABASE avidgeo WITH ENCODING='UTF8' OWNER=avidgeo;
ALTER DATABASE avidgeo SET search_path=public,postgis,contrib;

ALTER DEFAULT PRIVILEGES IN SCHEMA public 
GRANT ALL ON TABLES TO avidgeo WITH GRANT OPTION;

ALTER DEFAULT PRIVILEGES IN SCHEMA public
GRANT ALL ON SEQUENCES TO avidgeo WITH GRANT OPTION;

ALTER DEFAULT PRIVILEGES IN SCHEMA public
GRANT ALL ON FUNCTIONS TO avidgeo WITH GRANT OPTION;
                            
                        

Spatially enable database

We like to install PostGIS in its own schema call postgis and other extensions in a schema called contrib. To run these commands, use pgAdmin or psql. Make sure to connect to the database first.

(only super users can install C-based extensions so need to login as a superuser like postgres)

From pgAdmin or psql

                            
CREATE SCHEMA postgis;
GRANT ALL ON SCHEMA postgis TO public;
CREATE SCHEMA contrib;
GRANT ALL ON SCHEMA contrib TO public;
CREATE EXTENSION postgis SCHEMA postgis;
CREATE EXTENSION hstore SCHEMA contrib;
                            
                        

Verify installation

Exit out of psql and reconnect, then run

SELECT postgis_full_version();

Output should look something like

postgis_full_version
---------------------------------------------------------------------------------
POSTGIS="2.1.4 r12966" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March
2012" GDAL="GDAL 1.11.0, released 2014/04/16" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER
                        

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


Picture from Wikipedia: http://en.wikipedia.org/wiki/Equirectangular_projection

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

Web Mercator

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


Picture from Wikipedia: http://en.wikipedia.org/wiki/Web_Mercator

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)

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

Working with OpenStreetMap data (OSM)

We'll use osm2pgsql for loading since it works on windows as well and easiest to compile on other OS.

Binaries available: Osm2pgsql. If on windows use the one marked Cygwin.

Getting OSM data

Open Street Map data usually loaded in web mercator

Loading data with OSM2PGSQL

hstore needs to be installed before you can use --hstore-all or --hstore

-W switch prompts for password. if you are missing the default.style file get from here

                            
osm2pgsql 
    -d avidgeo -H localhost -U postgres -P 5432 -W -S default.style 
    --hstore boston_massachusetts.osm.pbf
                            
                        

Create functional indexes and views

							
CREATE INDEX idx_planet_osm_point_fgeom  
ON planet_osm_point  
USING gist( ( ST_Transform(way,26986)::geometry(POINT,26986) ) );

CREATE INDEX idx_planet_osm_point_fgeog 
ON planet_osm_point  
USING gist( ( ST_Transform(way,4326)::geography(POINT,4326) ) );

CREATE OR REPLACE VIEW vw_planet_osm_point AS 
SELECT *,  
	ST_Transform(way,26986)::geometry(POINT,26986) As geom,
	ST_Transform(way,4326)::geography(POINT,4326) As geog
FROM planet_osm_point;
					

Geography and Hstore queries

5 closest restaurants


WITH loc as ( SELECT ST_GeogFromText('POINT(-71.0824 42.3659)') As geog )
SELECT name,  tags->'cuisine' As cuisine,
  ST_Distance(pt.geog, loc.geog) As dist 
FROM vw_planet_osm_point AS pt INNER JOIN loc ON 
  ST_Dwithin(pt.geog, loc.geog, 500)
  AND pt.tags ? 'cuisine'
  ORDER BY dist Limit 5;
          name          |   cuisine   |     dist
------------------------+-------------+---------------
 Fuji at Kendall        | japanese    |  70.455413701
 Tatte Bakery and CafΘ  | bakery      | 109.979903628
 Voltage                | coffee_shop |   183.5050868
 Kika Tapas             | regional    | 281.472198034
 Commonwealth Cambridge | american    | 288.972055815

All restaurants within a region

--takes 230 ms, replacing geog with geom reduces to 22 ms

SELECT pt.osm_id, pt.name, 
	tags->'cuisine' As cuisine, amenity, pt.geog
 FROM vw_planet_osm_point AS pt INNER JOIN towns As t 
  ON (t.town = 'CAMBRIDGE' AND ST_Covers(t.geog, pt.geog))
  WHERE (pt.tags ? 'cuisine' 
        or amenity IN('restaurant', 'fast_food', 'cafe', 'pub')) ;

Exporting data with pgsql2pshp

Respects Postgres environment variables PGHOST, PGPORT, PGUSER or takes them as switch inputs

Output a table or view

-g to specify geometry column. only needed if more than one

pgsql2shp -f pois -g geom avidgeo vw_planet_osm_point

Output a query

SQL should be quoted and all on same line

pgsql2shp -u avidgeo -P g3o -p 5432 -f cuisine avidgeo "SELECT pt.name, pt.geom, pt.tags->'cuisine'::varchar(50) AS cuisine FROM vw_planet_osm_point AS pt WHERE pt.tags ? 'cuisine'"

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
                    

TIGER geocoder: Loading data

Tools to load data

  • unzip if on Linux/Mac, 7-zip if on Windows
  • shp2pgsql
  • wget on Linux/Mac. GNUWin32 wget on Windows.

Create staging folder

Create a folder called gisdata and within that temp

These are used for staging. Edit the tiger.loader_variables.staging_fold to have path to location.

Create copy of a template record

Here I use windows. If on Linux or Mac use sh.

                            
INSERT INTO tiger.loader_platform (
    os, declare_sect, pgbin, wget, unzip_command, psql, path_sep, loader, 
    environ_set_command, county_process_command
)
SELECT 
    'mydb', declare_sect, pgbin, wget, unzip_command, psql, path_sep, 
    loader, environ_set_command, county_process_command
FROM tiger.loader_platform
WHERE os = 'windows';
                            
                        

Use pgAdmin or favorite Postgres admin tool to edit the declare_sect.

Generate nation files load script

                            
\a \t \o c:/gisdata/nation_script.bat
SELECT Loader_Generate_Nation_Script('mydb');
\o
                            
                        

Then run the generated script from psql command line.

Tip 1: To reduce load skip block groups

If you don't need tab blocks, block groups, and tracts for statistics, don't load.

                            
UPDATE tiger.loader_lookuptables 
SET load = false 
WHERE table_name IN ('bg','tabblock','tract');
                            
                        

Generate State load script

                            
\a \t \o c:/gisdata/state_script.bat
SELECT Loader_Generate_Script('{MA}'::text[], 'mydb');
\o
                            
                        

Tip 2: Don't need whole state, edit down to county

                            
SELECT cntyidfp, name, namelsad
FROM county 
WHERE ST_Intersects(
        the_geom,
        ST_GeomFromText('POINT(-71.11678 42.36587)',4269)
);
                            
                        
 cntyidfp |   name    |     namelsad
----------+-----------+------------------
 25017    | Middlesex | Middlesex County
                        

If you don't need a whole state of data, edit the generated script.

In wget calls for faces, edges, featnames, addr, change 25 to 25017 (if just want portion of Middlesex county).

Then run the generated script from command line.

Install missing indexes and analyze

SELECT Install_Missing_Indexes();

Don't forget to update stats.

analyze verbose;

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

THE END

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