PostgreSQLisms and PostGIS

• A row can be a field.
• Rows can be aggregated, thus many rows can be a field.
• Ordered Aggregates: PostgreSQL by default allows aggregates to have an ORDER BY clause.
• All aggregates are window aggregates unless the author of the function forbids it.
• Short-hand casting operator ::
• Operators that mean different things based on types
• A function in the FROM is a lateral join and can be chained (ala cascading functions)
• Ordinality for numbering results of set returning functions
• Extension mixins

A row can be a field ..

						SELECT p.name, p
 FROM stlouis.points AS p
 WHERE name = 'Smoking Barrels BBQ';
					

-[ RECORD 1 ]-------------------------------------------------------------------------------
name | Smoking Barrels BBQ
p    | (2553435113,0101000020E6100000B53E9B0B129256C0EBF5381D124A4340,"Smoking Barrels BBQ",,,,,,,,
		"{""amenity"": ""restaurant"", ""cuisine"": ""barbecue"", ""website"
		": ""https://www.smokingbarrels.net/"", ""addr:city"": ""Saint Louis"", 
		""addr:street"": ""South Kingshighway Boulevard"", ""addr:postcode"" ...

					

In PostGIS, a row can be a geojson feature

						-- requires PostGIS 3.1+
SELECT p.name, ST_AsGeoJSON(p)::json AS p_geojson
 FROM stlouis.points AS p
 WHERE name = 'Smoking Barrels BBQ';
					

-[ RECORD 1 ]---------------------------------------------------------------------------------------------
name      | Smoking Barrels BBQ
p_geojson | {"type": "Feature", "geometry": {"type":"Point","coordinates":[-90.2823514,38.5786778]}, 
	  | "properties": {"osm_id": 2553435113,
	  | "name": "Smoking Barrels BBQ", "barrier": null, "highway": null, "ref": null,
	  | "address": null, "is_in": null, "place": null, 
	  | "man_made": null, 
	  | "other_tags": {"amenity": "restaurant", "cuisine": "barbecue",
	  | "website": "https://www.smokingbarrels.net/", "addr:city": "Saint Louis",
	  | "addr:street": "South Kingshighway Boulevard", "addr:postcode": "63109"...}}}
					

You can aggregate rows

						SELECT array_agg(p) AS p_array
 FROM stlouis.points AS p
 WHERE name ILIKE '%BBQ%';
					

						SELECT json_agg(p) AS p_json
 FROM st_louis.points AS p
 WHERE name ILIKE '%BBQ%';
					

You can aggregate a subset of columns

						SELECT json_agg(p) AS p_json
 FROM (SELECT osm_id, name FROM st_louis.points) AS p
 WHERE name ILIKE '%BBQ%';
					

[{"osm_id":2553435113,"name":"Smoking Barrels BBQ"}, 
 {"osm_id":4489679876,"name":"Sugarfire Smokehouse BBQ"}, 
 {"osm_id":4974243121,"name":"Pappy's BBQ"}, 
 {"osm_id":8993197283,"name":"The BBQ & Whiskey Saloon"}, 
 {"osm_id":9244847732,"name":"Bootleggin' BBQ Tavern"}, 
 {"osm_id":9710331857,"name":"Ravin’s BBQ"}]

This works too but probably doesn't give you what you want

						SELECT json_agg( (p.osm_id, p.name )) AS p_json
 FROM st_louis.points AS p
 WHERE name ILIKE '%BBQ%';
					

[{"f1":2553435113,"f2":"Smoking Barrels BBQ"}, 
 {"f1":4489679876,"f2":"Sugarfire Smokehouse BBQ"}, 
 {"f1":4974243121,"f2":"Pappy's BBQ"}, 
 {"f1":8993197283,"f2":"The BBQ & Whiskey Saloon"}, 
 {"f1":9244847732,"f2":"Bootleggin' BBQ Tavern"}, 
 {"f1":9710331857,"f2":"Ravin’s BBQ"}]

Just keep it in the FROM

						SELECT json_agg( p_sub ) AS p_json
 FROM st_louis.points AS p
     CROSS JOIN LATERAL ( SELECT p.osm_id, p.name ) AS p_sub
 WHERE p.name ILIKE '%BBQ%';
					

[{"osm_id":2553435113,"name":"Smoking Barrels BBQ"}, 
 {"osm_id":4489679876,"name":"Sugarfire Smokehouse BBQ"}, 
 {"osm_id":4974243121,"name":"Pappy's BBQ"}, 
 {"osm_id":8993197283,"name":"The BBQ & Whiskey Saloon"}, 
 {"osm_id":9244847732,"name":"Bootleggin' BBQ Tavern"}, 
 {"osm_id":9710331857,"name":"Ravin’s BBQ"}]

ST_AsGeoJSON and json_agg dancing

Together they can create a geojson FeatureCollection.

						SELECT json_build_object(
    'type', 'FeatureCollection',
    'features', json_agg(ST_AsGeoJSON(p)::json)
    )
 FROM stlouis.points AS p
 WHERE name ILIKE '%BBQ%';
					

PostGIS can aggregate rows

Several ST_AsGeoBuf, ST_FlatGeoBuf, ST_AsMVT pg_tileserv is essentially a wrapper around ST_AsMVT. Check out https://flatgeobuf.org/ to learn how to use FlatGeoBuf in your mapping apps.

						SELECT ST_AsFlatGeoBuf(p) AS fgeobuf
 FROM stlouis.points AS p 
WHERE name ILIKE '%BBQ%';
					

Ordered Aggregates

When you care about order, you can apply an ORDER BY to most aggregates

						SELECT array_agg(p ORDER BY name) AS points_by_name
 FROM stlouis.points AS p 
WHERE name ILIKE '%BBQ%';
					

PostGIS Ordered Aggregates

You want your path ordered by time right?

						SELECT name, ST_MakeLine(p.geom ORDER BY time) AS path
 FROM gps_points AS p 
GROUP BY name
ORDER BY name;
					

PostGIS ST_MakeLine also a window aggregate

Or maybe you want to create 2-3 point linestrings consisting of visit point before and visit point after

						SELECT name, 
	ST_MakeLine( p.geom ) 
		OVER( PARTITION BY name ORDER BY time 
			ROWS BETWEEN 1 PRECEDING AND 1 FOLLOWING ) AS path
 FROM gps_points AS p 
ORDER BY name, time;
					

Short-hand casting: The basics

SQL Standard casting syntax

SELECT name, CAST(geom As jsonb) AS geom
FROM st_louis.lines
WHERE name ILIKE 'Crescent%' LIMIT 1;

PostgreSQL casting is so common that there is short-hand for it. Also allows extensions to define casting behavior

SELECT name, geom::jsonb
FROM st_louis.lines
WHERE name ILIKE 'Crescent%' LIMIT 1;

Short-hand cast chaining

If there is no direct cast from one dataype to another you can chain

SELECT name, rast::geometry::geography
FROM dems
WHERE LIMIT 1;

Operators specific to types

What does subtracting a text do for jsonb?

						SELECT osm_id, osm_id - 1 AS id, 
    to_jsonb(p) - 'geom' - 'other_tags' 
		- 'ref' - 'osm_id' - 'man_made' - 'place'
		-  'address' - 'barrier' - 'highway' AS p_json
 FROM st_louis.points AS p
 WHERE name ILIKE '%BBQ%';
					

   osm_id   |     id     |                       p_json
------------+------------+-----------------------------------------------------
 2553435113 | 2553435112 | {"name": "Smoking Barrels BBQ", "is_in": null}
 4489679876 | 4489679875 | {"name": "Sugarfire Smokehouse BBQ", "is_in": null}
 4974243121 | 4974243120 | {"name": "Pappy's BBQ", "is_in": null}
 8993197283 | 8993197282 | {"name": "The BBQ & Whiskey Saloon", "is_in": null}
 9244847732 | 9244847731 | {"name": "Bootleggin' BBQ Tavern", "is_in": null}
 9710331857 | 9710331856 | {"name": "RavinÆs BBQ", "is_in": null}
(6 rows)

The ISO-SQL || operator concatenates text but it works with jsonb too in PostgreSQL to add or update.

						SELECT name || ' Good' AS name2,
    to_jsonb(p) || jsonb_build_object('name', p.name || ' Not Good') AS p_json
 FROM (SELECt osm_id, name 
     FROM stlouis.points) AS p
 WHERE name ILIKE '%BBQ%';
					

            name2             |                               p_json
-------------------------------+---------------------------------------------------------------------
 Smoking Barrels BBQ Good      | {"name": "Smoking Barrels BBQ Not Good", "osm_id": 2553435113}
 Sugarfire Smokehouse BBQ Good | {"name": "Sugarfire Smokehouse BBQ Not Good", "osm_id": 4489679876}
 Pappy's BBQ Good              | {"name": "Pappy's BBQ Not Good", "osm_id": 4974243121}
 The BBQ & Whiskey Saloon Good | {"name": "The BBQ & Whiskey Saloon Not Good", "osm_id": 8993197283}
 Bootleggin' BBQ Tavern Good   | {"name": "Bootleggin' BBQ Tavern Not Good", "osm_id": 9244847732}
 RavinÆs BBQ Good              | {"name": "RavinÆs BBQ Not Good", "osm_id": 9710331857}
(6 rows)					
					

PostGIS geometry / geography operators

KNN distance operator has meaning both in PostGIS and pg_vector extension

Geometry

						SELECT p.name, p.geom <-> ST_Point(-90.188, 38.62658, 4326) AS dist_deg
 FROM st_louis.points AS p
 WHERE name ILIKE '%BBQ%'
 ORDER BY dist_deg LIMIT 6;
					

Geography

						SELECT p.name, p.geog <-> ST_Point(-90.188, 38.62658, 4326)::geography AS dist_m
 FROM st_louis.points AS p
 WHERE name ILIKE '%BBQ%'
 ORDER BY dist_m LIMIT 6;
					

           name           |       dist_m
--------------------------+--------------------
 Sugarfire Smokehouse BBQ | 443.06224064314273
 Pappy's BBQ              | 3274.4233355063307
 RavinÆs BBQ              |  5384.221956799874
 Bootleggin' BBQ Tavern   | 6355.8718317511275
 The BBQ & Whiskey Saloon |  6638.181704382918
 Smoking Barrels BBQ      |    9777.2332126643
(6 rows)
				
					

What's a LATERAL?

For each barbecue place, list the two closest streets to it.

SELECT p.name, s.name AS street, p.geog AS location, s.dist_m
FROM st_louis.points AS p 
 CROSS JOIN LATERAL 
 (SELECT l.name, l.geom,  p.geog <-> l.geog AS dist_m
	FROM st_louis.lines 
	 WHERE highway > '' ORDER BY dist_m 
		LIMIT 2) AS s
 WHERE p.name ILIKE '%BBQ%';

LATERAL keyword can be omitted for functions

SELECT p.name, kv.*
FROM st_louis.points AS p 
	CROSS JOIN LATERAL jsonb_each_text(p.other_tags) AS kv;

Is equivalent to

SELECT p.name, kv.*
FROM st_louis.points AS p, jsonb_each_text(p.other_tags) AS kv;

Function LEFT JOIN laterals

SELECT p.name, kv.*
FROM st_louis.points AS p 
	LEFT JOIN LATERAL jsonb_each_text(p.other_tags) AS kv;

Is equivalent to

SELECT p.name, kv.*
FROM st_louis.points AS p LEFT JOIN jsonb_each_text(p.other_tags) AS kv ON true;

Explode and Implode ala PostGIS + raster

Chaining of functions in a lateral way
PostgreSQL / GDAL / PostGIS bump

SELECT d.district, gv.val AS elev_ft
/** PostGIS aggregate function, 
 will collect so one geometry per pixel value per district**/				
 , ST_Union(gv.geom) As geom
FROM dem
 INNER JOIN districts AS d ON ST_Intersects(dem.rast, d.geom)
 /** not a set returning function, returns 
  one row per dem.rast,geom intersection, 
	only consider elevation higher than 150 m and convert to ft**/
     , ST_Reclass(ST_Clip(dem.rast,d.geom,1),
                '150-400:328-1312', '16BUI' ) AS re(rast)
 /** set returning function returns at 
  least one row per pixel value per district **/
	 , ST_DumpAsPolygons( re.rast ) AS gv 
GROUP BY d.district, gv.val;

Use Ordinality with set returning functions

Use Ordinality to number results from set returning functions such as unnest, PostGIS ST_SubDivide that don't give you numbers.

SELECT town,  f.ord, ST_NPoints(f.geom)
FROM shps.towns_polym, ST_SubDivide(geom, 100) WITH ordinality f(geom,ord)
WHERE town IN('BOSTON', 'CAMBRIDGE');
            	

Before had 2 rows town | st_npoints -----------+------------ BOSTON | 1893 CAMBRIDGE | 235 (2 rows)

After have 68 rows, no geometry has more than 100 points town | ord | st_npoints -----------+-----+------------ BOSTON | 1 | 89 BOSTON | 2 | 62 : BOSTON | 22 | 97 : BOSTON | 64 | 6 CAMBRIDGE | 1 | 40 : CAMBRIDGE | 4 | 63 (68 rows)

Extension mixins

You'll often see PostGIS extensions paling around with these

• ogr_fdw
• pgrouting
• http
• postgres_fdw
• h3
• mobilitydb
• and the list goes on and on

Data wrangling with ogr_fdw

If you have all sorts of data of both a spatial and non-spatial flavor to tame, make sure you have ogr_fdw foreign data wrapper in your tool belt.

• For windows users, it's part of PostGIS bundle spatial extensions on application stackbuilder.
• For CentOS/Red Hat/Scientific etc, it's available via yum.postgresql.org
• For Debian based systems (Debian, Ubuntu) etc, it's available via apt.postgresql.org
• For others, if you have PostGIS with GDAL support, just need postgresql dev package to compile. Download the source https://github.com/pramsey/pgsql-ogr-fdw

Why is ogr_fdw so hot?

You have the combined power of GDAL, PostgreSQL, and any PostgreSQL extension you want (including PostGIS) working seamlessly together. So many kinds of data you can query and take advantage of PostgreSQL functions and any extension functions and types such as PostGIS, hstore, built-in json and jsonb to tame your data.

• Spreadsheets
• ODBC datasources
• OSM files (OSM, PBF)
• ESRI Shapefiles
• Spatial web services
• Many more

Enable it in your database

Binaries available at least on OSGeoLive, Windows PostGIS Bundle, apt.postgresql.org, yum.postgresql.org

CREATE EXTENSION ogr_fdw;

PostgreSQL + GDAL (OGR) ~ PostGIS = OGR_FDW
PostgreSQL Foreign Data Wrapper

Doesn't require PostGIS to use, but will expose spatial columns as PostGIS geometry if PostGIS is installed.

I call this the PostgreSQL/OGR/PostGIS bump.

Check what formats you can read

SELECT name  
	FROM unnest(ogr_fdw_drivers()) AS f(name) ORDER BY f.name;

name
----------------
AmigoCloud
:
CAD
:
CSV
:
FlatGeobuf
Geoconcept
GeoJSON
GeoJSONSeq
GeoRSS
GML
GMLAS
GPKG
GPSBabel
GPX
GTFS
HTTP
:
JML
JP2OpenJPEG
KML
:
MapInfo File
MapML
MBTiles
Memory
MSSQLSpatial
MVT
MySQL
:
ODBC
ODS
:
WFS
XLS
XLSX
(80 rows)

Link in a whole folder of ESRI Shapefiles and Dbase files

CREATE SERVER svr_shp FOREIGN DATA WRAPPER ogr_fdw
OPTIONS (datasource 'C:/fdw_data/massgis/shps',
 format 'ESRI Shapefile'
);
CREATE SCHEMA shps;
IMPORT FOREIGN SCHEMA ogr_all
FROM SERVER svr_shp INTO shps;

\dE shps.*

                    List of relations
 Schema |        Name         |     Type      |  Owner
--------+---------------------+---------------+----------
 shps   | biketrails_arc      | foreign table | postgres
 shps   | towns_arc           | foreign table | postgres
 shps   | towns_poly          | foreign table | postgres
 shps   | towns_poly_areacode | foreign table | postgres
 shps   | towns_polym         | foreign table | postgres
 shps   | towns_pop           | foreign table | postgres
 shps   | zipcodes_nt_poly    | foreign table | postgres
(7 rows)

Query your geometry_columns catalog

SELECT f_table_name As tbl, f_geometry_column As geom, srid, type
FROM geometry_columns 
WHERE f_table_schema = 'shps'
ORDER BY tbl;

       tbl        | geom | srid    |    type
------------------+------+--------+------------
 biketrails_arc   | geom |  26986 | LINESTRING
 towns_arc        | geom |  26986 | LINESTRING
 towns_poly       | geom |  26986 | POLYGON
 towns_polym      | geom |  26896 | POLYGON
 zipcodes_nt_poly | geom |  26986 | POLYGON
(5 rows)

You can fix bad geometries right in shape file with power of PostGIS

Make sure the user that postgres runs under has edit/delete rights to the folder holding the shape files.

UPDATE shps.towns_polym
    SET geom = ST_MakeValid(geom)
WHERE NOT ST_IsValid(geom)
RETURNING town;

NOTICE:  Ring Self-intersection at or near point 241494.43330000341 890709.87110000104
NOTICE:  Ring Self-intersection at or near point 306590.87370000035 822452.56080000103
NOTICE:  Ring Self-intersection at or near point 273304.93349999934 802752.31069999933

Total query runtime: 320 msec
town
-----
QUINCY
YARMOUTH
TISBURY

Web Files

If your GDAL is compiled with curl support (you see HTTP in list of driver options) , you can link files from the web using the /vsicurl switch. Use /vsizip//vsicurl if they are zipped as well.

CREATE SERVER svr_osm
   FOREIGN DATA WRAPPER ogr_fdw
  OPTIONS (datasource 
 '/vsizip//vsicurl/https://postgis.us/downloads/foss4gna2024/ST_Louis_MO.osm.zip',
	format 'OSM',
   updateable 'false');
   
CREATE SCHEMA IF NOT EXISTS osm;
IMPORT FOREIGN SCHEMA ogr_all 
FROM SERVER svr_osm INTO osm;
SELECT f_table_name, f_geometry_column, srid, type 
	FROM geometry_columns WHERE f_table_schema = 'osm';

  f_table_name   | f_geometry_column | srid |        type
------------------+-------------------+------+--------------------
 points           | geom              | 4326 | POINT
 lines            | geom              | 4326 | LINESTRING
 multilinestrings | geom              | 4326 | MULTILINESTRING
 multipolygons    | geom              | 4326 | MULTIPOLYGON
 other_relations  | geom              | 4326 | GEOMETRYCOLLECTION
(5 rows)

-- requires CREATE EXTENSION hstore; to cast other_tags to hstore 
-- and hstore extension has function hstore_to_jsonb that will cast hstore to jsonb
-- 56789 rows affected, 15.6 secs execution time.
-- 3602 is epsg stateplane meters missouri east
CREATE TABLE st_louis_polys AS 
SELECT osm_id, name, geom::geography As geog, 
	ST_Transform(geom,3602)::geometry(MULTIPOLYGON,3602) AS geom,
    other_tags::hstore::jsonb AS other_tags, building
FROM osm.multipolygons;

Even spreadsheets

Each workbook is considered a server and each sheet a table

CREATE SERVER svr_fedex FOREIGN DATA WRAPPER ogr_fdw
OPTIONS (datasource 'C:/fdw_data/Fedex2016.xls',
 format 'XLS'
);
-- link only 1 spreadsheet preserve headers
IMPORT FOREIGN SCHEMA ogr_all  LIMIT TO (Fedex_Rates_IP)
	FROM SERVER svr_fedex INTO public OPTIONS (launder_column_names 'false');
	
SELECT * FROM fedex_rates_ip;

 fid |     Type     | Weight | Zone A  | Zone B  | Zone C  | Zone D  | Zone E  | Zone F  | Zone G  | Zone H  | Zone I  | Zone J  | Zone K  | Zone L  | Zone M  | Zone N  | Zone O  | Zone Puerto Rico
-----+--------------+--------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+------------------
   2 | IntlPriority |      0 |   40.25 |    41.5 |      43 |   54.75 |   116.5 |      52 |    51.5 |   61.25 |    49.5 |    63.5 |    62.5 |      82 |  119.25 |      61 |   63.25 |            36.65
   3 | IntlPriority |     -1 |   66.25 |   67.75 |   62.25 |   74.25 |     132 |      68 |   68.25 |   85.75 |   66.25 |    84.5 |   82.25 |    99.5 |   136.5 |   79.75 |    85.5 |
   4 | IntlPriority |     -2 |   70.25 |    73.5 |   65.75 |   77.25 |  156.25 |      73 |   73.25 |   87.25 |    73.5 |    90.5 |   87.75 |  129.25 |  160.75 |      86 |      96 |

-- unpivot a subset of columns and keep others 
WITH fkv AS (
SELECT f."Type" As type, f."Weight" As weight,  
	each(row_to_jsonb(f) - '{fid,Type,Weight}'::text[]) AS kv 
from fedex_rates_ip AS f)
SELECT type, weight, (kv).key AS zone, (kv).value::numeric As price
FROM fkv;

     type     | weight |       zone       |  price
--------------+--------+------------------+---------
 IntlPriority |      0 | Zone A           |   40.25
 IntlPriority |      0 | Zone B           |    41.5
 IntlPriority |      0 | Zone C           |      43
 IntlPriority |      0 | Zone D           |   54.75
:

Even CSV files

You can point at a single CSV file or a whole folder of CSV files. Each file is considered a table.

Folder of CSV files

CREATE SERVER svr_census FOREIGN DATA WRAPPER ogr_fdw
OPTIONS (datasource 'C:/fdw_data/census',
 format 'CSV'
);

IMPORT FOREIGN SCHEMA ogr_all
FROM SERVER svr_census INTO public;

Single file

CREATE SERVER svr_census_income FOREIGN DATA WRAPPER ogr_fdw
OPTIONS (datasource '/fdw_data/census/income.csv',
 format 'CSV'
);

IMPORT FOREIGN SCHEMA ogr_all
FROM SERVER svr_census_income INTO public;

Even other relational databases

Format for SQL Server 'ODBC:your_user/your_password@yourDSN,table1,table2'. ODBC can be slow with a lot of tables (more than 150) so filter list if you have over 200 tables

CREATE SERVER svr_sqlserver FOREIGN DATA WRAPPER ogr_fdw
OPTIONS (datasource 'ODBC:pguser/whatever@MSSQLTest,dbo.IssueLog,dbo.IssueNotes',
 format 'ODBC'
);
CREATE SCHEMA IF NOT EXISTS ss;
IMPORT FOREIGN SCHEMA "dbo." 
	FROM SERVER svr_sqlserver INTO ss;

\dE ss.*

                 List of relations
 Schema |      Name      |     Type      |  Owner
--------+----------------+---------------+----------
 ss     | dbo_issuelog   | foreign table | postgres
 ss     | dbo_issuenotes | foreign table | postgres
(2 rows)

http extension: stupid idea or pure genius

pgsql-http now availabe on apt.postgresql.org

Debian / Ubuntu: Replace 17 with your PostgreSQL version

apt install postgresql-17-http

CREATE EXTENSION http;

Embrace the web without having to get out of your database.

Straight from the web to your database

Note the use of the PostGIS ST_GeomFromGeoJSON working in conjunction with the http_get function packaged in http extension as well as built-in PostgreSQL json functions.

CREATE TABLE community_improvement_districts AS 
SELECT (je->>'id')::integer AS id, 
	je->'properties'->>'Name' AS name, 
	ST_GeomFromGeoJSON(je->>'geometry') AS geom 
FROM http_get(
	'https://www.stlouis-mo.gov/data/upload/data-files/CID.geojson'
	) AS h, 
	jsonb_array_elements(h.content::jsonb->'features') AS je;

FIN

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