### ### Planet PostGIS

# Welcome to Planet PostGIS

## May 23, 2024

### Crunchy Data

#### Converting DMS to PostGIS Point Geometry

I love taking random spatial data and turning it into maps. Any location data can be put into PostGIS in a matter of minutes. Often when I’m working with data that humans collected, like historic locations or things that have not yet traditionally been done with computational data, I’ll find traditional Degrees, Minutes, Seconds (DMS) data. To get this into PostGIS and QGIS, you’ll need to convert this data to a different system for decimal degrees. There’s probably proprietary tools that will do this for you, but we can easily write our own code to do it. Let’s walk through a quick example today.

Let’s say I found myself with a list of coordinates, that look like this:

``````38°58′17″N 95°14′05″W
``````

(this is the location of my town’s haunted hotel 👻)

This format of writing geographic coordinates is called DMS, Degrees, Minutes, Seconds (DMS). If you remember from 4th grade geography lessons, that is the latitude on the left there, representing N or S of the equator and longitude East or West of the Prime Meridian.

## WKT & XY coordinates

PostGIS, and most computational spatial systems, work with a geographic system that is akin to an XY grid of the entire planet. Because it is XY, it is a longitude, latitude (X first) system.

PostGIS utilizes with two kinds of geometry values:

• WKT (Well-known text) where a point would look like this `POINT(-126.4 45.32)`
• WKB (Well-known binary) where a point would look like this `0101000000000000000000F03F000000000000F03`

Most often you’ll see the binary used to represent stored data and you can use a function, `st_astext`, to view or query it as text.

## Converting coordinates to decimal degrees

To convert our traditional coordinates into decimals or WKT, we can use decimal math like this:

``````({long_degree}+({long_minutes}/60)+({long_seconds}/3600)
``````

So for our location:

``````-- starting location
38°58′17″N 95°14′05″W

-- formula
38+(58/60)+(17/3600), 95+(14/60)+(05/3600)

-- switch the order since this is X first
-- make the Western quad negative
-- getting this result

-95.2472222, 38.9713888
``````

### Regex Function for Making PostGIS Points out of DMS

If you have one location like this, you probably have a lot, so we’ll need a more sophisticated solution for our whole data set. You know if you need something done right, you ask Paul Ramsey. Paul worked with me on getting this function written that will convert DMS to PostGIS friendly (binary geometry) point data.

``````CREATE OR REPLACE FUNCTION dms_to_postgis_point(dms_text TEXT)
RETURNS geometry AS
\$\$
DECLARE
dms TEXT[] := regexp_match(dms_text, '(\d+)\D+(\d+)\D+(\d+)\D+([NS])\D+(\d+)\D+(\d+)\D+(\d+)\D+([EW])');
lat float8;
lon float8;
BEGIN
lat := dms[1]::float8 + dms[2]::float8/60 + dms[3]::float8/3600;
lon := dms[5]::float8 + dms[6]::float8/60 + dms[7]::float8/3600;
IF upper(dms[4]) = 'S' THEN
lat := -1 * lat;
END IF;

IF upper(dms[8]) = 'W' THEN
lon := -1 * lon;
END IF;

RETURN ST_Point(lon, lat, 4326);
END;
\$\$
LANGUAGE 'plpgsql'
IMMUTABLE
STRICT;
``````

Let’s do a quick test with our original point:

``````SELECT st_astext(dms_to_postgis_point('38°58′17″N 95°14′05″W'));

st_astext
---------------------------------------------
POINT(-95.23472222222222 38.97138888888889)
(1 row)
``````

Great, that works.

### Creating a new column with your geometry

Now we can use built-in PostGIS functions to add a new geom column and run the function on our old lat_long column.

``````ALTER TABLE my_table ADD COLUMN geom geometry(Point);

UPDATE my_table SET geom = dms_to_postgis_point(lat_long);
``````

### Conclusion

PostGIS is just packed with so many cool functions to make sure you can turn anything into maps. Hope this helps you get started if you’re using traditional lat long data.

## March 19, 2024

### Crunchy Data

#### Inside PostGIS: Calculating Distance

Calculating distance is a core feature of a spatial database, and the central function in many analytical queries.

• "How many houses are within the evacuation radius?"
• "Which responder is closest to the call?"
• "How many more miles until the school bus needs routine maintenance?"

PostGIS and any other spatial database let you answer these kinds of questions in SQL, using ST_Distance(geom1, geom2) to return a distance, or ST_DWithin(geom1, geom2, radius) to return a true/false result within a tolerance.

``````SELECT ST_Distance(
'LINESTRING (150 300, 226 274, 320 280, 370 320, 390 370)'::geometry,
'LINESTRING (140 180, 250 230, 350 200, 390 240, 450 200)'::geometry
);
``````

It all looks very simple, but under the covers there is a lot of machinery around getting a result fast for different kinds of inputs.

## Distance Under the Covers

Distance should be easy! After all, we learn how to calculate distance in middle school! The Pythagorean Theorem tells us that the square of the hypotenuse of a right triangle is the sum of the squares of the two other sides.

So, problem solved, right?

Not so fast. Pythagorus gives us the distance between two points, but objects in spatial databases like PostGIS can be much more complex.

How would I calculate the distance between two complex polygons?

## Brute Force

The straight-forward solution is to just find the distance between every possible combination of edges in the two polygons, and return the minimum of that set.

This is a "quadratic" algorithm, what computer scientists call `O(n^2)`, because the amount of work it generates is proportional to the square of the number of inputs. As the inputs get big, the amount of work gets very very very big.

Fortunately, there are better ways.

## Projection and Pruning

The distance implementation in PostGIS has two major code paths:

• For disjoint (non-overlapping) inputs, an optimized calculation; and,
• For overlapping inputs, the brute force calculation.

Disjoint inputs are handled with a clever simplification of the problem space. Because the inputs are disjoint, it is possible to construct a line between the centers of the two inputs.

If every edge in each object is projected down onto the line, it becomes possible to perform a sort of those edges, such that edges that are near on the line are also near in the sorted lists, and near in space.

Starting from the mid-point of each object it is relatively inexpensive to quickly prune away large numbers of edges that are definitely not the nearest edges, leaving a much smaller number of potential targets that need to have their distance calculated.

The cost of creating the projected segments is just `O(n)`, but the cost of the sort step is `O(n*log(n))` so the overall cost of the algorithm is `O(n*log(n))`.

This is all well and good, but what if the inputs do overlap? Then the algorithm falls back to brute-force and `O(n^2)`. Is there any way to avoid that?

## Linear Time Spatial Trees

The project-and-prune approach is very clever, but it is possible to generate a spatially searchable representation of the edges even faster, by using the fact that edges in a LineString or LinearRing are highly spatial autocorrelated:

• The end point of one edge is always the start point of the next.
• The edges mostly don't cross each other.

Basically, the edges are already spatially pre-sorted. That means it is possible to build a decent tree structure from them incurring any non-linear computational cost.

Start with the edges in sorted order. The bounds of the edges form the leaf nodes of a spatial tree. Merge neighboring leaf nodes, now you have the first level of interior nodes. Continue until you have only one node left, that is your root node. The cost is `O(n) + O(0.5n) + O(0.25n) ...` which is to say in aggregate, `O(n)`.

Ordinarily, building a spatial tree would be expected to cost about `O(n*log(n))`, so this is a nice win.

The `CIRC_NODE` tree used to accelerate distance calculation for the `geography` type is built using this process.

## Overlapping Inputs and Distance Calculation

There is no guarantee that a tree-indexed approach will crack the overlapping polygon problem.

Disjoint polygons are very amenable to distance searching trees, because it is easy to discard whole branches of the tree that are definitionally too far away to contain candidate edges.

As inputs begin to overlap, it becomes harder to discard large portions of the trees, and as a result a lot of computation is spent traversing the tree, even if a moderate proportion of candidates can be discarded from the lower branches of the tree.

## Next Steps

The distance calculation in PostGIS has not been touched in many years, for good reason: it's really important, so any re-write has to be definitely an improvement on the existing code, over all known (and unknown) use cases.

However, there is some already built and tested code, in the code base, which has never been turned on, the RECT_TREE.

Like the `CIRC_NODE` tree in geography, this implementation is based on building a tree from spatially coherent inputs. Unlike the `CIRC_NODE` tree, it has not been proven to be faster than the existing implementation in all cases.

A next development step will be to revive this implementation, evaluate it for implementation efficiency, and test effectiveness:

• Can it exceed the current sort-and-prune strategy for disjoint polygons?
• Can it exceed brute-force for overlapping polygons?

## March 06, 2024

### Crunchy Data

#### Connecting QGIS to Postgres and PostGIS

QGIS, the Quantum Geographic Information System, is an open-source graphical user interface for map making. QGIS works with a wide variety of file types and has robust support for integrating with Postgres and PostGIS. Today I just wanted to step through getting QGIS connected to a Postgres database and the basic operations that let you connect the two systems.

## Connecting QGIS to Postgres

Connecting QGIS to Postgres is very similar to any other GUI or application, you’ll need the database host, login, and password details. This is the same process for a local connection or remote one (like Crunchy Bridge). You’ll connect the first time through the Browser option listed PostgreSQL and Add New Connection.

By default, QGIS will store your passwords as plain text in a file. If you’re just working with a local database and don’t have anything special in there, that may not be a problem. But if you’re working with a larger production database shared by lots of users, you’ll want to opt for a higher level of protection for the password. In your PostgreSQL connections box, you’ll see a way to add Configurations for the password. Here you can create a master password, store your database credentials, and they’ll be encrypted and only decrypted with your master password.

## Using QGIS to load data

QGIS is a great way to get spatial data into Postgres and PostGIS. You can use any file type supported by QGIS including vector types like shapefiles (shp), GeoJSON, and even csv files on your local machine. To load data into QGIS, you’ll first go to Layer —> Add Layer and choose the type of file you have.

For this sample I have a county map of the state of Kansas. Maps like this are often freely available for download from government agencies.

Once my layer is in, I can toggle on to show labels which will add any label data for your geometry.

Now that I have data in QGIS I can save this to a Postgres database. This will allow me to work with this data later. Go to the DB manager icon, and choose Import Layer.

There are several settings here, like choosing the primary key, the origin and destination SRID. QGIS will even suggest that adding an index for your geometry column is a good idea and will build an index in your database for you.

QGIS works both ways, so if you already have a dataset to work from, you can just use that data as your source. In that case, you’ll start from the Layer — Add Layer option. You’ll either need to specify the database connection you want, or add a new one here. You’ll be able to open all the tables in your database and choose which ones to add as a layer in your map viewer.

Depending on the file you get, QGIS may or may not be super happy with it. You might see a warning icon next to your file. The main issues that QGIS will be warning you about are:

• There’s no spatial reference id

The spatial reference ID is an important quirk when dealing with geospatial data. You can default to 4326 if you don’t have a better option.

To find a spatial reference id, or see if it is set:

``````SELECT ST_SRID(geom) FROM my_table_name LIMIT 1;
``````

And to update it:

``````SELECT UpdateGeometrySRID('my_table_name','geom',4326);
``````
• There’s no geometry column

Assuming you have points, lines, or polygon data and it is just in the wrong data type, you can create a new column for the geometry and point your data to that new column.

``````ALTER TABLE my_table_name

UPDATE my_table_name
SET geom = ST_SetSRID(ST_MakePoint(my_column, my_column), 4326);
``````
• There’s no primary key

Relational databases rely on a primary key to tie other data together. If there’s already an id column, you can just create a primary key index on it like this. If there’s not a unique column like that, you might have to do a bit more work on the data to get this fixed.

``````alter table my_table add primary key (id_column);
``````

## Writing and Saving SQL

One really cool feature of QGIS is that you can write SQL directly against your Postgres database and view the results as spatial geometries. You can save queries for use later as well. You can also use QGIS to create a layer based on your query results.

Here’s a sample query where I’m joining two data sources, one my geometry of Kansas counties that I loaded earlier. Second population data by county. I’m selecting just the geometry column and with the load option can have QGIS add that query result as a layer.

QGIS will also let you save a query as “view”. This is a database specific term that will save your query results as a table for use later. This can also be loaded as a layer in QGIS projects. There’s an overview of using views in the post Postgres Subquery Powertools. Views are a great idea if you’re using only a small subset of data in your QGIS map but you are storing a larger dataset as well.

Here’s an example of a map I made using 4 different query layers, one for each different population density.

Don’t forget that QGIS works in stacked layers, so your new SQL query layers will have to be on top of your base map or they won’t be visible.

## Saving QGIS projects in Postgres

You can also save your QGIS project in your Postgres database. This is under the Project — Save to options. This can be a good idea if you want others on your team to have access to projects or if you don’t want your QGIS projects stored locally.

## Final notes

There’s a great video of this and more from PostGIS Day 2020 called QGIS and PostGIS.

• You can use QGIS to load shape files and other file types into Postgres
• You can use QGIS to create maps from existing Postgres/PostGIS data sources
• You can write queries in QGIS against your Postgres data and show the results as geometry layers
• You can write join queries in QGIS and join your geometry fields with other attribute data or tables in your database
• You can save all of your project work for QGIS in the Postgres database

## February 13, 2024

### Crunchy Data

#### PostGIS Clustering with K-Means

Clustering points is a common task for geospatial data analysis, and PostGIS provides several functions for clustering.

We previously looked at the popular DBSCAN spatial clustering algorithm, that builds clusters off of spatial density.

This post explores the features of the PostGIS ST_ClusterKMeans function. K-means clustering is having a moment, as a popular way of grouping very high-dimensional LLM embeddings, but it is also useful in lower dimensions for spatial clustering.

ST_ClusterKMeans will cluster 2-dimensional and 3-dimensional data, and will also perform weighted clustering on points when weights are provided in the "measure" dimension of the points.

## Some Points to Cluster

To try out K-Means clustering we need some points to cluster, in this case the 1:10M populated places from Natural Earth.

``````ogr2ogr \
-f PostgreSQL \
-nln popplaces \
-lco GEOMETRY_NAME=geom \
PG:'dbname=postgres' \
ne_10m_populated_places_simple.shp
``````

## Planar Cluster

A simple clustering in 2D space looks like this, using 10 as the number of clusters:

``````CREATE TABLE popplaces_geographic AS
SELECT geom, pop_max, name,
ST_ClusterKMeans(geom, 10) OVER () AS cluster
FROM popplaces;
``````

Note that pieces of Russia are clustered with Alaska, and Oceania is split up. This is because we are treating the longitude/latitude coordinates of the points as if they were on a plane, so Alaska is very far away from Siberia.

For data confined to a small area, effects like the split at the dateline do not matter, but for our global example, it does. Fortunately there is a way to work around it.

## Geocentric Cluster

We can convert the longitude/latitude coordinates of the original data to a geocentric coordinate system using ST_Transform. A "geocentric" system is one in which the origin is the center of the Earth, and positions are defined by their X, Y and Z distances from that center.

In a geocentric system, positions on either side of the dateline are still very close together in space, so it's great for clustering global data without worrying about the effects of the poles or date line. For this example we will use EPSG:4978 as our geocentric system.

Here are the coordinates of New York, converted to geocentric.

``````SELECT ST_AsText(ST_Transform(ST_PointZ(74.0060, 40.7128, 0, 4326), 4978), 1);
``````
``````POINT Z (1333998.5 4654044.8 4138300.2)
``````

And here is the cluster operation performed in geocentric space.

``````CREATE TABLE popplaces_geocentric AS
SELECT geom, pop_max, name,
ST_ClusterKMeans(
ST_Transform(
ST_Force3D(geom),
4978),
10) OVER () AS cluster
FROM popplaces;
``````

The results look very similar to the planar clustering, but you can see the "whole world" effect in a few places, like how Australia and all the islands of Oceania are now in one cluster, and how the dividing point between the Siberia and Alaska clusters has moved west across the date line.

It's worth noting that this clustering has been performed in three dimensions (since geocentric coordinates require an X, Y and Z), even though we are displaying the results in two dimensions.

## Weighted Cluster

In addition to naïve k-means, ST_ClusterKMeans can carry out weighted k-means clustering, to push the cluster locations around using extra information in the "M" dimension (the fourth coordinate) of the input points.

Since we have a "populated places" data set, it makes sense to use population as a weight for this example. The weighted algorithm requires strictly positive weights, so we filter out the handful of records that are non-positive.

``````CREATE TABLE popplaces_geocentric_weighted AS
SELECT geom, pop_max, name,
ST_ClusterKMeans(
ST_Force4D(
ST_Transform(ST_Force3D(geom), 4978),
mvalue => pop_max
),
10) OVER () AS cluster
FROM popplaces
WHERE pop_max > 0;
``````

Again, the differences are subtle, but note how India is now a single cluster, how the Brazil cluster is now biased towards the populous eastern coast, and how North America is now split into east and west.

## February 08, 2024

### PostGIS Development

#### PostGIS Patch Releases

The PostGIS development team is pleased to provide bug fix releases for 3.4.2, 3.3.6, 3.2.7, 3.1.11, 3.0.11, and 2.5.11. Please refer to the links above for more information about the issues resolved by these releases.

## February 06, 2024

### Paul Ramsey

#### Building the PgConf.Dev Programme

Update: The programme is now public.

The programme for pgconf.dev in Vancouver (May 28-31) has been selected, the speakers have been notified, and the whole thing should be posted on the web site relatively soon.

I have been on programme committees a number of times, but for regional and international FOSS4G events, never for a PostgreSQL event, and the parameters were notably different.

The parameter that was most important for selecting a programme this year was the over 180 submissions, versus the 33 available speaking slots. For FOSS4G conferences, it has been normal to have between two- and three-times as many submissions as slots. To have almost six-times as many made the process very difficult indeed.

Why only 33 speaking slots? Well, that’s a result of two things:

• Assuming no more than modest growth over the last iteration of PgCon, puts attendence at around 200, which is the size of our plenary room. 200 attendees implies no more than 3 tracks of content.
• Historically, PostgreSQL events use talks of about 50 minutes in length, within a one hour slot. Over three tracks and two days, that gives us around 33 talks (with slight variations depending on how much time is in plenary, keynotes or lightning talks).

The content of those 33 talks falls out from being the successor to PgCon. PgCon has historically been the event attended by all major contributors. There is an invitation-only contributors round-table on the pre-event day, specifically for the valuable face-to-face synch-up.

Given only 33 slots, and a unique audience that contains so many contributors, the question of what pgconf.dev should “be” ends up focussed around making the best use of that audience. pgconf.dev should be a place where users, developers, and community organizers come together to focus on Postgres development and community growth.

That’s why in addition to talks about future development directions there are talks about PostgreSQL coding concepts, and patch review, and extensions. High throughput memory algorithms are good, but so is the best way to write a technical blog entry.

Getting from 180+ submissions to 33 selections (plus some stand-by talks in case of cancellations) was a process that consumed three calls of over 2 hours each and several hours of reading every submitted abstract.

The process was shepherded by the inimitable Jonathan Katz.

• A first phase of just coding talks as either “acceptable” or “not relevant”. Any talks that all the committee members agreed was “not relevant” were dropped from contention.
• A second phase where each member picked 40 talks from the remaining set into a kind of “personal program”. The talks with just one program member selecting them were then reviewed one at a time, and that member would make the case for them being retained, or let them drop.
• A winnow looking for duplicate topic talks and selecting the strongest, or encouraging speakers to collaborate.
• A third “personal program” phase, but this time narrowing the list to 33 talks each.
• A winnow of the most highly ranked talks, to make sure they really fit the goal of the programme and weren’t just a topic we all happened to find “cool”.
• A talk by talk review of all the remaining talks, ensuring we were comfortable with all choices, and with the aggregate make up of the programme.

The programme committee was great to work with, willing to speak up about their opinions, disagree amicably, and come to a consensus.

Since we had to leave 150 talks behind, there’s no doubt lots of speakers who are sad they weren’t selected, and there’s lots of talks that we would have taken if we had more slots.

If you read all the way to here, you must be serious about coming, so you need to register and book your hotel right away. Spaces are, really, no kidding, very limited.

## January 30, 2024

### Crunchy Data

#### JSON and SVG from PostGIS into Google Sheets

At PostGIS Day 2023, one of our speakers showed off a really cool demo for getting JSON and SVGs in and out of Postgres / PostGIS and into Google Sheets. Brian Timoney put together several open source projects in such a cool way that I just had to try it myself. If you want to see his demo video, it is on YouTube. With Brian’s blessing, I’m writing up some additional details with a few of the sample code bits for those of you that want to try this or something similar for your own projects.

So what do we have here? We have Postgres data that is coming into Google sheets in real time, tied to a custom SVG.

Before we dive in, an overview of things that I’ll cover to make this happen :

• Running pg_featureserv from inside a database
• Serving JSON data from Postgres with pg_featureserv
• Functions to create SVG imasges
• Serving SVG images from Postgres with pg_featureserv
• JSON & SVGs delivered to Google Sheets

## pg_svg functions

Brian wrote some special stuff for the goal and shot data in his examples to get it just right so if you want to play with a sample of that, here’s some data for you.

Sample data
``````-- Regular season goals of Colorado Avalanche for 2021-2022 season
--

-- information derived from public NHL API
-- shot locations normalized for offensive end of ice
-- hex_id was derived for a demo for the 2022 PostGIS Day demo

SET client_encoding = 'UTF8';

CREATE TABLE public.goals (
rec_num integer,
game_id text,
this_event text,
this_event_code text,
players text,
playerid text,
player_role text,
x_coord text,
y_coord text,
team_scored text,
period_num text,
period_type text,
plot_x numeric,
plot_y numeric,
plot_pt public.geometry(Point,32613),
hex_id numeric
);

ALTER TABLE public.goals OWNER TO postgres;

INSERT INTO public.goals VALUES (26, '2021020005', 'Goal', 'COL24', 'Jack Johnson', '8471677', 'Scorer', '-77.0', '1.0', 'Colorado Avalanche', '1', 'REGULAR', 77.0, -1.0, '01010000000000000000405340000000000000F0BF', 258);
INSERT INTO public.goals VALUES (5997, '2021020982', 'Goal', 'SJS204', 'Darren Helm', '8471794', 'Scorer', '-76.0', '5.0', 'Colorado Avalanche', '1', 'REGULAR', 76.0, -5.0, '0101000000000000000000534000000000000014C0', 257);
INSERT INTO public.goals VALUES (2432, '2021020415', 'Goal', 'COL45', 'Darren Helm', '8471794', 'Scorer', '-75.0', '-3.0', 'Colorado Avalanche', '1', 'REGULAR', 75.0, 3.0, '01010000000000000000C052400000000000000840', 258);
``````
Sample data for the rink (this is a DXF CAD file encoded as a geometry string)
``````  -- Rink map created from random DXF file found on internet
--
-- ** IMPORTANT: SRID of 32613 is fake!! It is just an arbitrary X/Y plane

SET client_encoding = 'UTF8';

CREATE TABLE public.therink ( id_num integer, geom
public.geometry(MultiLineString,32613) );

ALTER TABLE public.therink OWNER TO postgres;

INSERT INTO public.therink VALUES (1,

``````

The special visual functions you will see with the shots and the rink are done with some functions for Postgres in pg_svg. This isn’t really an extension, it is just a handy functions you can load into any Postgres database to help create SVGs. To load this in your database download it and run something like:

``````psql postgres://postgres:123155012sdfxxcwsdfweweff@p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5432/postgres < pg-svg-lib.sql
``````

Brian wrote two sample SVG shot chart functions for this project. One large and one small that could fit inside a spreadsheet if that’s your final data destination.

Big chart function
``````-- FUNCTION: postgisftw.big_chart(text)

-- pg_featureserv looks for functions in the 'postgisftw' schema

-- ** IMPORTANT  SRID 32613 is 'fake', just needed a planar projection to work with the arbitrary X/Y of the rink

CREATE OR REPLACE FUNCTION postgisftw.big_chart(
player_id text DEFAULT '8476455'::text) -- if no id, then Landeskog goals
RETURNS TABLE(svg text)
LANGUAGE 'plpgsql'
COST 100
STABLE STRICT PARALLEL UNSAFE
ROWS 1000

AS \$BODY\$
BEGIN
RETURN QUERY

-- we only want to display the offensive end of the rink
with half_rink as (
select st_intersection(therink.geom,ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from therink
),

goals as (
-- collect all of a player's goals into a single geometry
select 	ST_SetSRID(st_collect(geom)::geometry, 32613) as geom
from
(	SELECT
st_intersection(ST_SetSRID(goals.plot_pt,32613),ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from postgisftw.goals
WHERE playerid ILIKE player_id || '%'

)q1
),
-- get player name, team total number of goals for title display
playerinfo AS (

select UPPER(q1.players) as this_player,min(q1.team_scored) as this_team
, count(q1.playerid)::text as num_goals from
(
select a.players,a.playerid,a.team_scored
from postgisftw.goals a WHERE playerid ILIKE player_id || '%'
) q1
group by q1.players

),

-- make the SVG document

shapes AS (

--rink styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#2E9AFE',
'stroke-width', 0.5::text )
)
svg FROM half_rink

UNION ALL
-- goals styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#F5A9A9',
'stroke-width', 0.5::text,
'fill','#FF0040'
)
)
svg FROM goals

UNION ALL
-- player name, team, and total number of goals
-- create an arbitrary point underneath the rink and reasonably centered
SELECT NULL, svgText(ST_SetSRID(ST_MakePoint(50,-46),32613),this_player||'  ('||this_team||') -  '||num_goals||' goals',
style => svgStyle(  'fill', '#585858', 'text-anchor', 'middle', 'font', 'bold 3px sans-serif' ) )
svg
from playerinfo

)
-- create the final viewbox
-- Martin Davis uses a sensible default of expanding the collected geometry
-- since our rink geometry is static, I hard coded a viewbox
SELECT svgDoc( array_agg( shapes.svg ),
viewbox => '-2.1 -44.5 104.4 93'
--viewbox => svgViewbox( ST_Expand( ST_Extent(geom), 2))
) AS svg FROM shapes
;
END;
\$BODY\$;

ALTER FUNCTION postgisftw.big_chart(text)
OWNER TO postgres;
``````

The big chart will be centered in the browser with no height and width specification in the first line. In addition, we have added a title line with player name, team, and number of goals.

Cell size chart function
``````  -- FUNCTION: postgisftw.cell_chart(text)

-- DROP FUNCTION IF EXISTS postgisftw.cell_chart(text);

-- pg_featureserv looks for functions in the 'postgisftw' schema

-- TO GET OUR CHART TO DISPLAY REASONABLY IN Google Sheets we need three adjustments:
--  1) 'shrink' the geometries to 40% of their default size applying ST_SCALE
--  2) Use ST_TRANSLATE to slide the upper left corner of the rink to (0,0) -- as much as practical
--  3) Add a hard coded set of height and width values to the final SVG document

-- ** IMPORTANT  SRID 32613 is 'fake', just needed a planar projection to work with the arbitrary X/Y of the rink

CREATE OR REPLACE FUNCTION postgisftw.cell_chart(
rec_num_id text DEFAULT '4668'::text -- if no id, then Landeskog goal
)
RETURNS TABLE(svg text)
LANGUAGE 'plpgsql'
COST 100
STABLE STRICT PARALLEL UNSAFE
ROWS 1000

AS \$BODY\$
BEGIN
RETURN QUERY

with half_rink as (

select st_translate(     -- Scale and Translate for Google Sheets use case
st_scale(q1.geom,0.4,0.4)
,-2,-18
) as geom from
(
-- we only want to display the offensive end of the rink
select st_intersection(therink.geom,ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from therink
) q1
),
goals as (
select
st_translate(   -- Scale and Translate for Google Sheets use case
st_scale(q2.geom,0.4,0.4)
,-2,-18 )
as geom from
(
-- collect all of a player's goals into a single geometry
SELECT ST_SetSRID(st_collect(geom)::geometry, 32613) as geom
from
(	SELECT
st_intersection(ST_SetSRID(goals.plot_pt,32613),ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from postgisftw.goals
WHERE rec_num = rec_num_id::INTEGER -- uses the rec_num id fed into the function
)q1
)q2
),

shapes AS (

-- Rink SVG + styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#2E9AFE',
'stroke-width', 0.5::text )
)
svg FROM half_rink

UNION ALL

-- goals SVG + styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#F5A9A9',
'stroke-width', 0.5::text,
'fill','#FF0040'
)
)
svg FROM goals

)
--IMPORTANT we are hard-coding the extent of the document + adding height and width explicitly
--Google Sheets needs the height and width of the document specified (apparently)

SELECT svgDoc( array_agg( shapes.svg ),
'0 0 43 36" width="43mm" height="36mm '  --**WARNING -- hard coded values
) AS svg FROM shapes
;
END;
\$BODY\$;

ALTER FUNCTION postgisftw.cell_chart(text)
OWNER TO postgres;
``````

The small chart is scaled down, and the content is shifted to the origin of the SVG coordinate space ( 0,0 is in the upper left corner). In the small chart, you can see we have added width and height parameters in millimeters.

## pg_featureserv = JSON and SVG URLs from Postgres

pg_featureserv is a project that will run a separate lightweight Go server on top of your Postgres database to expose JSON, Geojson, and (newly) SVGs. pg_featureserv requires data to have a spatial reference id (SRID). You can just use a stand-in SRID and that is what is in the example data and functions.

pg_featureserv can be run as its own server and you can also run it from inside your database using the Crunchy Bridge Container Apps feature, based off of the podman project. To build a pg_featureserv container in a database, you’ll run something like this:

``````SELECT run_container('-dt -p 5433:5433/tcp -e DATABASE_URL="postgres://postgres:Rb3bZ1VZ7dZIUiFiy62J0OHVZybYROJjoDId@p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5432/postgres" -e PGFS_SERVER_HTTPPORT=5433 docker.io/pramsey/pg_featureserv:latest');
``````

pg_featureserv will then be running in a web browser at a link like `http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433`. To get data from pg_featureserv, you make requests with URLs for the data you want. Here’s a couple of samples:

JSON:

``````http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/collections/postgisftw.goals/items.json
``````

SVG:

``````http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/functions/postgisftw.cell_chart/items.svg
``````

You can further qualify URLs with player details and other queries in the url strings, like this:

``````http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/functions/postgisftw.cell_chart/items.svg?player_id=8471677
``````

Here’s your resulting SVG file displayed at a browser URL.

## Google sheets import of JSON and SVG

To load the JSON data into a Google Sheet, the fastest approach is to use Apps Scripts to write some JavaScript to process the JSON into an array of rows containing an array of columns. In a Google Sheet, go to Extensions > App Scripts. Create a new blank script, and replace the generated code with the following:

``````function ImportJSON(url) {
const response = UrlFetchApp.fetch(url)
const jsonString = response.getContentText()
const data = JSON.parse(jsonString).features.map(feature => ({
point_x: feature.geometry.coordinates[0],
point_y: feature.geometry.coordinates[1],
...feature.properties,
}))

const columns = [
'players',
'team_scored',
'period_num',
'this_event',
'this_event_code',
'playerid',
'game_id',
'rec_num',
]

const rows = data.map(item => columns.map(column => item[column]))

return [columns].concat(rows)
}
``````

Please note this is a specialized version of the function to only select specific columns and reorder them. To make it more general, replace the columns variable instantiation with something like the following:

``````const columns = Object.keys(data[0])
``````

Consider renaming the file to something like `ImportJSON.gs` then save the file as the last step. This will provide an `ImportJSON()` function we can use within the Google Sheet. In your Google Sheet, enter into a cell like A1 and use the something like this to combine bring in the JSON data with a limit or filters:

``````=ImportJSON("http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/collections/postgisftw.goals/items.json?limit=50")
``````

This will load the JSON data into the Google Sheet starting at A1.

To make a SVG ready for Google Sheet, use the Google Sheets `IMAGE` function with the pg_featureserv URL. Concatenating the URL lets us tie the JSON row to the right SVG data.

``````=IMAGE(CONCATENATE("http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/functions/postgisftw.cell_chart/items.svg?rec_num_id=", B1))
``````

Here’s a view of our final spreadsheet

## Open source for the win

Crunchy Data supports both pg_featureserv and pg_svg until I saw this talk, I had no idea they would work together. Martin obviously did because he wrote them both! The best part is that our fully managed Crunchy Bridge supported all of this, so it was super easy for me to set this up on a test sersver and tear it down when I'm finished testing.

Thanks to Jay Zawrotny for the javascript assistance. A huge thank you to Brian Timoney for letting me experiment with his hockey code and SVGs.

### Crunchy Data

#### JSON and SVG from PostGIS into Google Sheets

At PostGIS Day 2023, one of our speakers showed off a really cool demo for getting JSON and SVGs in and out of Postgres / PostGIS and into Google Sheets. Brian Timoney put together several open source projects in such a cool way that I just had to try it myself. If you want to see his demo video, it is on YouTube. With Brian’s blessing, I’m writing up some additional details with a few of the sample code bits for those of you that want to try this or something similar for your own projects.

So what do we have here? We have Postgres data that is coming into Google sheets in real time, tied to a custom SVG.

Before we dive in, an overview of things that I’ll cover to make this happen :

• Running pg_featureserv from inside a database
• Serving JSON data from Postgres with pg_featureserv
• Functions to create SVG imasges
• Serving SVG images from Postgres with pg_featureserv
• JSON & SVGs delivered to Google Sheets

## pg_svg functions

Brian wrote some special stuff for the goal and shot data in his examples to get it just right so if you want to play with a sample of that, here’s some data for you.

Sample data
``````-- Regular season goals of Colorado Avalanche for 2021-2022 season
--

-- information derived from public NHL API
-- shot locations normalized for offensive end of ice
-- hex_id was derived for a demo for the 2022 PostGIS Day demo

SET client_encoding = 'UTF8';

CREATE TABLE public.goals (
rec_num integer,
game_id text,
this_event text,
this_event_code text,
players text,
playerid text,
player_role text,
x_coord text,
y_coord text,
team_scored text,
period_num text,
period_type text,
plot_x numeric,
plot_y numeric,
plot_pt public.geometry(Point,32613),
hex_id numeric
);

ALTER TABLE public.goals OWNER TO postgres;

INSERT INTO public.goals VALUES (26, '2021020005', 'Goal', 'COL24', 'Jack Johnson', '8471677', 'Scorer', '-77.0', '1.0', 'Colorado Avalanche', '1', 'REGULAR', 77.0, -1.0, '01010000000000000000405340000000000000F0BF', 258);
INSERT INTO public.goals VALUES (5997, '2021020982', 'Goal', 'SJS204', 'Darren Helm', '8471794', 'Scorer', '-76.0', '5.0', 'Colorado Avalanche', '1', 'REGULAR', 76.0, -5.0, '0101000000000000000000534000000000000014C0', 257);
INSERT INTO public.goals VALUES (2432, '2021020415', 'Goal', 'COL45', 'Darren Helm', '8471794', 'Scorer', '-75.0', '-3.0', 'Colorado Avalanche', '1', 'REGULAR', 75.0, 3.0, '01010000000000000000C052400000000000000840', 258);
``````
Sample data for the rink (this is a DXF CAD file encoded as a geometry string)
``````  -- Rink map created from random DXF file found on internet
--
-- ** IMPORTANT: SRID of 32613 is fake!! It is just an arbitrary X/Y plane

SET client_encoding = 'UTF8';

CREATE TABLE public.therink ( id_num integer, geom
public.geometry(MultiLineString,32613) );

ALTER TABLE public.therink OWNER TO postgres;

INSERT INTO public.therink VALUES (1,

``````

The special visual functions you will see with the shots and the rink are done with some functions for Postgres in pg_svg. This isn’t really an extension, it is just a handy functions you can load into any Postgres database to help create SVGs. To load this in your database download it and run something like:

``````psql postgres://postgres:123155012sdfxxcwsdfweweff@p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5432/postgres < pg-svg-lib.sql
``````

Brian wrote two sample SVG shot chart functions for this project. One large and one small that could fit inside a spreadsheet if that’s your final data destination.

Big chart function
``````-- FUNCTION: postgisftw.big_chart(text)

-- pg_featureserv looks for functions in the 'postgisftw' schema

-- ** IMPORTANT  SRID 32613 is 'fake', just needed a planar projection to work with the arbitrary X/Y of the rink

CREATE OR REPLACE FUNCTION postgisftw.big_chart(
player_id text DEFAULT '8476455'::text) -- if no id, then Landeskog goals
RETURNS TABLE(svg text)
LANGUAGE 'plpgsql'
COST 100
STABLE STRICT PARALLEL UNSAFE
ROWS 1000

AS \$BODY\$
BEGIN
RETURN QUERY

-- we only want to display the offensive end of the rink
with half_rink as (
select st_intersection(therink.geom,ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from therink
),

goals as (
-- collect all of a player's goals into a single geometry
select 	ST_SetSRID(st_collect(geom)::geometry, 32613) as geom
from
(	SELECT
st_intersection(ST_SetSRID(goals.plot_pt,32613),ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from postgisftw.goals
WHERE playerid ILIKE player_id || '%'

)q1
),
-- get player name, team total number of goals for title display
playerinfo AS (

select UPPER(q1.players) as this_player,min(q1.team_scored) as this_team
, count(q1.playerid)::text as num_goals from
(
select a.players,a.playerid,a.team_scored
from postgisftw.goals a WHERE playerid ILIKE player_id || '%'
) q1
group by q1.players

),

-- make the SVG document

shapes AS (

--rink styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#2E9AFE',
'stroke-width', 0.5::text )
)
svg FROM half_rink

UNION ALL
-- goals styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#F5A9A9',
'stroke-width', 0.5::text,
'fill','#FF0040'
)
)
svg FROM goals

UNION ALL
-- player name, team, and total number of goals
-- create an arbitrary point underneath the rink and reasonably centered
SELECT NULL, svgText(ST_SetSRID(ST_MakePoint(50,-46),32613),this_player||'  ('||this_team||') -  '||num_goals||' goals',
style => svgStyle(  'fill', '#585858', 'text-anchor', 'middle', 'font', 'bold 3px sans-serif' ) )
svg
from playerinfo

)
-- create the final viewbox
-- Martin Davis uses a sensible default of expanding the collected geometry
-- since our rink geometry is static, I hard coded a viewbox
SELECT svgDoc( array_agg( shapes.svg ),
viewbox => '-2.1 -44.5 104.4 93'
--viewbox => svgViewbox( ST_Expand( ST_Extent(geom), 2))
) AS svg FROM shapes
;
END;
\$BODY\$;

ALTER FUNCTION postgisftw.big_chart(text)
OWNER TO postgres;
``````

The big chart will be centered in the browser with no height and width specification in the first line. In addition, we have added a title line with player name, team, and number of goals.

Cell size chart function
``````  -- FUNCTION: postgisftw.cell_chart(text)

-- DROP FUNCTION IF EXISTS postgisftw.cell_chart(text);

-- pg_featureserv looks for functions in the 'postgisftw' schema

-- TO GET OUR CHART TO DISPLAY REASONABLY IN Google Sheets we need three adjustments:
--  1) 'shrink' the geometries to 40% of their default size applying ST_SCALE
--  2) Use ST_TRANSLATE to slide the upper left corner of the rink to (0,0) -- as much as practical
--  3) Add a hard coded set of height and width values to the final SVG document

-- ** IMPORTANT  SRID 32613 is 'fake', just needed a planar projection to work with the arbitrary X/Y of the rink

CREATE OR REPLACE FUNCTION postgisftw.cell_chart(
rec_num_id text DEFAULT '4668'::text -- if no id, then Landeskog goal
)
RETURNS TABLE(svg text)
LANGUAGE 'plpgsql'
COST 100
STABLE STRICT PARALLEL UNSAFE
ROWS 1000

AS \$BODY\$
BEGIN
RETURN QUERY

with half_rink as (

select st_translate(     -- Scale and Translate for Google Sheets use case
st_scale(q1.geom,0.4,0.4)
,-2,-18
) as geom from
(
-- we only want to display the offensive end of the rink
select st_intersection(therink.geom,ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from therink
) q1
),
goals as (
select
st_translate(   -- Scale and Translate for Google Sheets use case
st_scale(q2.geom,0.4,0.4)
,-2,-18 )
as geom from
(
-- collect all of a player's goals into a single geometry
SELECT ST_SetSRID(st_collect(geom)::geometry, 32613) as geom
from
(	SELECT
st_intersection(ST_SetSRID(goals.plot_pt,32613),ST_SetSRID(ST_MakeBox2D(ST_Point(-0.1, 42.55), ST_Point(101, -42.55)), 32613))
as geom
from postgisftw.goals
WHERE rec_num = rec_num_id::INTEGER -- uses the rec_num id fed into the function
)q1
)q2
),

shapes AS (

-- Rink SVG + styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#2E9AFE',
'stroke-width', 0.5::text )
)
svg FROM half_rink

UNION ALL

-- goals SVG + styling
SELECT geom, svgShape( geom,
style => svgStyle(  'stroke', '#F5A9A9',
'stroke-width', 0.5::text,
'fill','#FF0040'
)
)
svg FROM goals

)
--IMPORTANT we are hard-coding the extent of the document + adding height and width explicitly
--Google Sheets needs the height and width of the document specified (apparently)

SELECT svgDoc( array_agg( shapes.svg ),
'0 0 43 36" width="43mm" height="36mm '  --**WARNING -- hard coded values
) AS svg FROM shapes
;
END;
\$BODY\$;

ALTER FUNCTION postgisftw.cell_chart(text)
OWNER TO postgres;
``````

The small chart is scaled down, and the content is shifted to the origin of the SVG coordinate space ( 0,0 is in the upper left corner). In the small chart, you can see we have added width and height parameters in millimeters.

## pg_featureserv = JSON and SVG URLs from Postgres

pg_featureserv is a project that will run a separate lightweight Go server on top of your Postgres database to expose JSON, Geojson, and (newly) SVGs. pg_featureserv requires data to have a spatial reference id (SRID). You can just use a stand-in SRID and that is what is in the example data and functions.

pg_featureserv can be run as its own server and you can also run it from inside your database using the Crunchy Bridge Container Apps feature, based off of the podman project. To build a pg_featureserv container in a database, you’ll run something like this:

``````SELECT run_container('-dt -p 5433:5433/tcp -e DATABASE_URL="postgres://postgres:Rb3bZ1VZ7dZIUiFiy62J0OHVZybYROJjoDId@p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5432/postgres" -e PGFS_SERVER_HTTPPORT=5433 docker.io/pramsey/pg_featureserv:latest');
``````

pg_featureserv will then be running in a web browser at a link like `http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433`. To get data from pg_featureserv, you make requests with URLs for the data you want. Here’s a couple of samples:

JSON:

``````http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/collections/postgisftw.goals/items.json
``````

SVG:

``````http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/functions/postgisftw.cell_chart/items.svg
``````

You can further qualify URLs with player details and other queries in the url strings, like this:

``````http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/functions/postgisftw.cell_chart/items.svg?player_id=8471677
``````

Here’s your resulting SVG file displayed at a browser URL.

## Google sheets import of JSON and SVG

To load the JSON data into a Google Sheet, the fastest approach is to use Apps Scripts to write some JavaScript to process the JSON into an array of rows containing an array of columns. In a Google Sheet, go to Extensions > App Scripts. Create a new blank script, and replace the generated code with the following:

``````function ImportJSON(url) {
const response = UrlFetchApp.fetch(url)
const jsonString = response.getContentText()
const data = JSON.parse(jsonString).features.map(feature => ({
point_x: feature.geometry.coordinates[0],
point_y: feature.geometry.coordinates[1],
...feature.properties,
}))

const columns = [
'players',
'team_scored',
'period_num',
'this_event',
'this_event_code',
'playerid',
'game_id',
'rec_num',
]

const rows = data.map(item => columns.map(column => item[column]))

return [columns].concat(rows)
}
``````

Please note this is a specialized version of the function to only select specific columns and reorder them. To make it more general, replace the columns variable instantiation with something like the following:

``````const columns = Object.keys(data[0])
``````

Consider renaming the file to something like `ImportJSON.gs` then save the file as the last step. This will provide an `ImportJSON()` function we can use within the Google Sheet. In your Google Sheet, enter into a cell like A1 and use the something like this to combine bring in the JSON data with a limit or filters:

``````=ImportJSON("http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/collections/postgisftw.goals/items.json?limit=50")
``````

This will load the JSON data into the Google Sheet starting at A1.

To make a SVG ready for Google Sheet, use the Google Sheets `IMAGE` function with the pg_featureserv URL. Concatenating the URL lets us tie the JSON row to the right SVG data.

``````=IMAGE(CONCATENATE("http://p.aqlunx3nqvebfh2bgffjj364ey.db.postgresbridge.com:5433/functions/postgisftw.cell_chart/items.svg?rec_num_id=", B1))
``````

Here’s a view of our final spreadsheet

## Open source for the win

Crunchy Data supports both pg_featureserv and pg_svg until I saw this talk, I had no idea they would work together. Martin obviously did because he wrote them both! The best part is that our fully managed Crunchy Bridge supported all of this, so it was super easy for me to set this up on a test sersver and tear it down when I'm finished testing.

Thanks to Jay Zawrotny for the javascript assistance. A huge thank you to Brian Timoney for letting me experiment with his hockey code and SVGs.

## January 07, 2024

### Paul Ramsey

#### PgConf.Dev @ Vancouver, May 28-31

This year, the global gathering of PostgreSQL developers has a new name, and a new location (but more-or-less the same dates) … pgcon.org is now pgconf.dev!

Some important points right up front:

• The call for papers is closing in one week! If you are planning to submit, now is the time!
• The hotel scene in Vancouver is competitive, so if you put off booking accomodations… don’t do that! Book a room right away.
• The venue capacity is 200. That’s it, so once we have 200 registrants, we are full for this year. Register now.

I first attended pgcon.org in 2011, when I was invited to keynote on the topic of PostGIS. Speaking in front of an audience of PostgreSQL luminaries was really intimidating, but also gratifying and empowering. Notwithstanding my imposter syndrome, all those super clever developers thought our little geospatial extension was… kind of clever.

I kept going to PgCon as regularly as I was able over the years, and was never disappointed. The annual gathering of the core developers of PostgreSQL necessarily includes content and insignts that you simply can not come across elsewhere, all compactly in one smallish conference, and the hallway track is amazing.

PostgreSQL may be a global development community, but the power of personal connection is not to be denied. Getting to meet and talk with core developers helped me understand where the project was going, and gave me the confidence to push ahead with my (very tiny) contributions.

This year, the event is in Vancouver! Still in Canada, but a little more directly connected to international air hubs than Ottawa was.

Also, this year I am honored to get a chance to serve on the program committee! We are looking for technical talks from across the PostgreSQL ecosystem, as well as about happenings in core. PostgreSQL is so much larger than just the core, and spreading the word about how you are building on PostgreSQL is important (and I am not just saying that as an extension author).

I hope to see you all there!

## December 19, 2023

### Paul Ramsey

#### Data Science is Getting Ducky

For a long time, a big constituency of users of PostGIS has been people with large data analytics problems that crush their desktop GIS systems. Or people who similarly find that their geospatial problems are too large to run in R. Or Python.

These are data scientists or adjacent people. And when they ran into those problems, the first course of action would be to move the data and parts of the workload to a “real database server”.

This all made sense to me.

But recently, something transformative happened – Crunchy Data upgraded my work laptop to a MacBook Pro.

Suddenly a GEOS compile that previously took 20 minutes, took 45 seconds.

I now have processing power on my local laptop that previously was only available on a server. The MacBook Pro may be a leading indicator of this amount of power, but the trend is clear.

What does that mean for default architectures and tooling?

Well, for data science, it means that a program like DuckDB goes from being a bit of a curiosity, to being the default tool for handling large data processing workloads.

What is DuckDB? According to the web site, it is “an in-process SQL OLAP database management system”. That doesn’t sound like a revolution in data science (it sounds really confusing).

But consider what DuckDB rolls together:

• A column-oriented processing engine that makes the most efficient possible use of the processors in modern computers. Parallelism to ensure all CPUs are made use of, and low-level optimizations to ensure each tick of those processors pushes as much data through the pipe as possible.
• Wide ranging support for different data formats, so that integration can take place on-the-fly without requiring translation or sometimes even data download steps.

Having those things together makes it a data science power tool, and removes a lot of the prior incentive that data scientists had to move their data into “real” databases.

When they run into the limits of in-memory analysis in R or Python, they will instead serialize their data to local disk and use DuckDB to slam through the joins and filters that were blowing out their RAM before.

They will also take advantage of DuckDB’s ability to stream remote data from data lake object stores.

What, stream multi-gigabyte JSON files? Well, yes that’s possible, but it’s not where the action is.

The CPU is not the only laptop component that has been getting ridiculously powerful over the past few years. The network pipe that connects that laptop to the internet has also been getting both wider and lower latency with every passing year.

As the propect of streaming data for analysis has come into view, the formats for remote data have also evolved. Instead of JSON, which is relatively fluffy, and hard to efficiently filter, the Parquet format is becoming a new standard for data lakes.

Parquet is a binary format, that organizes the data into blocks for efficient subsetting and processing. A DuckDB query to a properly organized Parquet time series file might easily pull only records for 2 of 20 columns, and 1 day of 365, reducing a multi-gigabyte download to a handful of megabytes.

The huge rise in available local computation, and network connectivity is going to spawn some new standard architectures.

Imagine a “two tier” architecture where tier one is an HTTP object store and tier two is a Javascript single page app? The COG Explorer has already been around for a few years, and it’s just such a two tier application.

(For fun, recognize that an architecture where the data are stored in an access-optimized format, and access is via primitive file-system requests, while all the smarts are in the client-side visualization software is… the old workstation GIS model. Everything old is new again.)

The technology is fresh, but the trendline is pretty clear. See Kyle Barrron’s talk about GeoParquet and DeckGL for a taste of where we are going.

Meanwhile, I expect that a lot of the growth in PostGIS / PostgreSQL we have seen in the data science field will level out for a while, as the convenience of DuckDB takes over a lot of workloads.

The limitations of Parquet (efficient remote access limited to a handful of filter variables being the primary one, as will cojoint spatial/non-spatial filter and joins) will still leave use cases that require a “real” database, but a lot of people who used to reach for PostGIS will be reaching for Duck, and that is going to change a lot of architectures, some for the better, and some for the worse.

### Crunchy Data

#### PostGIS Clustering with DBSCAN

A common problem in geospatial analysis is extracting areas of density from point fields. PostGIS has four window clustering functions that take in geometries and return cluster numbers (or NULL for unclustered inputs), which apply different algorithms to the problem of grouping the geometries in the input partitions.

The ST_ClusterDBSCAN function in PostGIS is a quick and easy way to extract clusters from point data. DBSCAN specifically works with density and is well suited for population or density type spatial data. To demonstrate `ST_ClusterDBSCAN` I'm going to work with the geographic names data, specifically the schools, and show how we can quickly create a U.S. population density map.

## Geographic Names Data

Let's explore clustering using geographic names data.

Create a table to hold the data. Note that the table is generating the points automatically from the longitude/latitude (EPSG:4326) and transforming into a planar projection for the USA (EPSG:5070).

``````CREATE TABLE geonames (
geonameid integer,
name text,
asciiname text,
alternatenames text,
latitude float8,
longitude float8,
fclass char,
fcode text,
country text,
cc2 text,
population bigint,
elevation integer,
dem text,
timezone text,
modification date,
geom geometry(point, 5070)
GENERATED ALWAYS AS
(ST_Transform(ST_Point(longitude, latitude, 4326),5070)) STORED
);
``````

Now load the table. Note the super fun use of `PROGRAM` to pull data directly from the web and feed a `COPY`.

``````COPY geonames
WITH (FORMAT CSV, DELIMITER E'\t', HEADER false);
``````

(This trick only works using the `postgres` superuser, since it involves calling a program and writing to system disk. If you do not have superuser access, download and unzip the `US.TXT` file by hand and load it using `COPY` from the file.)

Finally, add a spatial index to the `geom` column.

``````CREATE INDEX geonames_geom_x
ON geonames
USING GIST (geom);
``````

## Schools

There are 434 distinct feature codes in the `geonames` table. We will restrict our analysis to just the 205,848 schools, with an `fcode` of `SCH`.

``````SELECT Count(DISTINCT fcode) FROM geonames;
SELECT Count(fcode) FROM geonames WHERE fcode = 'SCH';
``````

Schools are an interesting feature to analyze because there's a nice strong correlation between the number of schools and the population. There's a lot of schools! But they are not uniformly distributed.

If we zoom into the midwest, the concentration of schools in populated places pops out. We can use PostGIS to turn this distribution difference into a data set of populated places!

## Clustering on Schools

The DBSCAN clustering algorithm is a "density based spatial clustering of applications with noise". The PostGIS ST_ClusterDBSCAN implementation is a window function that takes three parameters:

• The geometries to be analyzed for clusters.
• A 'eps' distance tolerance. Geometries must be within this distance to be added to a cluster.
• A 'minpoints' count. If a point is within the 'eps' distance of 'minpoints' cluster members, it is a "core member" of the cluster.

An input geometry is added to a cluster if it is either:

• A "core" geometry, that is within eps distance of at least minpoints input geometries (including itself); or
• A "border" geometry, that is within eps distance of a core geometry.

How does this play out in practice?

If we zoom further into Chicago, around the suburban/exurban transition, the schools are about 1000 meters apart, sometimes more sometimes less, transitioning out to 2000 meters and more in the exurbs.

For our clusters, we will use:

• A `eps` distance of 200m
• A `minpoints` of 5
• A partition on the state code (`admin1`) to cut down on the number of cluster numbers.
``````CREATE TABLE geonames_sch AS
SELECT ST_ClusterDBScan(geom, 2000, 5)
OVER (PARTITION BY admin1) AS cluster, *
FROM geonames
WHERE fcode = 'SCH';
``````

The result looks like this, with each cluster given a distinct color, and un-clustered schools rendered transparent.

The smaller clusters look a little arbitrary, but if we zoom in, we can see that even small population centers have been surfaced with this analytical technique.

Here is Kanakee, Illinois, neatly identified as a populated place by its cluster of schools.

## Clusters to Points

Now that we have clusters, getting a populated place point is as simple as using the ST_Centroid function.

``````CREATE TABLE geonames_popplaces AS
SELECT ST_Centroid(ST_Collect(geom))::geometry(Point, 5070) AS geom,
Count(*) AS school_count,
FROM geonames_sch
``````

We have completed the analysis, converting the density difference in school locations into a set of derived populated place points.

Now for the whole population cluster map!

## Quick recap

• Create a table `ST_ClusterDBScan`
• Set an `eps` for distance tolerance
• Set a `minpoints` to reduce density
• Partition on a different field to cut down on the number of cluster numbers.
• Create a final table using the `ST_Centroid`

### Crunchy Data

#### PostGIS Clustering with DBSCAN

A common problem in geospatial analysis is extracting areas of density from point fields. PostGIS has four window clustering functions that take in geometries and return cluster numbers (or NULL for unclustered inputs), which apply different algorithms to the problem of grouping the geometries in the input partitions.

The ST_ClusterDBSCAN function in PostGIS is a quick and easy way to extract clusters from point data. DBSCAN specifically works with density and is well suited for population or density type spatial data. To demonstrate `ST_ClusterDBSCAN` I'm going to work with the geographic names data, specifically the schools, and show how we can quickly create a U.S. population density map.

## Geographic Names Data

Let's explore clustering using geographic names data.

Create a table to hold the data. Note that the table is generating the points automatically from the longitude/latitude (EPSG:4326) and transforming into a planar projection for the USA (EPSG:5070).

``````CREATE TABLE geonames (
geonameid integer,
name text,
asciiname text,
alternatenames text,
latitude float8,
longitude float8,
fclass char,
fcode text,
country text,
cc2 text,
population bigint,
elevation integer,
dem text,
timezone text,
modification date,
geom geometry(point, 5070)
GENERATED ALWAYS AS
(ST_Transform(ST_Point(longitude, latitude, 4326),5070)) STORED
);
``````

Now load the table. Note the super fun use of `PROGRAM` to pull data directly from the web and feed a `COPY`.

``````COPY geonames
WITH (FORMAT CSV, DELIMITER E'\t', HEADER false);
``````

(This trick only works using the `postgres` superuser, since it involves calling a program and writing to system disk. If you do not have superuser access, download and unzip the `US.TXT` file by hand and load it using `COPY` from the file.)

Finally, add a spatial index to the `geom` column.

``````CREATE INDEX geonames_geom_x
ON geonames
USING GIST (geom);
``````

## Schools

There are 434 distinct feature codes in the `geonames` table. We will restrict our analysis to just the 205,848 schools, with an `fcode` of `SCH`.

``````SELECT Count(DISTINCT fcode) FROM geonames;
SELECT Count(fcode) FROM geonames WHERE fcode = 'SCH';
``````

Schools are an interesting feature to analyze because there's a nice strong correlation between the number of schools and the population. There's a lot of schools! But they are not uniformly distributed.

If we zoom into the midwest, the concentration of schools in populated places pops out. We can use PostGIS to turn this distribution difference into a data set of populated places!

## Clustering on Schools

The DBSCAN clustering algorithm is a "density based spatial clustering of applications with noise". The PostGIS ST_ClusterDBSCAN implementation is a window function that takes three parameters:

• The geometries to be analyzed for clusters.
• A 'eps' distance tolerance. Geometries must be within this distance to be added to a cluster.
• A 'minpoints' count. If a point is within the 'eps' distance of 'minpoints' cluster members, it is a "core member" of the cluster.

An input geometry is added to a cluster if it is either:

• A "core" geometry, that is within eps distance of at least minpoints input geometries (including itself); or
• A "border" geometry, that is within eps distance of a core geometry.

How does this play out in practice?

If we zoom further into Chicago, around the suburban/exurban transition, the schools are about 1000 meters apart, sometimes more sometimes less, transitioning out to 2000 meters and more in the exurbs.

For our clusters, we will use:

• A `eps` distance of 2000m
• A `minpoints` of 5
• A partition on the state code (`admin1`) to cut down on the number of cluster numbers.
``````CREATE TABLE geonames_sch AS
SELECT ST_ClusterDBScan(geom, 2000, 5)
OVER (PARTITION BY admin1) AS cluster, *
FROM geonames
WHERE fcode = 'SCH';
``````

The result looks like this, with each cluster given a distinct color, and un-clustered schools rendered transparent.

The smaller clusters look a little arbitrary, but if we zoom in, we can see that even small population centers have been surfaced with this analytical technique.

Here is Kanakee, Illinois, neatly identified as a populated place by its cluster of schools.

## Clusters to Points

Now that we have clusters, getting a populated place point is as simple as using the ST_Centroid function.

``````CREATE TABLE geonames_popplaces AS
SELECT ST_Centroid(ST_Collect(geom))::geometry(Point, 5070) AS geom,
Count(*) AS school_count,
FROM geonames_sch
``````

We have completed the analysis, converting the density difference in school locations into a set of derived populated place points.

Now for the whole population cluster map!

## Quick recap

• Create a table `ST_ClusterDBScan`
• Set an `eps` for distance tolerance
• Set a `minpoints` to reduce density
• Partition on a different field to cut down on the number of cluster numbers.
• Create a final table using the `ST_Centroid`

## December 12, 2023

### Paul Ramsey

#### Keynote @ FOSS4G NA 2023

Preparing the keynote for FOSS4G North America this year felt particularly difficult. I certainly sweated over it.

• Audience was a problem. I wanted to talk about my usual thing, business models and economics, but the audience was going to be a mash of people new to the topic and people who has seen my spiel multiple times.
• Length was a problem. Out of an excess of faith in my abilities, the organizers gave me a full hour long slot! That is a very long time to keep people’s attention and try to provide something interesting.

The way it all ended up was:

• Cadging some older content from keynotes about business models, to bring new folks up to speed.
• Mixing in some only slightly older content about cloud models.
• Adding in some new thoughts about the way everyone can work together to make open source more sustainable (or at least less extractive) over the long term.

Here’s this year’s iteration.

The production of this kind of content is involved. The goal is to remain interesting over a relatively long period of time.

I have become increasingly opinionated about how to do that.

• No freestyling. Blathering over bullet points is unfair to your audience. The aggregate time of an audience of 400 is very large. 5 minutes of your “um” and “ah” translates into 33 hours of dead audience time.
• Get right to it. No mini-resume, no talking about your employer (unless you are really sneaky about it, like me 😉), this is about delivering ideas and facts that are relevant to the audience. Your introducer can handle your bona fides.
• Have good content. The hardest part! (?) Do you have something thematic you can bookend the start and end with? Are there some interesting facts that much of the audience does not know yet? Are there some unappreciated implications? This is, presumably, why you were asked to keynote, so hopefully not too, too hard. This is the part that I worry over the most, because I really have no faith that what I have to say is actually going be interesting to an audience, no matter how much I gussy it up.
• Work from a text. The way to avoid blather is to know exactly what you are going to say. At 140 words-per-minute speaking pace, a 55 minute talk is 7700 words, which coincidentally (not) is exactly how long my keynote text is.
• Write a speech, not an article. You will have to say all those words! Avoid complicated sentence constructions. Keep sentences short. Take advantage of parallel constructions to make a point, drive a narrative, force a conclusion. (see?) Repeat yourself. Repeat yourself.
• Perform, don’t read. Practice reading out loud. Get used to leaving longer gaps and get comfortable with silence. Practice modulating your voice. Louder, softer. Faster, slower. Drop. The. Hammer. Watch a gifted speaker like Barack Obama deliver a text. He isn’t ad libbing, he’s performing a prepared text. See what he does to make that sound spontaneous and interesting.
• Visuals as complements, not copies. Your slides should complement and amplify your content, not recapitulate it. In the limit, you could do all-text slides, which just give the three-word summary of your current main point. (This classic Lessig talk is my favourite example.)
• Visuals as extra channel. Keep changing up the visual. Use the slide notes space to get a feel for how long each slide should be up. (Hint, about 50 words on average.) Keeping slide duration low also helps in terms of using the per-slide speaker notes as low-end teleprompter (increase notes font size! reduce slide preview size!) from which you deliver your performance.

I originally started scripting talks because it allowed me to smooth out the quality of my talks. With a script, it wasn’t a crapshoot whether I had a good ad lib delivery or a bad one, I had a nice consistent level. From there, leveraging up to take advantage of the format to increase the talk quality was a natural step. Speakers like Lessig and Damian Conway remain my guide posts.

If you liked the keynote video and want to use the materials, the slides are available here under CC BY.

## December 07, 2023

### Boston GIS (Regina Obe, Leo Hsu)

#### PostGIS Day 2023 Summary

PostGIS Day 2023 videos came out recently. PostGIS Day conference is always my favorite conference of the year because you get to see what people are doing all over the world, and it always has many many new tricks for using PostgreSQL and PostGIS family of extensions you had never thought of. Most importantly it's virtual, which makes it much easier for people to fit in their schedules than an on site conference. We really need more virtual conferences in the PostgreSQL community. Many many thanks to Crunchy Data for putting this together again, in particular to Elizabeth Christensen who did the hard behind the scenes work of corraling all the presenters and stepping in to give a talk herself, and my PostGIS partner in development Paul Ramsey who did the MC'ing probably with very little sleep, but still managed to be very energetic. Check out Elizabeth's summary of the event. Many of her highlights would have been mine too, so I'm going to skip those.

Continue reading "PostGIS Day 2023 Summary"

## December 06, 2023

### Crunchy Data

#### PostGIS Day 2023 Summary

We hosted our annual PostGIS day a couple weeks ago with some great talks on a big variety of topics within open-source GIS. Here is a summary of the themes I saw take shape across the day’s events that will point you towards the recordings, depending on your interests. A full playlist of PostGIS Day 2023 is available on our YouTube channel.

## The Ideal GIS Stack

If you’ve spent time with developers this year you know that folks love to tell you the details and reasoning behind their tech stack and the GIS community is no different. PostGIS is really the engine behind the modern day open-source GIS installations and several presenters came to talk about their GIS architecture and preferred toolchains, backed by PostGIS.

The PRAM stack

Paul asked Rhys if the PRAM stack name is flattery or if he’s being trolled. I think you’ll see that Rhys’ talk is almost all flattery. He digs into using PostGIS with the ogr_fdw, pgRouting, PostgREST, Sqitch, and pgTAP for projects in the utilities industry. Rhys’ talk had some of the best screengrabs including this gem.

Modern GIS Analytics

Matt Forrest from Carto also had a great talk on the analytics workflows and his take on a modern GIS analytics stack. He had details on using Geoparquet, DuckDB, dbt, and H3. What connects everything, in Matt’s opinion is, SQL. He had some great slides, including this one (bonus points for putting PostGIS at the center of everything):

PostgREST & making PostGIS your modern REST API

PostgREST turns your Postgres database, and PostGIS, into a REST API and it is pretty cool. Krishna has a great technical overview of PostgREST and how to get started with this, including some of the tricks in working with the Swagger API Platform.

## Movement, Mobility, and Emergency Services

We had several really good talks from people working in the field to solve issues with getting people and things where they need to go. We had two speakers in the emergency services sector.

Laure-Hélène Bruneton from CamptoCamp talked about her work at NexSIS emergency operations management in France with her talk “Custom Road Network Contraction for Routing”. She digs into working with a large routing dataset and some tips for reducing size and making it more performant.

Randal Hale came to talk about his work with 911 in rural Tennessee, comparing using just the Geopackage files for 911 to using PostGIS. Without even storing data in a database, Randal is able to work with data, QGIS, and SQL all through a Geopackage file.

Vicky Vergara, the primary developer behind pgRouting came to talk about a special pgRouting project she has been developing for a UN initiative on how close people live to hospitals. She demonstrated how this data is accessible with open source tools and open data and you can see how important something like this would be for a developing nation.

Ford Motor Company’s research and development team presented on using PostGIS and pgRouting in their BlueCruise Hands-Free Highway Driving technology. I love that PostGIS is literally behind the wheel! Brendan Farrell came to talk about “Mapping Where the Data is Not”. This is when you’re dealing with missing pieces of data, denoting that, and using complementary geometry.

Initially released in early 2022, MobilityDB, has been getting more and more attention and hands-on love in 2023 and we ended up with two talks digging into some details. MobilityDB is an extension that is built on top of PostGIS and Postgres and expands the capabilities even more for moving geospatial objects. The main project leader, Prof Esteban Zimanyi, gave a great overview of how this fits in. Following Zimanyi’s talk was Wendell Turner with “Air Traffic Analysis with PostGIS and MobilityDB”. He used airplane, airport, and weather data to do demo research on everyone’s favorite topic, flight delays.

## Leveling up your SQL & PostGIS skills

We’re always blessed to have some of the best and brightest come to share their expertise. One great thing about PostGIS day is that it’s kind of a mix of hearing stories, learning about projects and tools with demos of technical skills outside of your day to day life.

I presented a talk on being a “Spatial DBA in a Pinch” as I’ve found lots of folks end up taking care of a database even though they didn’t really set out to do that. This is a basic talk, but if you’re new to PostGIS, there are handy tips about creating roles, looking at queries, and crating basic indexes. Paul added some really great tips on PostGIS performance too, which expand on some of the basics I presented.

Regina Obe always brings the coolest PostGIS examples and this year she did not disappoint her fans. She showed off a bunch of Star Wars graphics with her talk on “PostGIS surprise extensions”. If you want to show your kid some PostgreSQL over the holidays, Regina has all the code ready to go so you can run this yourself.

We had a couple other really great technical deep dives. Benjamin Trigona-Harany talked about processing airplane flight data in “Trajectory Analysis Using PostGIS” (with some tips on getting free airline and flight path data). Crunchy Data’s own Martin Davis discussed some new PostGIS features for handling “Simple polygonal coverages”, including validation, union and simplification (which he calls the “killer app” for coverages).

## The role as GIS professionals

We got some nice breaks from technical talks to spend time thinking about our role as GIS professionals in this big wide world and what all of it means. Bonnie McClain came to talk to us on how “We are part of the infrastructure, not above it.” Bonnie has been doing amazing things in the GIS world by telling stories, looking at urban data through the lens of both GIS and human geography. Her talk digs into some of the things she uncovered on a recent project where she uses GIS data to build a Global and Healthy Sustainable Cities Indicator.

Brian Timoney talked about “Refactoring the Way We Talk About SQL”. This talk reviews open source software, value, and how we talk about our roles. He also gives an incredible demo of getting data straight from PostGIS through the pg_svg SVG extension and pg_featureserv and then into a spreadsheet.

Brian is a uniquely talented speaker with a love for this line of work that’s rarely communicated into words. His recent talk from FOSS4G NA on “You Can’t Get There From Here Alone” is also worth a mention here; it’s excellent (I had to hold back tears when I saw this in person).

## Thank you!

Thanks to everyone who participated this year by speaking, coming to the event, chatting with us, or waiting until the videos are up on YouTube to catch on this year’s PostGIS Day. We had event attendees from more than 54 countries! We’ll be posting a call for papers next September so keep an eye out for that.

### Crunchy Data

#### PostGIS Day 2023 Summary

We hosted our annual PostGIS day a couple weeks ago with some great talks on a big variety of topics within open-source GIS. Here is a summary of the themes I saw take shape across the day’s events that will point you towards the recordings, depending on your interests. A full playlist of PostGIS Day 2023 is available on our YouTube channel.

## The Ideal GIS Stack

If you’ve spent time with developers this year you know that folks love to tell you the details and reasoning behind their tech stack and the GIS community is no different. PostGIS is really the engine behind the modern day open-source GIS installations and several presenters came to talk about their GIS architecture and preferred toolchains, backed by PostGIS.

The PRAM stack

Paul asked Rhys if the PRAM stack name is flattery or if he’s being trolled. I think you’ll see that Rhys’ talk is almost all flattery. He digs into using PostGIS with the ogr_fdw, pgRouting, PostgREST, Sqitch, and pgTAP for projects in the utilities industry. Rhys’ talk had some of the best screengrabs including this gem.

Modern GIS Analytics

Matt Forrest from Carto also had a great talk on the analytics workflows and his take on a modern GIS analytics stack. He had details on using Geoparquet, DuckDB, dbt, and H3. What connects everything, in Matt’s opinion is, SQL. He had some great slides, including this one (bonus points for putting PostGIS at the center of everything):

PostgREST & making PostGIS your modern REST API

PostgREST turns your Postgres database, and PostGIS, into a REST API and it is pretty cool. Krishna has a great technical overview of PostgREST and how to get started with this, including some of the tricks in working with the Swagger API Platform.

## Movement, Mobility, and Emergency Services

We had several really good talks from people working in the field to solve issues with getting people and things where they need to go. We had two speakers in the emergency services sector.

Laure-Hélène Bruneton from CamptoCamp talked about her work at NexSIS emergency operations management in France with her talk “Custom Road Network Contraction for Routing”. She digs into working with a large routing dataset and some tips for reducing size and making it more performant.

Randal Hale came to talk about his work with 911 in rural Tennessee, comparing using just the Geopackage files for 911 to using PostGIS. Without even storing data in a database, Randal is able to work with data, QGIS, and SQL all through a Geopackage file.

Vicky Vergara, the primary developer behind pgRouting came to talk about a special pgRouting project she has been developing for a UN initiative on how close people live to hospitals. She demonstrated how this data is accessible with open source tools and open data and you can see how important something like this would be for a developing nation.

Ford Motor Company’s research and development team presented on using PostGIS and pgRouting in their BlueCruise Hands-Free Highway Driving technology. I love that PostGIS is literally behind the wheel! Brendan Farrell came to talk about “Mapping Where the Data is Not”. This is when you’re dealing with missing pieces of data, denoting that, and using complementary geometry.

Initially released in early 2022, MobilityDB, has been getting more and more attention and hands-on love in 2023 and we ended up with two talks digging into some details. MobilityDB is an extension that is built on top of PostGIS and Postgres and expands the capabilities even more for moving geospatial objects. The main project leader, Prof Esteban Zimanyi, gave a great overview of how this fits in. Following Zimanyi’s talk was Wendell Turner with “Air Traffic Analysis with PostGIS and MobilityDB”. He used airplane, airport, and weather data to do demo research on everyone’s favorite topic, flight delays.

## Leveling up your SQL & PostGIS skills

We’re always blessed to have some of the best and brightest come to share their expertise. One great thing about PostGIS day is that it’s kind of a mix of hearing stories, learning about projects and tools with demos of technical skills outside of your day to day life.

I presented a talk on being a “Spatial DBA in a Pinch” as I’ve found lots of folks end up taking care of a database even though they didn’t really set out to do that. This is a basic talk, but if you’re new to PostGIS, there are handy tips about creating roles, looking at queries, and crating basic indexes. Paul added some really great tips on PostGIS performance too, which expand on some of the basics I presented.

Regina Obe always brings the coolest PostGIS examples and this year she did not disappoint her fans. She showed off a bunch of Star Wars graphics with her talk on “PostGIS surprise extensions”. If you want to show your kid some PostgreSQL over the holidays, Regina has all the code ready to go so you can run this yourself.

We had a couple other really great technical deep dives. Benjamin Trigona-Harany talked about processing airplane flight data in “Trajectory Analysis Using PostGIS” (with some tips on getting free airline and flight path data). Crunchy Data’s own Martin Davis discussed some new PostGIS features for handling “Simple polygonal coverages”, including validation, union and simplification (which he calls the “killer app” for coverages).

## The role as GIS professionals

We got some nice breaks from technical talks to spend time thinking about our role as GIS professionals in this big wide world and what all of it means. Bonnie McClain came to talk to us on how “We are part of the infrastructure, not above it.” Bonnie has been doing amazing things in the GIS world by telling stories, looking at urban data through the lens of both GIS and human geography. Her talk digs into some of the things she uncovered on a recent project where she uses GIS data to build a Global and Healthy Sustainable Cities Indicator.

Brian Timoney talked about “Refactoring the Way We Talk About SQL”. This talk reviews open source software, value, and how we talk about our roles. He also gives an incredible demo of getting data straight from PostGIS through the pg_svg SVG extension and pg_featureserv and then into a spreadsheet.

Brian is a uniquely talented speaker with a love for this line of work that’s rarely communicated into words. His recent talk from FOSS4G NA on “You Can’t Get There From Here Alone” is also worth a mention here; it’s excellent (I had to hold back tears when I saw this in person).

## Thank you!

Thanks to everyone who participated this year by speaking, coming to the event, chatting with us, or waiting until the videos are up on YouTube to catch on this year’s PostGIS Day. We had event attendees from more than 54 countries! We’ll be posting a call for papers next September so keep an eye out for that.

## November 20, 2023

### PostGIS Development

#### PostGIS Patch Releases

The PostGIS development team is pleased to provide bug fix and performance enhancements 3.4.1, 3.3.5, 3.2.6, 3.1.10, 3.0.10 for the 3.4, 3.3, 3.2, 3.1, 3.0 stable branches.

## September 11, 2023

### Crunchy Data

#### Random Geometry Generation with PostGIS

A user on the postgis-users had an interesting question today: how to generate a geometry column in PostGIS with random points, linestrings, or polygons?

Random data is important for validating processing chains, analyses and reports. The best way to test a process is to feed it inputs!

## Random Points

Random points is pretty easy -- define an area of interest and then use the PostgreSQL `random()` function to create the X and Y values in that area.

``````CREATE TABLE random_points AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
)
SELECT ST_Point(width  * (random() - 0.5) + origin_x,
height * (random() - 0.5) + origin_y,
4326)::Geometry(Point, 4326) AS geom,
id
FROM bounds,
generate_series(0, 100) AS id
``````

Filling a target shape with random points is a common use case, and there's a special function just for that, `ST_GeneratePoints()`. Here we generate points inside a circle created with `ST_Buffer()`.

``````CREATE TABLE random_points AS
SELECT ST_GeneratePoints(
ST_Buffer(
ST_Point(0, 0, 4326),
50),
100) AS geom
``````

If you have PostgreSQL 16, you can use the brand new `random_normal()` function to generate coordinates with a central tendency.

``````CREATE TABLE random_normal_points AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
)
SELECT ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326)::Geometry(Point, 4326) AS geom,
id
FROM bounds,
generate_series(0, 100) AS id
``````

For PostgreSQL versions before 16, here is a user-defined version of `random_normal()`.
``````CREATE OR REPLACE FUNCTION random_normal(
mean double precision DEFAULT 0.0,
stddev double precision DEFAULT 1.0)
RETURNS double precision AS
\$\$
DECLARE
u1 double precision;
u2 double precision;
z0 double precision;
z1 double precision;
BEGIN
u1 := random();
u2 := random();

z0 := sqrt(-2.0 * ln(u1)) * cos(2.0 * pi() * u2);
z1 := sqrt(-2.0 * ln(u1)) * sin(2.0 * pi() * u2);

RETURN mean + (stddev * z0);
END;
\$\$ LANGUAGE plpgsql;
``````

## Random Linestrings

Linestrings are a little harder, because they involve more points, and aesthetically we like to avoid self-crossings of lines.

Two-point linestrings are pretty easy to generate with `ST_MakeLine()` -- just generate twice as many random points, and use them as the start and end points of the linestrings.

``````CREATE TABLE random_2point_lines AS
WITH bounds AS (
SELECT 0 AS origin_x, 80 AS width,
0 AS origin_y, 80 AS height
)
SELECT ST_MakeLine(
ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326),
ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326))::Geometry(LineString, 4326) AS geom,
id
FROM bounds,
generate_series(0, 100) AS id
``````

Multi-point random linestrings are harder, at least while avoiding self-intersections, and there are a lot of potential approaches. While a recursive CTE could probably do it, an imperative approach using PL/PgSQL is more readable.

The `generate_random_linestring()` function starts with an empty linestring, and then adds on new segments one at a time, changing the direction of the line with each new segment.

Here is the full `generate_random_linestring()` definition.
``````CREATE OR REPLACE FUNCTION generate_random_linestring(
start_point geometry(Point))
RETURNS geometry(LineString, 4326) AS
\$\$
DECLARE
num_segments integer := 10; -- Number of segments in the linestring
deviation_max float := radians(45); -- Maximum deviation
random_point geometry(Point);
deviation float;
direction float := 2 * pi() * random();
segment_length float := 5; -- Length of each segment (adjust as needed)
i integer;
result geometry(LineString) := 'SRID=4326;LINESTRING EMPTY';
BEGIN
FOR i IN 1..num_segments LOOP
-- Generate a random angle within the specified deviation
deviation := 2 * deviation_max * random() - deviation_max;
direction := direction + deviation;

-- Calculate the coordinates of the next point
random_point := ST_Point(
ST_X(start_point) + cos(direction) * segment_length,
ST_Y(start_point) + sin(direction) * segment_length,
ST_SRID(start_point)
);

-- Add the point to the linestring

-- Update the start point for the next segment
start_point := random_point;

END LOOP;

RETURN result;
END;
\$\$
LANGUAGE plpgsql;
``````

We can use the `generate_random_linestring()` function now to turn random start points (created in the usual way) into fully random squiggly lines!

``````CREATE TABLE random_lines AS
WITH bounds AS (
SELECT 0 AS origin_x, 80 AS width,
0 AS origin_y, 80 AS height
)
SELECT id,
generate_random_linestring(
ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326))::Geometry(LineString, 4326) AS geom
FROM bounds,
generate_series(1, 100) AS id;
``````

## Random Polygons

At the simplest level, a set of random boxes is a set of random polygons, but that's pretty boring, and easy to generate using `ST_MakeEnvelope()`.

``````CREATE TABLE random_boxes AS
WITH bounds AS (
SELECT 0 AS origin_x, 80 AS width,
0 AS origin_y, 80 AS height
)
SELECT ST_MakeEnvelope(
random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
random_normal(origin_x, width/4),
random_normal(origin_y, height/4)
)::Geometry(Polygon, 4326) AS geom,
id
FROM bounds,
generate_series(0, 20) AS id

``````

But more interesting polygons have curvy and convex shapes, how can we generate those?

### Random Polygons with Concave Hull

One way is to extract a polygon from a set of random points, using `ST_ConcaveHull()`, and then applying an "erode and dilate" effect to make the curves more pleasantly round.

We start with a random center point for each polygon, and create a circle with `ST_Buffer()`.

Then use `ST_GeneratePoints()` to fill the circle with some random points -- not too many, so we get a nice jagged result.

Then use `ST_ConcaveHull()` to trace a "boundary" around those points.

Then apply a negative buffer, to erode the shape.

And finally a positive buffer to dilate it back out again.

Generating multiple hulls involves stringing together all the above operations with CTEs or subqueries.

Here is the full query to generate multiple polygons with the concave hull method.
``````CREATE TABLE random_hulls AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
),
polypts AS (
SELECT ST_Point(random_normal(origin_x, width/2),
random_normal(origin_y, width/2),
4326)::Geometry(Point, 4326) AS geom,
polyid
FROM bounds,
generate_series(1,10) AS polyid
),
pts AS (
SELECT ST_GeneratePoints(ST_Buffer(geom, width/5), 20) AS geom,
polyid
FROM bounds,
polypts
)
SELECT ST_Multi(ST_Buffer(
ST_Buffer(
ST_ConcaveHull(geom, 0.3),
-2.0),
3.0))::Geometry(MultiPolygon, 4326) AS geom,
polyid
FROM pts;
``````

### Random Polygons with Voronoi Polygons

Another approach is to again start with random points, but use the Voronoi diagram as the basis of the polygon.

Generate random points in the circle.

Use the `ST_VoronoiPolygons()` function to generate polygons that subdivide the space using the random points as seeds.

Filter just the polygons that are fully contained in the originating circle.

And then use `ST_Union()` to merge those polygons into a single output shape.

Generating multiple hulls again involves stringing together the above operations with CTEs or subqueries.

Here is the full query to generate multiple polygons with the Voronoi method.
``````CREATE TABLE random_delaunay_hulls AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
),
polypts AS (
SELECT ST_Point(random_normal(origin_x, width/2),
random_normal(origin_y, width/2),
4326)::Geometry(Point, 4326) AS geom,
polyid
FROM bounds,
generate_series(1,20) AS polyid
),
voronois AS (
SELECT ST_VoronoiPolygons(
ST_GeneratePoints(ST_Buffer(geom, width/5), 10)
) AS geom,
ST_Buffer(geom, width/5) AS geom_clip,
polyid
FROM bounds,
polypts
),
cells AS (
SELECT (ST_Dump(geom)).geom, polyid, geom_clip
FROM voronois
)
SELECT ST_Union(geom)::Geometry(Polygon, 4326) AS geom, polyid
FROM cells
WHERE ST_Contains(geom_clip, geom)
GROUP BY polyid;
``````

### Crunchy Data

#### Random Geometry Generation with PostGIS

A user on the postgis-users had an interesting question today: how to generate a geometry column in PostGIS with random points, linestrings, or polygons?

Random data is important for validating processing chains, analyses and reports. The best way to test a process is to feed it inputs!

## Random Points

Random points is pretty easy -- define an area of interest and then use the PostgreSQL `random()` function to create the X and Y values in that area.

``````CREATE TABLE random_points AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
)
SELECT ST_Point(width  * (random() - 0.5) + origin_x,
height * (random() - 0.5) + origin_y,
4326)::Geometry(Point, 4326) AS geom,
id
FROM bounds,
generate_series(0, 100) AS id
``````

Filling a target shape with random points is a common use case, and there's a special function just for that, `ST_GeneratePoints()`. Here we generate points inside a circle created with `ST_Buffer()`.

``````CREATE TABLE random_points AS
SELECT ST_GeneratePoints(
ST_Buffer(
ST_Point(0, 0, 4326),
50),
100) AS geom
``````

If you have PostgreSQL 16, you can use the brand new `random_normal()` function to generate coordinates with a central tendency.

``````CREATE TABLE random_normal_points AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
)
SELECT ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326)::Geometry(Point, 4326) AS geom,
id
FROM bounds,
generate_series(0, 100) AS id
``````

For PostgreSQL versions before 16, here is a user-defined version of `random_normal()`.
``````CREATE OR REPLACE FUNCTION random_normal(
mean double precision DEFAULT 0.0,
stddev double precision DEFAULT 1.0)
RETURNS double precision AS
\$\$
DECLARE
u1 double precision;
u2 double precision;
z0 double precision;
z1 double precision;
BEGIN
u1 := random();
u2 := random();

z0 := sqrt(-2.0 * ln(u1)) * cos(2.0 * pi() * u2);
z1 := sqrt(-2.0 * ln(u1)) * sin(2.0 * pi() * u2);

RETURN mean + (stddev * z0);
END;
\$\$ LANGUAGE plpgsql;
``````

## Random Linestrings

Linestrings are a little harder, because they involve more points, and aesthetically we like to avoid self-crossings of lines.

Two-point linestrings are pretty easy to generate with `ST_MakeLine()` -- just generate twice as many random points, and use them as the start and end points of the linestrings.

``````CREATE TABLE random_2point_lines AS
WITH bounds AS (
SELECT 0 AS origin_x, 80 AS width,
0 AS origin_y, 80 AS height
)
SELECT ST_MakeLine(
ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326),
ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326))::Geometry(LineString, 4326) AS geom,
id
FROM bounds,
generate_series(0, 100) AS id
``````

Multi-point random linestrings are harder, at least while avoiding self-intersections, and there are a lot of potential approaches. While a recursive CTE could probably do it, an imperative approach using PL/PgSQL is more readable.

The `generate_random_linestring()` function starts with an empty linestring, and then adds on new segments one at a time, changing the direction of the line with each new segment.

Here is the full `generate_random_linestring()` definition.
``````CREATE OR REPLACE FUNCTION generate_random_linestring(
start_point geometry(Point))
RETURNS geometry(LineString, 4326) AS
\$\$
DECLARE
num_segments integer := 10; -- Number of segments in the linestring
deviation_max float := radians(45); -- Maximum deviation
random_point geometry(Point);
deviation float;
direction float := 2 * pi() * random();
segment_length float := 5; -- Length of each segment (adjust as needed)
i integer;
result geometry(LineString) := 'SRID=4326;LINESTRING EMPTY';
BEGIN
FOR i IN 1..num_segments LOOP
-- Generate a random angle within the specified deviation
deviation := 2 * deviation_max * random() - deviation_max;
direction := direction + deviation;

-- Calculate the coordinates of the next point
random_point := ST_Point(
ST_X(start_point) + cos(direction) * segment_length,
ST_Y(start_point) + sin(direction) * segment_length,
ST_SRID(start_point)
);

-- Add the point to the linestring

-- Update the start point for the next segment
start_point := random_point;

END LOOP;

RETURN result;
END;
\$\$
LANGUAGE plpgsql;
``````

We can use the `generate_random_linestring()` function now to turn random start points (created in the usual way) into fully random squiggly lines!

``````CREATE TABLE random_lines AS
WITH bounds AS (
SELECT 0 AS origin_x, 80 AS width,
0 AS origin_y, 80 AS height
)
SELECT id,
generate_random_linestring(
ST_Point(random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
4326))::Geometry(LineString, 4326) AS geom
FROM bounds,
generate_series(1, 100) AS id;
``````

## Random Polygons

At the simplest level, a set of random boxes is a set of random polygons, but that's pretty boring, and easy to generate using `ST_MakeEnvelope()`.

``````CREATE TABLE random_boxes AS
WITH bounds AS (
SELECT 0 AS origin_x, 80 AS width,
0 AS origin_y, 80 AS height
)
SELECT ST_MakeEnvelope(
random_normal(origin_x, width/4),
random_normal(origin_y, height/4),
random_normal(origin_x, width/4),
random_normal(origin_y, height/4)
)::Geometry(Polygon, 4326) AS geom,
id
FROM bounds,
generate_series(0, 20) AS id

``````

But more interesting polygons have curvy and convex shapes, how can we generate those?

### Random Polygons with Concave Hull

One way is to extract a polygon from a set of random points, using `ST_ConcaveHull()`, and then applying an "erode and dilate" effect to make the curves more pleasantly round.

We start with a random center point for each polygon, and create a circle with `ST_Buffer()`.

Then use `ST_GeneratePoints()` to fill the circle with some random points -- not too many, so we get a nice jagged result.

Then use `ST_ConcaveHull()` to trace a "boundary" around those points.

Then apply a negative buffer, to erode the shape.

And finally a positive buffer to dilate it back out again.

Generating multiple hulls involves stringing together all the above operations with CTEs or subqueries.

Here is the full query to generate multiple polygons with the concave hull method.
``````CREATE TABLE random_hulls AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
),
polypts AS (
SELECT ST_Point(random_normal(origin_x, width/2),
random_normal(origin_y, width/2),
4326)::Geometry(Point, 4326) AS geom,
polyid
FROM bounds,
generate_series(1,10) AS polyid
),
pts AS (
SELECT ST_GeneratePoints(ST_Buffer(geom, width/5), 20) AS geom,
polyid
FROM bounds,
polypts
)
SELECT ST_Multi(ST_Buffer(
ST_Buffer(
ST_ConcaveHull(geom, 0.3),
-2.0),
3.0))::Geometry(MultiPolygon, 4326) AS geom,
polyid
FROM pts;
``````

### Random Polygons with Voronoi Polygons

Another approach is to again start with random points, but use the Voronoi diagram as the basis of the polygon.

Generate random points in the circle.

Use the `ST_VoronoiPolygons()` function to generate polygons that subdivide the space using the random points as seeds.

Filter just the polygons that are fully contained in the originating circle.

And then use `ST_Union()` to merge those polygons into a single output shape.

Generating multiple hulls again involves stringing together the above operations with CTEs or subqueries.

Here is the full query to generate multiple polygons with the Voronoi method.
``````CREATE TABLE random_delaunay_hulls AS
WITH bounds AS (
SELECT 0 AS origin_x,
0 AS origin_y,
80 AS width,
80 AS height
),
polypts AS (
SELECT ST_Point(random_normal(origin_x, width/2),
random_normal(origin_y, width/2),
4326)::Geometry(Point, 4326) AS geom,
polyid
FROM bounds,
generate_series(1,20) AS polyid
),
voronois AS (
SELECT ST_VoronoiPolygons(
ST_GeneratePoints(ST_Buffer(geom, width/5), 10)
) AS geom,
ST_Buffer(geom, width/5) AS geom_clip,
polyid
FROM bounds,
polypts
),
cells AS (
SELECT (ST_Dump(geom)).geom, polyid, geom_clip
FROM voronois
)
SELECT ST_Union(geom)::Geometry(Polygon, 4326) AS geom, polyid
FROM cells
WHERE ST_Contains(geom_clip, geom)
GROUP BY polyid;
``````

## September 10, 2023

### Boston GIS (Regina Obe, Leo Hsu)

#### Why People care about PostGIS and Postgres and FOSS4GNA

Paul Ramsey and I recently had a Fireside chat with Path to Cituscon. Checkout the Podcast Why People care about PostGIS and Postgres. There were a surprising number of funny moments and very insightful stuff.

It was a great fireside chat but without the fireplace. We covered the birth and progression of PostGIS for the past 20 years and the trajectory with PostgreSQL. We also learned of Paul's plans to revolutionize PostGIS which was new to me. We covered many other side-line topics, like QGIS whose birth was inspired by PostGIS. We covered pgRouting and mobilitydb which are two other PostgreSQL extension projects that extend PostGIS.

We also managed to fall into the Large Language Model conversation of which Paul and I are on different sides of the fence on.

## August 15, 2023

### PostGIS Development

#### PostGIS 3.4.0

The PostGIS Team is pleased to release PostGIS 3.4.0! This version works with versions PostgreSQL 12-16, GEOS 3.6 or higher, and Proj 6.1+. To take advantage of all features, GEOS 3.12+ is needed. To take advantage of all SFCGAL features, SFCGAL 1.4.1+ is needed.

### 3.4.0

This release is a major release, it includes bug fixes since PostGIS 3.3.4 and new features.

## August 05, 2023

### PostGIS Development

#### PostGIS 3.4.0rc1

The PostGIS Team is pleased to release PostGIS 3.4.0rc1! Best Served with PostgreSQL 16 Beta2 and GEOS 3.12.0.

This version requires PostgreSQL 12 or higher, GEOS 3.6 or higher, and Proj 6.1+. To take advantage of all features, GEOS 3.12+ is needed. To take advantage of all SFCGAL features, SFCGAL 1.4.1+ is needed.

### 3.4.0rc1

This release is a release candidate of a major release, it includes bug fixes since PostGIS 3.3.4 and new features.

## July 29, 2023

### PostGIS Development

#### PostGIS 3.4.0beta2

The PostGIS Team is pleased to release PostGIS 3.4.0beta2! Best Served with PostgreSQL 16 Beta2 and GEOS 3.12.0.

This version requires PostgreSQL 12 or higher, GEOS 3.6 or higher, and Proj 6.1+. To take advantage of all features, GEOS 3.12+ is needed. To take advantage of all SFCGAL features, SFCGAL 1.4.1+ is needed.

### 3.4.0beta2

This release is a beta of a major release, it includes bug fixes since PostGIS 3.3.4 and new features.

## July 28, 2023

### PostGIS Development

#### PostGIS 3.3.4 Patch Release

The PostGIS development team is pleased to announce bug fix release 3.3.4, mostly focused on Topology fixes.

## July 14, 2023

### PostGIS Development

#### PostGIS 3.4.0beta1

The PostGIS Team is pleased to release PostGIS 3.4.0beta1! Best Served with PostgreSQL 16 Beta2 and GEOS 3.12.0.

This version requires PostgreSQL 12 or higher, GEOS 3.6 or higher, and Proj 6.1+. To take advantage of all features, GEOS 3.12+ is needed. To take advantage of all SFCGAL features, SFCGAL 1.4.1+ is needed.

### 3.4.0beta1

This release is a beta of a major release, it includes bug fixes since PostGIS 3.3.3 and new features.

## July 08, 2023

### Paul Ramsey

#### MapScaping Podcast: Pg_EventServ

Last month I got to record a couple podcast episodes with the MapScaping Podcast’s Daniel O’Donohue. One of them was on the benefits and many pitfalls of putting rasters into a relational database, and the other was about real-time events and pushing data change information out to web clients!

TL;DR: geospatial data tends to be more “visible” to end user clients, so communicating change to multiple clients in real time can be useful for “common operating” situations.

I also recorded a presentation about pg_eventserv for PostGIS Day 2022.

## June 03, 2023

### Boston GIS (Regina Obe, Leo Hsu)

#### PostGIS Bundle 3.3.3 for Windows with MobilityDB

I recently released PostGIS 3.3.3. bundle for Windows which is available on application stackbuilder and OSGeo download site for PostgreSQL 11 - 15. If you are running PostgreSQL 12 or above, you get an additional bonus extension MobilityDB which is an extension that leverages PostGIS geometry and geography types and introduces several more spatial-temporal types and functions specifically targeted for managing objects in motion.

What kind of management, think of getting the average speed a train is moving at a segment in time or collisions in time, without any long SQL code. Just use a function on the trip path, and viola. Think about storing GPS data very compactly in a singe row /column with time and being able to ask very complex questions with very little SQL. True PostGIS can do some of this using geometry with Measure (geometryM) geometry types, but you have to deal with that craziness of converting M back to timestamps, which mobilitydb temporal types automatically encode as true PostgreSQL timestamp types.

Anita Graser, of QGIS and Moving Pandas fame, has written several posts about it such as: Visualizing Trajectories with QGIS and mobilitydb and Detecting close encounters using MobilityDB 1.0.

Continue reading "PostGIS Bundle 3.3.3 for Windows with MobilityDB"

## May 30, 2023

### Crunchy Data

#### SVG Images from Postgres

PostGIS excels at storing, manipulating and analyzing geospatial data. At some point it's usually desired to convert raw spatial data into a two-dimensional representation to utilize the integrative capabilities of the human visual cortex. In other words, to see things on a map.

PostGIS is a popular backend for mapping technology, so there are many options to choose from to create maps. Data can be rendered to a raster image using a web map server like GeoServer or MapServer; it can be converted to GeoJSON or vector tiles via servers such as `pg_featureserv` and `pg_tileserv` and then shipped to a Web browser for rendering by a library such as OpenLayers, MapLibre or Leaflet; or a GIS application such as QGIS can connect to the database and create richly-styled maps from spatial queries.

What these options have in common is that they require external tools which need to be installed, configured and maintained in a separate environment. This can introduce unwanted complexity to a geospatial architecture.

This post presents a simple way to generate maps entirely within the database, with no external infrastructure required.

## SVG for the win

A great way to display vector data is to use the Scalable Vector Graphic (SVG) format. It provides rich functionality for displaying and styling geometric shapes. SVG is widely supported by web browsers and other tools.

By including CSS and Javascript it's possible to add advanced styling, custom popups, dynamic behaviour and interaction with other web page elements.

## Introducing `pg-svg`

Generating SVG "by hand" is difficult. It requires detailed knowledge of the SVG specification, and constructing a complex text format in SQL is highly error-prone. While PostGIS has had the function `ST_AsSVG` for years, it only produces the SVG path data attribute value. Much more is required to create a fully-styled SVG document.

The PL/pgSQL library `pg-svg` solves this problem! It makes it easy to convert PostGIS data into styled SVG documents. The library provides a simple API as a set of PL/pgSQL functions which allow creating an SVG document in a single SQL query. Best of all, this installs with a set of functions, nothing else required!

## Example map of US high points

The best way to understand how `pg-svg` works is to see an example. We'll create an SVG map of the United States showing the highest point in each state. The map has the following features:

• All 50 states are shown, with Alaska and Hawaii transformed to better fit the map
• States are labeled, and filled with a gradient
• High points are shown at their location by triangles whose color and size indicate the height of the high point.

The resulting map looks like this (to see tooltips open the raw image):

The SQL query to generate the map is here. It can be downloaded and run using `psql`:

``````psql -A -t -o us-highpt.svg  < us-highpt-svg.sql
``````

The SVG output `us-highpt.svg` can be viewed in any web browser.

## How it Works

Let's break the query down to see how the data is prepared and then rendered to SVG. A dataset of US states in geodetic coordinate system (WGS84, SRID = 4326) is required. We used the Natural Earth states and provinces data available here. It is loaded into a table `ne.admin_1_state_prov` with the following command:

``````shp2pgsql -c -D -s 4326 -i -I ne_10m_admin_1_states_provinces.shp ne.admin_1_state_prov | psql
``````

The query uses the SQL `WITH` construct to organize processing into simple, modular steps. We'll describe each one in turn.

### Select US state features

First, the US state features are selected from the Natural Earth boundaries table `ne.admin_1_state_prov`.

``````us_state AS (SELECT name, abbrev, postal, geom
``````

### Make a US state map

Next, the map is made more compact by realigning the far-flung states of Alaska and Hawaii.
This is done using PostGIS affine transformation functions. The states are made more proportionate using `ST_Scale`, and moved closer to the lower 48 using `ST_Translate`. The scaling is done around the location of the state high point, to make it easy to apply the same transformation to the high point feature.

``````,us_map AS (SELECT name, abbrev, postal,
-- transform AK and HI to make them fit map
CASE WHEN name = 'Alaska' THEN
ST_Translate(ST_Scale(
ST_Intersection( ST_GeometryN(geom,1), 'SRID=4326;POLYGON ((-141 80, -141 50, -170 50, -170 80, -141 80))'),
'POINT(0.5 0.75)', 'POINT(-151.007222 63.069444)'::geometry), 18, -17)
WHEN name = 'Hawaii' THEN
ST_Translate(ST_Scale(
ST_Intersection(geom, 'SRID=4326;POLYGON ((-161 23, -154 23, -154 18, -161 18, -161 23))'),
'POINT(3 3)', 'POINT(-155.468333 19.821028)'::geometry), 32, 10)
ELSE geom END AS geom
FROM us_state)
``````

### High Points of US states

Data for the highest point in each state is provided as an inline table of values:

``````,high_pt(name, state, hgt_m, hgt_ft, lon, lat) AS (VALUES
('Denali',              'AK', 6198, 20320,  -151.007222,63.069444)
,('Mount Whitney',       'CA', 4421, 14505,  -118.292,36.578583)
...
,('Britton Hill',        'FL',  105,   345,  -86.281944,30.988333)
)
``````

### Prepare High Point symbols

The next query does several things:

• translates the `lon` and `lat` location for Alaska and Hawaii high points to match the transformation applied to the state geometry
• computes the `symHeight` attribute for the height of the high point triangle symbol
• assigns a fill color value to each high point based on the height
• uses `ORDER BY` to sort the high points by latitude, so that their symbols overlap correctly when rendered
``````,highpt_shape AS (SELECT name, state, hgt_ft,
-- translate high points to match shifted states
CASE WHEN state = 'AK' THEN lon + 18
WHEN state = 'HI' THEN lon + 32
ELSE lon END AS lon,
CASE WHEN state = 'AK' THEN lat - 17
WHEN state = 'HI' THEN lat + 10
ELSE lat END AS lat,
(2.0 * hgt_ft) / 15000.0 + 0.5 AS symHeight,
CASE WHEN hgt_ft > 14000 THEN '#ffffff'
WHEN hgt_ft >  7000 THEN '#aaaaaa'
WHEN hgt_ft >  5000 THEN '#ff8800'
WHEN hgt_ft >  2000 THEN '#ffff44'
WHEN hgt_ft >  1000 THEN '#aaffaa'
ELSE '#558800'
END AS clr
FROM high_pt ORDER BY lat DESC)
``````

### Generate SVG elements

The previous queries transformed the raw data into a form suitable for rendering.
Now we get to see `pg-svg` in action! The next query generates the SVG text for each output image element, as separate records in a result set called `shapes`.

The SVG elements are generated in the order in which they are drawn - states and labels first, with high-point symbols on top. Let's break it down:

### SVG for states

The first subquery produces SVG shapes from the state geometries. The `svgShape` function produces an SVG shape element for any PostGIS geometry. It also provides optional parameters supporting other capabilities of SVG. Here `title` specifies that the state name is displayed as a tooltip, and `style` specifies the styling of the shape. Styling in SVG can be provided using properties defined in the Cascaded Style Sheets (CSS) specification. `pg-svg` provides the `svgStyle` function to make it easy to specify the names and values of CSS styling properties.

Note that the `fill` property value is a URL instead of a color specifier. This refers to an SVG gradient fill which is defined later.

The state geometry is also included in the subquery result set, for reasons which will be discussed below.

``````,shapes AS (
-- State shapes
SELECT geom, svgShape( geom,
title => name,
style => svgStyle(  'stroke', '#ffffff',
'stroke-width', 0.1::text,
'fill', 'url(#state)',
'stroke-linejoin', 'round' ) )
svg FROM us_map
``````

### SVG for state labels

Labels for state abbreviations are positioned at the point produced by the `ST_PointOnSurface` function. (Alternatively, `ST_MaximumInscribedCircle` could be used.) The SVG is generated by the `svgText` function, using the specified styling.

``````  UNION ALL
-- State names
SELECT NULL, svgText( ST_PointOnSurface( geom ), abbrev,
style => svgStyle(  'fill', '#6666ff', 'text-anchor', 'middle', 'font', '0.8px sans-serif' ) )
svg FROM us_map
``````

### SVG for high point symbols

The high point features are displayed as triangular symbols. We use the convenient `svgPolygon` function with a simple array of ordinates specifying a triangle based at the high point location, with height given by the previously computed `svgHeight` column. The title is provided for a tooltip, and the styling uses the computed `clr` attribute as the fill.

``````  UNION ALL
-- High point triangles
SELECT NULL, svgPolygon( ARRAY[ lon-0.5, -lat, lon+0.5, -lat, lon, -lat-symHeight ],
title => name || ' ' || state || ' - ' || hgt_ft || ' ft',
style => svgStyle(  'stroke', '#000000',
'stroke-width', 0.1::text,
'fill', clr  ) )
svg FROM highpt_shape
)
``````

### Produce final SVG image

The generated shape elements need to be wrapped in an `<svg>` document element. This is handled by the `svgDoc` function.

The viewable extent of the SVG data needs to be provided by the `viewbox` parameter. The most common case is to display all of the rendered data. An easy way to determine this is to apply the PostGIS `ST_Exrtent` aggregate function to the input data (this is why we included the `geom` column as well as the `svg` text column). We can include a border by enlarging the extent using the `ST_Expand` function. The function `svgViewBox` converts the PostGIS geometry for the extent into SVG format.

We also include a definition for an SVG linear gradient to be used as the fill style for the state features.

``````SELECT svgDoc( array_agg( svg ),
viewbox => svgViewbox( ST_Expand( ST_Extent(geom), 2)),
) AS svg FROM shapes;
``````

The output from `svgDoc` is a `text` value which can be used anywhere that SVG is supported.

## More to Explore

We've shown how the `pg-svg` SQL function library lets you easily generate map images from PostGIS data right in the database. This can be used as a simple ad-hoc way of visualizing spatial data. Or, it could be embedded in a larger system to automate repetitive map generation workflows.

Although SVG is a natural fit for vector data, there may be situations where producing a map as a bitmap (raster) image makes sense.
For a way of generating raster maps right in the database see this PostGIS Day 2022 presentation. This would be especially appealing where the map is displaying data stored using PostGIS raster data. It would also be possible to combine vector and raster data into a hybrid SVG/image output.

Although we've focussed on creating maps of geospatial data, SVG is often used for creating other kinds of graphics. For examples of using it to create geometric and mathematical designs see the `pg-svg` `demo` folder. Here's an image of a Lissajous knot generated by this SQL.

You could even use `pg-svg` to generate charts of non-spatial data (although this would be better handled by a more task-specific API).

Let us know if you find `pg-svg` useful, or if you have ideas for improving it!

### Crunchy Data

#### SVG Images from Postgres

PostGIS excels at storing, manipulating and analyzing geospatial data. At some point it's usually desired to convert raw spatial data into a two-dimensional representation to utilize the integrative capabilities of the human visual cortex. In other words, to see things on a map.

PostGIS is a popular backend for mapping technology, so there are many options to choose from to create maps. Data can be rendered to a raster image using a web map server like GeoServer or MapServer; it can be converted to GeoJSON or vector tiles via servers such as `pg_featureserv` and `pg_tileserv` and then shipped to a Web browser for rendering by a library such as OpenLayers, MapLibre or Leaflet; or a GIS application such as QGIS can connect to the database and create richly-styled maps from spatial queries.

What these options have in common is that they require external tools which need to be installed, configured and maintained in a separate environment. This can introduce unwanted complexity to a geospatial architecture.

This post presents a simple way to generate maps entirely within the database, with no external infrastructure required.

## SVG for the win

A great way to display vector data is to use the Scalable Vector Graphic (SVG) format. It provides rich functionality for displaying and styling geometric shapes. SVG is widely supported by web browsers and other tools.

By including CSS and Javascript it's possible to add advanced styling, custom popups, dynamic behaviour and interaction with other web page elements.

## Introducing `pg-svg`

Generating SVG "by hand" is difficult. It requires detailed knowledge of the SVG specification, and constructing a complex text format in SQL is highly error-prone. While PostGIS has had the function `ST_AsSVG` for years, it only produces the SVG path data attribute value. Much more is required to create a fully-styled SVG document.

The PL/pgSQL library `pg-svg` solves this problem! It makes it easy to convert PostGIS data into styled SVG documents. The library provides a simple API as a set of PL/pgSQL functions which allow creating an SVG document in a single SQL query. Best of all, this installs with a set of functions, nothing else required!

## Example map of US high points

The best way to understand how `pg-svg` works is to see an example. We'll create an SVG map of the United States showing the highest point in each state. The map has the following features:

• All 50 states are shown, with Alaska and Hawaii transformed to better fit the map
• States are labeled, and filled with a gradient
• High points are shown at their location by triangles whose color and size indicate the height of the high point.

The resulting map looks like this (to see tooltips open the raw image):

The SQL query to generate the map is here. It can be downloaded and run using `psql`:

``````psql -A -t -o us-highpt.svg  < us-highpt-svg.sql
``````

The SVG output `us-highpt.svg` can be viewed in any web browser.

## How it Works

Let's break the query down to see how the data is prepared and then rendered to SVG. A dataset of US states in geodetic coordinate system (WGS84, SRID = 4326) is required. We used the Natural Earth states and provinces data available here. It is loaded into a table `ne.admin_1_state_prov` with the following command:

``````shp2pgsql -c -D -s 4326 -i -I ne_10m_admin_1_states_provinces.shp ne.admin_1_state_prov | psql
``````

The query uses the SQL `WITH` construct to organize processing into simple, modular steps. We'll describe each one in turn.

### Select US state features

First, the US state features are selected from the Natural Earth boundaries table `ne.admin_1_state_prov`.

``````us_state AS (SELECT name, abbrev, postal, geom
``````

### Make a US state map

Next, the map is made more compact by realigning the far-flung states of Alaska and Hawaii.
This is done using PostGIS affine transformation functions. The states are made more proportionate using `ST_Scale`, and moved closer to the lower 48 using `ST_Translate`. The scaling is done around the location of the state high point, to make it easy to apply the same transformation to the high point feature.

``````,us_map AS (SELECT name, abbrev, postal,
-- transform AK and HI to make them fit map
CASE WHEN name = 'Alaska' THEN
ST_Translate(ST_Scale(
ST_Intersection( ST_GeometryN(geom,1), 'SRID=4326;POLYGON ((-141 80, -141 50, -170 50, -170 80, -141 80))'),
'POINT(0.5 0.75)', 'POINT(-151.007222 63.069444)'::geometry), 18, -17)
WHEN name = 'Hawaii' THEN
ST_Translate(ST_Scale(
ST_Intersection(geom, 'SRID=4326;POLYGON ((-161 23, -154 23, -154 18, -161 18, -161 23))'),
'POINT(3 3)', 'POINT(-155.468333 19.821028)'::geometry), 32, 10)
ELSE geom END AS geom
FROM us_state)
``````

### High Points of US states

Data for the highest point in each state is provided as an inline table of values:

``````,high_pt(name, state, hgt_m, hgt_ft, lon, lat) AS (VALUES
('Denali',              'AK', 6198, 20320,  -151.007222,63.069444)
,('Mount Whitney',       'CA', 4421, 14505,  -118.292,36.578583)
...
,('Britton Hill',        'FL',  105,   345,  -86.281944,30.988333)
)
``````

### Prepare High Point symbols

The next query does several things:

• translates the `lon` and `lat` location for Alaska and Hawaii high points to match the transformation applied to the state geometry
• computes the `symHeight` attribute for the height of the high point triangle symbol
• assigns a fill color value to each high point based on the height
• uses `ORDER BY` to sort the high points by latitude, so that their symbols overlap correctly when rendered
``````,highpt_shape AS (SELECT name, state, hgt_ft,
-- translate high points to match shifted states
CASE WHEN state = 'AK' THEN lon + 18
WHEN state = 'HI' THEN lon + 32
ELSE lon END AS lon,
CASE WHEN state = 'AK' THEN lat - 17
WHEN state = 'HI' THEN lat + 10
ELSE lat END AS lat,
(2.0 * hgt_ft) / 15000.0 + 0.5 AS symHeight,
CASE WHEN hgt_ft > 14000 THEN '#ffffff'
WHEN hgt_ft >  7000 THEN '#aaaaaa'
WHEN hgt_ft >  5000 THEN '#ff8800'
WHEN hgt_ft >  2000 THEN '#ffff44'
WHEN hgt_ft >  1000 THEN '#aaffaa'
ELSE '#558800'
END AS clr
FROM high_pt ORDER BY lat DESC)
``````

### Generate SVG elements

The previous queries transformed the raw data into a form suitable for rendering.
Now we get to see `pg-svg` in action! The next query generates the SVG text for each output image element, as separate records in a result set called `shapes`.

The SVG elements are generated in the order in which they are drawn - states and labels first, with high-point symbols on top. Let's break it down:

### SVG for states

The first subquery produces SVG shapes from the state geometries. The `svgShape` function produces an SVG shape element for any PostGIS geometry. It also provides optional parameters supporting other capabilities of SVG. Here `title` specifies that the state name is displayed as a tooltip, and `style` specifies the styling of the shape. Styling in SVG can be provided using properties defined in the Cascaded Style Sheets (CSS) specification. `pg-svg` provides the `svgStyle` function to make it easy to specify the names and values of CSS styling properties.

Note that the `fill` property value is a URL instead of a color specifier. This refers to an SVG gradient fill which is defined later.

The state geometry is also included in the subquery result set, for reasons which will be discussed below.

``````,shapes AS (
-- State shapes
SELECT geom, svgShape( geom,
title => name,
style => svgStyle(  'stroke', '#ffffff',
'stroke-width', 0.1::text,
'fill', 'url(#state)',
'stroke-linejoin', 'round' ) )
svg FROM us_map
``````

### SVG for state labels

Labels for state abbreviations are positioned at the point produced by the `ST_PointOnSurface` function. (Alternatively, `ST_MaximumInscribedCircle` could be used.) The SVG is generated by the `svgText` function, using the specified styling.

``````  UNION ALL
-- State names
SELECT NULL, svgText( ST_PointOnSurface( geom ), abbrev,
style => svgStyle(  'fill', '#6666ff', 'text-anchor', 'middle', 'font', '0.8px sans-serif' ) )
svg FROM us_map
``````

### SVG for high point symbols

The high point features are displayed as triangular symbols. We use the convenient `svgPolygon` function with a simple array of ordinates specifying a triangle based at the high point location, with height given by the previously computed `svgHeight` column. The title is provided for a tooltip, and the styling uses the computed `clr` attribute as the fill.

``````  UNION ALL
-- High point triangles
SELECT NULL, svgPolygon( ARRAY[ lon-0.5, -lat, lon+0.5, -lat, lon, -lat-symHeight ],
title => name || ' ' || state || ' - ' || hgt_ft || ' ft',
style => svgStyle(  'stroke', '#000000',
'stroke-width', 0.1::text,
'fill', clr  ) )
svg FROM highpt_shape
)
``````

### Produce final SVG image

The generated shape elements need to be wrapped in an `<svg>` document element. This is handled by the `svgDoc` function.

The viewable extent of the SVG data needs to be provided by the `viewbox` parameter. The most common case is to display all of the rendered data. An easy way to determine this is to apply the PostGIS `ST_Exrtent` aggregate function to the input data (this is why we included the `geom` column as well as the `svg` text column). We can include a border by enlarging the extent using the `ST_Expand` function. The function `svgViewBox` converts the PostGIS geometry for the extent into SVG format.

We also include a definition for an SVG linear gradient to be used as the fill style for the state features.

``````SELECT svgDoc( array_agg( svg ),
viewbox => svgViewbox( ST_Expand( ST_Extent(geom), 2)),
) AS svg FROM shapes;
``````

The output from `svgDoc` is a `text` value which can be used anywhere that SVG is supported.

## More to Explore

We've shown how the `pg-svg` SQL function library lets you easily generate map images from PostGIS data right in the database. This can be used as a simple ad-hoc way of visualizing spatial data. Or, it could be embedded in a larger system to automate repetitive map generation workflows.

Although SVG is a natural fit for vector data, there may be situations where producing a map as a bitmap (raster) image makes sense.
For a way of generating raster maps right in the database see this PostGIS Day 2022 presentation. This would be especially appealing where the map is displaying data stored using PostGIS raster data. It would also be possible to combine vector and raster data into a hybrid SVG/image output.

Although we've focussed on creating maps of geospatial data, SVG is often used for creating other kinds of graphics. For examples of using it to create geometric and mathematical designs see the `pg-svg` `demo` folder. Here's an image of a Lissajous knot generated by this SQL.

You could even use `pg-svg` to generate charts of non-spatial data (although this would be better handled by a more task-specific API).

Let us know if you find `pg-svg` useful, or if you have ideas for improving it!

## May 29, 2023

### PostGIS Development

#### PostGIS 3.3.3, 3.2.5, 3.1.9, 3.0.9 Patch Releases

The PostGIS development team is pleased to provide bug fixes and performance enhancements 3.3.3, 3.2.5, 3.1.9 and 3.0.9 for the 3.3, 3.2, 3.1, and 3.0 stable branches.

## May 24, 2023

### Paul Ramsey

#### Keynote @ CUGOS Spring Fling

Last month I was invited to give a keynote talk at the CUGOS Spring Fling, a delightful gathering of “Cascadia Users of Open Source GIS” in Seattle. I have been speaking about open source economics at FOSS4G conferences more-or-less every two years, since 2009, and took this opportunity to somewhat revisit the topics of my 2019 FOSS4GNA keynote.

If you liked the video and want to use the materials, the slides are available here under CC BY.

## May 17, 2023

### Paul Ramsey

#### MapScaping Podcast: Rasters and PostGIS

Last month I got to record a couple podcast episodes with the MapScaping Podcast’s Daniel O’Donohue. One of them was on the benefits and many pitfalls of putting rasters into a relational database, and it is online now!

TL;DR: most people think “put it in a database” is a magic recipe for: faster performance, infinite scalability, and easy management.

Where the database is replacing a pile of CSV files, this is probably true.

Where the database is replacing a collection of GeoTIFF imagery files, it is probably false. Raster in the database will be slower, will take up more space, and be very annoying to manage.

So why do it? Start with a default, “don’t!”, and then evaluate from there.

For some non-visual raster data, and use cases that involve enriching vectors from raster sources, having the raster co-located with the vectors in the database can make working with it more convenient. It will still be slower than direct access, and it will still be painful to manage, but it allows use of SQL as a query language, which can give you a lot more flexibility to explore the solution space than a purpose built data access script might.

There’s some other interesting tweaks around storing the actual raster data outside the database and querying it from within, that I think are the future of “raster in (not really in) the database”, listen to the episode to learn more!

## March 02, 2023

### Crunchy Data

#### Geocoding with Web APIs in Postgres

Geocoding is the process of taking addresses or location information and getting the coordinates for that location. Anytime you route a new location or look up a zip code, the back end is geocoding the location and then using the geocode inside other PostGIS functions to give you the routes, locations, and other data you asked for.

PostGIS comes equipped with an easy way to use the US Census data with the Tiger geocoder. Using the Tiger geocoder requires downloading large amounts of census data and in space-limited databases, this may not be ideal. Using a geocoding web API service can be a space saving solution in these cases.

I am going to show you how to set up a really quick function using `plpython3u` to hit a web service geocoder every time that we get a new row in the database.

### Installing plpython3u

The plpython3u extension comes with Crunchy Bridge or you can add it to your database. To get started run the following:

``````CREATE EXTENSION  plpython3u;
``````

### Creating a function to geocode addresses

In this example, I'll use the US census geocoding API as our web service, and build a function to geocode addresses based on that.

The function puts together parameters to hit the census geocoding API and then parses the resulting object, and returns a geometry:

``````CREATE OR REPLACE FUNCTION geocode(address text)
RETURNS geometry
AS \$\$
import requests
try:
r = requests.get(base_geocode, params = payload)
lon = coords['x']
lat = coords['y']
geom = f'SRID=4326;POINT({lon} {lat})'
except Exception as e:
plpy.notice(f'error: {e.message}')
geom = None
return geom
\$\$
LANGUAGE 'plpython3u';
``````

Using this function to geocode Crunchy Data's headquarters:

``````SELECT ST_AsText(geocode('162 Seven Farms Drive Charleston, SC 29492'));
``````

### Deploying this function for new data

But what if we want to automatically run this every time an address is inserted into a table? Let's say we have a table with a field ID, an address, and a point that we want to auto-populate on inserts.

``````CREATE TABLE addresses (
fid SERIAL PRIMARY KEY,
geom GEOMETRY(POINT, 4326)
);
``````

We can make use of a Postgres trigger to add the geocode before every insert! Triggers are a very powerful way to leverage built in functions to automatically transform your data as it enters the database, and this particular case is a great demo for them!

``````CREATE OR REPLACE FUNCTION add_geocode()
RETURNS trigger AS
\$\$
DECLARE
loc geometry;
BEGIN
NEW.geom = loc;
RETURN NEW;
END;
\$\$
LANGUAGE plpgsql;
CREATE TRIGGER update_geocode BEFORE INSERT ON addresses
FOR EACH ROW EXECUTE FUNCTION add_geocode();
``````

Now when running an insert, the value is automatically geocoded!

``````INSERT INTO addresses(address) VALUES ('415 Mission St, San Francisco, CA 94105');

-----+-----------------------------------------+----------------------------------------------------
1 | 415 Mission St, San Francisco, CA 94105 | 0101000020E610000097CD0E2B66995EC0BB004B2729E54240
``````

### Summary

If you’re space limited, using a web API based geocoder might be the way to go. Using a geocoder function with triggers on new row inserts will get you geocoded addresses in a snap.

## February 23, 2023

### Crunchy Data

#### Postgres Raster Query Basics

In geospatial terminology, a "raster" is a cover of an area divided into a uniform gridding, with one or more values assigned to each grid cell.

A "raster" in which the values are associated with red, green and blue bands might be a visual image. The rasters that come off the Landsat 7 earth observation satellite have eight bands: red, green, blue, near infrared, shortwave infrared, thermal, mid-infrared and panchromatic.

## Database Rasters

Working with raster data via SQL is a little counter-intuitive: rasters don't neatly fit the relational model the way vector geometries do. A table of parcels where one column is the geometry and the others are the owner name, address, and tax roll id makes sense. How should a raster fit into a table? As a row for every pixel? For every scan row? What other values should be associated with each row?

There is no clean relationship between "real world objects" and the database representation of a raster, because a raster has nothing to say about objects, it is just a collection of measurements.

We can squeeze rasters into the database, but doing so makes working with the data more complex. Before loading data, we need to enable PostGIS and the raster module.

``````CREATE EXTENSION postgis;
CREATE EXTENSION postgis_raster;
``````

For this example, we will load raster data for a "digital elevation model" (DEM), a raster with just one band, the elevation at each pixel.

Using the SRTM Tile Grabber I downloaded one tile of old SRTM data. Then using the `gdalinfo` utility, read out the metadata about the file.

``````wget https://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_12_03.zip
unzip srtm_12_03.zip
gdalinfo srtm_12_03.tif
``````

• The coordinate system of the data is WGS 84.
• The pixel size is 3 arc-seconds.
• The pixel type is `Int16`, so two bytes per pixel.

Knowing that, I can build a `raster2pgsql` call to load the data into a raster table.

``````raster2pgsql \
-I \                 # create a spatial index on the column
-s 4326 \            # use 4326 (WGS 84) as the spatial reference for the raster
-t 32x32 \           # tile the raster into 32 by 32 pixel tiles
srtm_12_03.tif | \   # name of the raster
psql dem             # target database connection
``````

Once loaded the raster table looks like this on a map.

And it looks like this in the database.

``````Table "public.srtm_12_03"

Column |  Type
--------+---------
rid    | integer
rast   | raster
Indexes:
"srtm_12_03_pkey" PRIMARY KEY, btree (rid)
"srtm_12_03_st_convexhull_idx" gist (st_convexhull(rast))
``````

It's a pretty boring table! Just a bunch of binary raster tiles and a unique key for each.

``````-- 29768 rows
SELECT Count(*) FROM srtm_12_03;
``````

Those binary raster tiles aren't just opaque blobs though, we can look inside them with the right functions. Here we get a summary of all the raster tiles in the table.

``````SELECT (ST_SummaryStatsAgg(rast, 1, true)).* FROM srtm_12_03;
``````
``````count  | 28966088
sum    | 20431360140
mean   | 705.3544869434907
stddev | 561.252765463607
min    | -291
max    | 4371
``````

## Tiles, Tiles, Tiles

Remember when we loaded the data with `raster2pgsql` we specified a "tile size" of 32 by 32 pixels? This has a number of implications.

• First, at 32x32, a tile of 2-byte `Int16` data like our DEM will take up 2048 bytes. This is small enough to fit in the database page size, which means the data will not end up stored in a side table by the TOAST subsystem that handles large row values.
• Second, our input file had 6000x6000 pixels, which is enough pixels to generate 35156 32x32 tiles. Our table only has 29768 rows, because `raster2pgsql` does not generate tiles when the contents are all "no data" pixels (as the DEM data is over the ocean).

The loaded data looks like this.

Notice how small each tile is. As a general rule, when working with raster data queries,

• the first step will be to efficiently find the relevant tile(s);
• the second step will be to use the tiles to find the answer you want.

Finding tiles efficiently means using spatial index, and the spatial index definition as we saw above is this:

``````"srtm_12_03_st_convexhull_idx" gist (st_convexhull(rast))
``````

This is a functional index, which means in order to access it, we need to copy the functional part: `st_convechull(rast)` when forming our query.

The ST_ConvexHull(raster) converts a raster tile into a polygon defining the boundary of the tile. When querying raster tables, you will use this function a great deal to convert rasters into polygons suitable for querying a spatial index.

## Point Query

The simplest raster query is to take a point, and find the value of the raster under that point.

Here is a point table with one point in it:

``````CREATE TABLE mappoint AS
SELECT ST_Point(-123.7273, 47.8467, 4326)::geometry(Point, 4326) AS geom,
1 AS fid;
``````

The nice thing about points is that they only hit one tile at a time. So we don't have to think too hard about what to do with our tile sets.

``````SELECT ST_Value(srtm.rast, pt.geom)
FROM srtm_12_03 srtm
JOIN mappoint pt
ON ST_Intersects(pt.geom, ST_ConvexHull(srtm.rast))
``````
`````` st_value
----------
1627
``````

Here we use the ST_ConvexHull(raster) function to get access to our spatial index on the raster table, and the ST_Intersects(geom, geom) function to test the condition.

The ST_Value(raster, geom) function reads the pixel value from the raster at the location of the point.

## Polygon Query

Summarizing rasters under polygons is more involved than reading point values, because polygons will frequently overlap multiple tiles, so you have to think in terms of "sets of raster tiles" instead of "the raster" when building your query.

• Then we will summarize the clipped tiles, to get the average elevation in our polygon.

The final complete query looks like this.

``````WITH clipped_tiles AS (
SELECT ST_Clip(srtm.rast, ply.geom) AS rast, srtm.rid
FROM srtm_12_03 srtm
JOIN mappoly ply
ON ST_Intersects(ply.geom, ST_ConvexHull(srtm.rast))
)
SELECT (ST_SummaryStatsAgg(rast, 1, true)).*
FROM clipped_tiles;
``````
``````count  | 362369
sum    | 388175193
mean   | 1071.2152336430545
stddev | 441.7982032761408
min    | 108
max    | 2374
``````

## Conclusion

Working with database rasters analytically can be challenging, particularly if you are used to thinking about them as single, unitary coverages. Remember to apply the basic rules of database rasters:

• Make sure to use ST_ConvexHull(raster) to drive a spatial index filter.
• Assemble the chips as needed to answer the query (you might need to use the ST_Union(raster) aggregate).
• Carry out your actual raster query.

## February 10, 2023

### Crunchy Data

#### Temporal Filtering in pg_featureserv with CQL

In a previous post we announced the CQL filtering capability in `pg_featureserv`. It provides powerful functionality for attribute and spatial querying of data in PostgreSQL and PostGIS.

Another important datatype which is often present in datasets is temporal. Temporal datasets contain attributes which are dates or timestamps. The CQL standard defines some special-purpose syntax to support temporal filtering. This allows `pg_featureserv` to take advantage of the extensive capabilities of PostgreSQL for specifying queries against time-valued attributes. This post in the CQL series will show some examples of temporal filtering in `pg_featureserv`.

## CQL Temporal filters

Temporal filtering in CQL is provided using temporal literals and conditions.

Temporal literal values may be dates or timestamps:

``````2001-01-01
2010-04-23T01:23:45
``````

Note: The temporal literal syntax is based on an early version of the OGC API Filter and CQL standard. The current draft CQL standard has a different syntax: `DATE('1969-07-20')` and `TIMESTAMP('1969-07-20T20:17:40Z')`. It also supports intervals: `INTERVAL('1969-07-16', '1969-07-24')`. A subsequent version of `pg_featureserv` will support this syntax as well.

Temporal conditions allow time-valued properties and literals to be compared via the standard boolean comparison operators `<`,`>`,`<=`,`>=`,`=`,`<>` and the `BETWEEN..AND` operator:

``````start_date >= 2001-01-01
event_time BETWEEN 2010-04-22T06:00 AND 2010-04-23T12:00
``````

The draft CQL standard provides dedicated temporal operators, such as `T_AFTER`, `T_BEFORE`, `T_DURING`, etc. A future version of `pg_featureserv` will likely provide these operators.

## Publishing Historical Tropical Storm tracks

We'll demonstrate temporal filters using a dataset with a strong time linkage: tracks of tropical storms (or hurricanes). There is a dataset of Historical Tropical Storm Tracks available here.

The data requires some preparation. It is stored as a set of records of line segments representing 6-hour long sections of storm tracks. To provide simpler querying we will model the data using a single record for each storm, with a line geometry showing the entire track and attributes for the start and end time for the track.

The data is provided in Shapefile format. As expected for a worldwide dataset, it is in the WGS84 geodetic coordinate system (lat/long). In PostGIS this common Spatial Reference System is assigned an identifier (SRID) of 4326.

The PostGIS `shp2pgsql` utility can be used to load the dataset into a spatial table called `trop_storm_raw`. The `trop_storm_raw` table is a temporary staging table allowing the raw data to be loaded and made available for the transformation phase of data preparation.

``````shp2pgsql -c -D -s 4326 -i -I -W LATIN1 "Historical Tropical Storm Tracks.shp" public.trop_storm_raw | psql -d database
``````

The options used are:

• `-c` - create a new table
• `-D` - use PostgreSQL dump format to load the data
• `-s` - specify the SRID of 4326
• `-i` - use 32-bit integers
• `-I` - create a GIST index on the geometry column (this is not strictly necessary, since this is just a temporary staging table)
• `-W` - specifies the encoding of the input attribute data in the DBF file

Next, create the table having the desired data model:

``````CREATE TABLE public.trop_storm (
btid int PRIMARY KEY,
name text,
wind_kts numeric,
pressure float8,
basin text,
time_start timestamp,
time_end timestamp,
geom geometry(MultiLineString, 4326)
);
``````

It's good practice to add comments to the table and columns. These will be displayed in the `pg_featureserv` Web UI.

``````COMMENT ON TABLE public.trop_storm IS 'This is my spatial table';
COMMENT ON COLUMN public.trop_storm.geom IS 'Storm track LineString';
COMMENT ON COLUMN public.trop_storm.name IS 'Name assigned to storm';
COMMENT ON COLUMN public.trop_storm.btid IS 'Id of storm';
COMMENT ON COLUMN public.trop_storm.wind_kts IS 'Maximum wind speed in knots';
COMMENT ON COLUMN public.trop_storm.pressure IS 'Minumum pressure in in millibars';
COMMENT ON COLUMN public.trop_storm.basin IS 'Basin in which storm occured';
COMMENT ON COLUMN public.trop_storm.time_start IS 'Timestamp of storm start';
COMMENT ON COLUMN public.trop_storm.time_end IS 'Timestamp of storm end';
``````

Now the power of SQL can be used to transform the raw data into the simpler data model. The track sections can be combined into single tracks with a start and end time using the following query.

• The original data represents the track sections as `MultiLineString`s with single elements. The element is extracted using `ST_GeometryN` so that the result of aggregating them using `ST_Collect` is a `MultiLineString`, not a `GeometryCollection`. (An alternative is to aggregate into a GeometryCollection and use `ST_CollectionHomogenize` to reduce it to a `MultiLineString`.)
• The final `ST_Multi` ensures that all tracks are stored as `MultiLineStrings`, as required by the type constraint on the `geom` column.
• the filter condition `time_end - time_start < '1 year'::interval` removes tracks spanning the International Date Line.
``````WITH data AS (
SELECT btid, name, wind_kts, pressure, basin, geom,
make_date(year::int, month::int, day::int) + ad_time::time AS obs_time
FROM trop_storm_raw ORDER BY obs_time
),
tracks AS (
SELECT btid,
MAX(name) AS name,
MAX(wind_kts) AS wind_kts,
MAX(pressure) AS pressure,
MAX(basin) AS basin,
MIN(obs_time) AS time_start,
MAX(obs_time) AS time_end,
ST_Multi( ST_LineMerge( ST_Collect( ST_GeometryN(geom, 1)))) AS geom
FROM data GROUP BY btid
)
INSERT INTO trop_storm
SELECT * FROM tracks WHERE time_end - time_start < '1 year'::interval;
``````

This is a small dataset, and `pg_featureserv` does not require one, but as per best practice we can create a spatial index on the geometry column:

``````CREATE INDEX trop_storm_gix ON public.trop_storm USING GIST ( geom );
``````

Once the `trop_storm` table is created and populated, it can be published in `pg_featureserv`. Issuing the following request in a browser shows the feature collection in the Web UI:

``````http://localhost:9000/collections.html
``````

``````http://localhost:9000/collections/public.trop_storm.html
``````

The dataset can be viewed using `pg_featureserv`'s built-in map viewer (note that to see all 567 records displayed it is probably necessary to increase the limit on the number of response features):

``````http://localhost:9000/collections/public.trop_storm/items.html?limit=1000
``````

## Querying by Time Range

That's a lot of storm tracks. It would be easier to visualize a smaller number of tracks. A natural way to subset the data is by querying over a time range. Let's retrieve the storms between the start of 2005 and the end of 2009. This is done by adding a `filter` parameter with a CQL expression against the dataset temporal property `time_start` (storms typically do not span the start of years). To query values lying between a range of times it is convenient to use the `BETWEEN` operator. The filter condition is `time_start BETWEEN 2005-01-01 AND 2009-12-31`. The full request is:

``````http://localhost:9000/collections/public.trop_storm/items.html?filter=time_start BETWEEN 2005-01-01 AND 2009-12-31&limit=100
``````

Submitting this query produces a result with 68 tracks:

## Querying by Time and Space

Temporal conditions can be combined with other kinds of filters. For instance, we can execute a spatio-temporal query by using a temporal condition along with a spatial condition. In this example, we query the storms which occurred in 2005 and after in Florida. The temporal condition is expressed as `time_start > 2005-01-01`.

The spatial condition uses the `INTERSECTS` predicate to test whether the line geometry of a storm track intersects a polygon representing the (simplified) coastline of Florida. The polygon is provided as a geometry literal using WKT. (For more information about spatial filtering with CQL in `pg_featureserv` see this blog post.)

``````POLYGON ((-81.4067 30.8422, -79.6862 25.3781, -81.1609 24.7731, -83.9591 30.0292, -85.2258 29.6511, -87.5892 29.9914, -87.5514 31.0123, -81.4067 30.8422))
``````

Putting these conditions together in a boolean expression using `AND`, the request to retrieve the desired tracks from `pg_featureserv` is:

``````http://localhost:9000/collections/public.trop_storm/items.html?filter=time_start > 2005-01-01 AND INTERSECTS(geom, POLYGON ((-81.4067 30.8422, -79.6862 25.3781, -81.1609 24.7731, -83.9591 30.0292, -85.2258 29.6511, -87.5892 29.9914, -87.5514 31.0123, -81.4067 30.8422)) )&limit=100
``````

This query produces a result with only 9 tracks, all of which cross Florida:

## Try it yourself!

CQL temporal filtering is included in the forthcoming `pg_featureserv` Version 1.3. But you can try it out now by downloading the latest build. Let us know what use cases you find for CQL temporal filtering! Crunchy Data offers full managed PostGIS in the Cloud, with Container apps to run pg_featureserv. Try it today.

## February 02, 2023

### Paul Ramsey

#### When Proj Grid-Shifts Disappear

Last week a user noted on the postgis-users list (paraphrase):

I upgraded from PostGIS 2.5 to 3.3 and now the results of my coordinate transforms are wrong. There is a vertical shift between the systems I’m using, but my vertical coordinates are unchanged.

Hmmm.

PostGIS gets all its coordinate reprojection smarts from the proj library. The user’s query looked like this:

``````SELECT ST_AsText(
ST_Transform('SRID=7405;POINT(545068 258591 8.51)'::geometry,
4979
));
``````

“We just use proj” is a lot less certain and stable an assertion than it appears on the surface. In fact, PostGIS “just uses proj” for proj versions from 4.9 all the way up to 9.2, and there has been a lot of change to the proj library over that sweep of releases.

• The API radically changed around proj version 6, and that required a major rework of how PostGIS called the library.
• The way proj dereferenced spatial reference identifiers into reprojection algorithms changed around then too (it got much slower) which required more changes in how we interacted with the library.
• More recently the way proj handled “transformation grids” changed, which was the source of the complaint.

## Exploring Proj

The first thing to do in debugging this “PostGIS problem” was to establish if it was in fact a PostGIS problem, or a problem in proj. There are commandline tools in proj to query what pipelines the system will use for a transform, and what the effect on coordinates will be, so I can take PostGIS right out of the picture.

We can run the query on the commandline:

``````echo 545068 258591 8.51 | cs2cs 'EPSG:7405' 'EPSG:4979'
``````

Which returns:

``````52d12'23.241"N  0d7'17.603"E 8.510
``````

So directly using proj we are seeing the same problem as in PostGIS SQL: no change in the vertical dimension, it goes in at 8.51 and comes out at 8.51. So the problem is not PostGIS, is is proj.

## Transformation Grids

Cartographic transformations are nice deterministic functions, they take in a longitude and latitude and spit out an X and a Y.

``````(x,y) = f(theta, phi)
(theta, phi) = finv(x, y)
``````

But not all transformations are cartographic transformations, some are transformation between geographic reference systems. And many of those are lumpy and kind of random.

For example, the North American 1927 Datum (NAD27) was built from classic survey techniques, starting from the “middle” (Kansas) and working outwards, chain by chain, sighting by sighting. The North American 1983 Datum (NAD83) was built with the assistance of the first GPS units. The accumulated errors of survey over distance are not deterministic, they are kind of lumpy and arbitrary. So the transformation from NAD27 to NAD83 is also kind of lumpy and arbitrary.

How do you represent lumpy and arbitrary transformations? With a grid! The grid says “if your observation falls in this cell, adjust it this much, in this direction”.

Here’s a picture of the horizontal deviations in the NADCON grid.

Transformations between vertical systems also frequently require a grid.

So what does this have to do with our bug? Well, the way proj gets its grids changed in version 7.

## Grid History

Proj grids have always been a bit of an outlier. They are much larger than just the source code is. They are localized in interest (someone in New Zealand probably doesn’t need European grids), not everyone needs all the grids. So historically they were distributed in zip files separately from the code.

This is all well and good, but software packagers wanted to provide a good “works right at install” experience to their end users, so they bundled up the grids into the proj packages.

As more and more people consumed proj via packages and software installers, the fact that the grids were “separate” from proj became invisible to the end users: they just download software and it works.

This was fine while the collection of grids was a manageable size. But it is not manageable any more.

In working through the GDALBarn project to improve proj, Even Roualt decided to find all the grids that various agencies had released for various places. It turns out, there are a lot more grids than proj previously bundled. Gigabytes more.

## Grid CDN

Simply distributing the whole collection of grids as a default with proj was not going to work anymore.

So for proj 7, Even proposed moving to a download-on-demand model for proj grids. If a transformation request requires a grid, proj will attempt to download the necessary grid from the internet, and save it in a local cache.

Now everyone can get the very best possible tranformation between system, everywhere on the globe, as long as they are connected to the internet.

## Turn It On!

Except… the network grid feature is not turned on by default! So for versions of proj higher than 7, the software ships with no grids, and the software won’t check for grids on the network… until you turn on the feature!

There are three ways to turn it on, I’m going to focus on the `PROJ_NETWORK` environment variable because it’s easy to toggle. Let’s look at the proj transformation pipeline from our original bug.

``````projinfo -s EPSG:7405 -t EPSG:4979
``````

The `projinfo` utility reads out all the possible transformation pipelines, in order of desirability (accuracy) and shows what each step is. Here’s the most desireable pipeline for our transform.

``````+proj=pipeline
+step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
+y_0=-100000 +ellps=airy
+step +proj=hgridshift +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif
+step +proj=vgridshift +grids=uk_os_OSGM15_GB.tif +multiplier=1
+step +proj=axisswap +order=2,1
``````

This transform actually uses two grids! A horizontal and a vertical shift. Let’s run the shift with the network explicitly turned off.

``````echo 545068 258591 8.51 | PROJ_NETWORK=OFF cs2cs 'EPSG:7405' 'EPSG:4979'

52d12'23.241"N  0d7'17.603"E 8.510
``````

Same as before, and the elevation value is unchanged. Now run with `PROJ_NETWORK=ON`.

``````echo 545068 258591 8.51 | PROJ_NETWORK=ON cs2cs 'EPSG:7405' 'EPSG:4979'

52d12'23.288"N  0d7'17.705"E 54.462
``````

Note that the horizontal and vertical results are different with the network, because we now have access to both grids, via the CDN.

## No Internet?

If you have no internet, how do you do grid shifted transforms? Well, much like in the old days of proj, you have to manually grab the grids you need. Fortunately there is a utility for that now that makes it very easy: projsync.

``````projsync --all
``````

``````projsync --bbox 2,49,2,49
``````

If you don’t want to turn on network access via the environment variable, you can hunt down the `proj.ini` file and flip the `network = on` variable.

## January 12, 2023

### Crunchy Data

#### Fun with Letters in PostGIS 3.3!

Working at Crunchy Data on the spatial team, I'm always looking for new features and fun things to show on live demos. I recently started playing around with `ST_Letters` and wanted to jot down some quick code samples for playing around with this feature, introduced in PostGIS 3.3. These examples are super easy to use, they don't need any data!

The screenshots shown below came from pgAdmin's geometry viewer and will also work with other query GUI tools like QGIS or DBeaver.

## `ST_Letters`

Here's a simple example to get started with `ST_Letters`. This will work on any Postgres database, running the PostGIS extension version 3.3+.

`Select ST_Letters('PostGIS');`

It's also possible to overlay letters on a map, just like any other polygon. Since the default for `ST_Letters` results in a polygon starting at the baseline at the origin of the chosen projection, with a maximum height of 100 "units" (from the bottom of the descenders to the tops of the capitals).

That's not ideal. We need a way to both move it and resize it.

First, we want to make a point in the middle of San Francisco in order to serve as a centroid for where we want to move the letters, and we also want to rescale the letters in order to approximately fit over the City of San Francisco. Using the formula for converting units in WGS84 to meters, 0.001 works approximately well enough to fit over the San Francisco Bay Area.

Next we use `ST_Translate` in order to move the letters from the top of the map to fit over the Bay Area. Finally, mostly because it looks cool, we use `ST_Rotate` to rotate the polygon 45 degrees.

``````WITH
san_fran_pt AS (
SELECT ST_Point(-122.48, 37.758, 4326) AS geom),
letters AS (
SELECT ST_Scale(ST_SetSRID(
ST_Letters('San Francisco'), 4326),
0.001, 0.001) AS geom),
letters_geom AS (
SELECT ST_Translate(
letters.geom,
ST_X(san_fran_pt.geom) - ST_X(ST_Centroid(letters.geom)),
ST_Y(san_fran_pt.geom) - ST_Y(ST_Centroid(letters.geom))
) AS geom
FROM letters, san_fran_pt
)
SELECT ST_Rotate(geom, -pi() / 4, ST_Centroid(geom))
FROM letters_geom;
``````

## `ST_ConcaveHull` demo'd with `ST_Letters`

A great use case for `ST_Letters` is for demoing PostGIS functions. In this post, I'm going to demo the function `ST_ConcaveHull`, which creates a concave polygon which encloses the vertices of a target geometry. `ST_ConcaveHull` was recently updated in PostGIS 3.3.0, in order to use GEOS 3.11, which makes the input parameters easier to understand and results in a large speed upgrade. Here's a short demo of how different parameters of `param_pctconvex` and `param_allow_holes` for `ST_ConcaveHull` operate on points generated by `ST_GeneratePoints` and `ST_Letters`.

First, let's generate a table of randomly generated points that fill in the letters in 'postgis'.

``````CREATE TABLE public.word_pts AS
WITH word AS (
SELECT ST_Letters('postgis') AS geom
),
letters AS ( -- dump letter multipolygons into individual polygons
SELECT (ST_Dump(word.geom)).geom
FROM word
)
SELECT
letters.geom AS polys,
ST_GeneratePoints(letters.geom, 100) AS pts
FROM letters;
``````

`SELECT pts FROM word_pts.pts`

Then, we set the convexity to a fairly high parameter (`param_pctconvex=0.75`, indicating a highly convex shape), and don't allow there to be holes in the shape (`param_allow_holes=false`)

``````SELECT ST_ConcaveHull(pts, 0.75, false) FROM word_pts;
``````

Doesn't look much like 'postgis'!

Next, we reduce the convexity, but don't allow holes in the shape.

``````SELECT ST_ConcaveHull(pts, 0.5, false) FROM word_pts;
``````

A little better, but still hard to recognize 'postgis'. What if we allowed holes?

``````SELECT ST_ConcaveHull(pts, 0.5, true) FROM word_pts;
``````

This starts to look a bit more like the word 'postgis', with the hole in 'p' being clear.

As we start to make the shape more concave, it begins to take on more and more recognizable as 'postgis'....until it doesn't and starts to look closer to modern art.

``````SELECT ST_ConcaveHull(pts, 0.35, true) FROM word_pts;
``````

``````SELECT ST_ConcaveHull(pts, 0.05, true) FROM word_pts;
``````

## Polygons too!

`ST_ConcaveHull` is also useful on multipolygons, and follows the same properties as demo'd on multipoints. It's important to note that if there are already holes in the existing multipolygon, setting `param_allow_holes=false` will still create convex polygons with "holes" in the middle, following the original polygon. The concave hulls will always contains the original polygons!

``````SELECT ST_ConcaveHull(ST_Letters('postgis'), 0.5, false);
``````

As the convexity decreases and holes are allowed, the shape looks more and more like the original polygons in the original table.

``````SELECT ST_ConcaveHull(ST_Letters('postgis'), 0.1, true);
``````

## `ST_TriangulatePolygon`

The last demo here is the function `ST_TriangulatePolygon`, new in PostGIS 3.3. This function computes the "best quality" triangulation of a polygon (and also works on multipolygons too!). This can be extremely useful for computing meshes of polygons in a quick and efficient manner.

``````SELECT ST_TriangulatePolygon(ST_Letters('postgis'));
``````

## Summary

`ST_Letters` provides a useful starting point for demoing functions on points and polygons. The new improvements in `ST_ConcaveHull` make it more useful for generating concave hulls of geometries and they are significantly more intuitive to use. `ST_TriangulatePolygon` can be useful for finding meshes of polygons and multipolygons. The team at Crunchy Data will continue to make important contributions to PostGIS in order to help our users create interesting and innovative open source solutions!

## January 06, 2023

### Crunchy Data

#### Timezone Transformation Using Location Data & PostGIS

Imagine your system captures event data from all over the world but the data all comes in UTC time giving no information about the local timing of an event. How can we quickly convert between the UTC timestamp and local time zone using GPS location? We can quickly solve this problem using PostgreSQL and PostGIS.

This example assumes you have a Postgres database running with PostGIS. If you’re new to PostGIS, see PostGIS for Newbies.

### Steps we will follow

1. Timezone Shape file Overview: For World Timezone shape file, I have been following a really nice project by Evan Siroky, Timezone Boundary Builder. We’ll download the timezones-with-oceans.shapefile.zip from this location.
2. Load Shape file: Using shp2pgsql, convert shape file to sql to create timezones_with_ocean table.
3. PostgreSQL internal view pg_timezone_names: Understand pg_timezone_names view.
4. Events table and insert sample data: Create events table and insert sample data.
5. Transformation Query: Transform event UTC timestamp to event local timestamp.

### Overview of data relationship

Below is an overview of the data relationship and join conditions we will be using.

### Timezone Shape file Overview

A “shape file” commonly refers to a collection of files with .shp, .shx, .dbf, and other extensions on a common prefix name, which in our case is combined-shapefile-with-oceans.*. **combined-shapefile-with-oceans contains polygons with the boundaries of the world's timezones. With this data we can start our process.

We will be using shp2pgsql to generate sql file from shape file to create public.timezones_with_ocean and inserts data in a table. The table contains fields gid, tzid and geometry.

Export the Host, user, password variables

``````export PGHOST=p.<pgcluster name>.db.postgresbridge.com
export PGUSER=<DBUser>
``````

Create sql file from shape file

``````shp2pgsql -s 4326  "combined-shapefile.shp" public.timezones_with_oceans   > timezone_shape.sql
``````

Create public.timezones_with_ocean and load timestamp data

``````psql -d timestamp -f timezone_shape.sql
``````

Query a bit of sample data

``````SELECT tzid, ST_AsText(geom), geom FROM public.timezone_with_oceans limit 10;
``````

Visualize Sample data

Using PgAdmin highlight geom column and click on eye icon visualize the geometry on map showing below.

### PostgreSQL internal view pg_timezone_names

PostgreSQL provides a view of pg_timezone_names with a list of time zone names recognized by SET TIMEZONE. By default, PostgreSQL also provides their associated abbreviations, UTC offsets, and daylight-savings status, which our clients need to know.

pg_timezone_names view columns description

Column Type Description
name text Time zone name
abbrev text Time zone abbreviation
utc_offset interval Offset from UTC (positive means east of Greenwich)
is_dst bool True if currently observing daylight savings

pg_timezone sample data

### Events table and insert sample data

Now that we have the timezone shape file loaded, we can create an event table, load sample transaction data, and apply a timestamp conversion transformation query.

``````CREATE TABLE IF NOT EXISTS public.events
(
event_id bigint NOT NULL,
eventdatetime timestamp without time zone NOT NULL,
event_type varchar(25) not null,
latitude double precision NOT NULL,
longitude double precision NOT NULL,
CONSTRAINT events_pkey PRIMARY KEY (event_id)
);

INSERT INTO public.events(
event_id, eventdatetime, event_type, latitude, longitude)
VALUES (10086492,'2021-08-17 23:17:05','Walking',34.894089,-86.51148),
(50939,'2021-08-19 10:27:12','Hiking',34.894087,-86.511484),
(10086521,'2021-09-09 19:32:37','Swiming',34.642584,-86.761291),
(22465493,'2021-09-30 11:43:34','Swiming',33.611151,-86.799522),
(22465542,'2021-11-26 22:40:44.197','Swiming',34.64259,-86.761452),
(22465494,'2021-09-30 11:43:34','Hiking',33.611151,-86.799522),
(10087348,'2021-07-01 13:42:15','Swiming',25.956098,-97.535303),
(22466679,'2021-09-01 12:25:06','Hiking',25.956112,-97.535304),
(22466685,'2021-09-02 13:41:07','Swiming',25.956102,-97.535305),
(10088223,'2021-11-29 13:19:53','Hiking',25.956097,-97.535303),
(22246192,'2021-06-16 22:21:23','Walking',37.083726,-113.577984),
(9844188,'2021-06-23 20:18:43','Swiming',37.1067,-113.561401),
(22246294,'2021-06-25 21:50:06','Walking',37.118719,-113.598038),
(22246390,'2021-07-01 18:15:54','Hiking',37.109579,-113.562923),
(9844332,'2021-07-04 19:11:13','Walking',37.251538,-113.614708),
(9845242,'2021-11-04 13:25:40.425','Swiming',37.251542,-113.614699),
(84843,'2021-11-23 14:33:20','Swiming',37.251541,-113.614698),
(22247674,'2021-12-21 14:31:15','Swiming',37.251545,-113.614691),
(22246714,'2021-08-09 14:46:51','Swiming',37.109597,-113.562912),
(9845116,'2021-10-18 14:59:51','Swiming',37.082777,-113.554991);
``````

Sample Event Data

### Transformation Query

Now we can convert the UTC timestamp to the local time for an event. Using PostGIS function St_Intersects, we can find the timezone_with_oceans.geom polygon in which an event point lies. This gives the name of the timezone where the event occurred. To create our transformation query:

• First we create the location geometry using Longitude and Latitude from the events table.

• Using PostGIS function St_Intersects, we will find common points between timezone_with_oceans.geom and an event’s location geometry giving us information on where the event occurred.

• Join pg_timezone_names to timezone_with_oceans on name and tzid respectively, to retrieve abbrev, utc_offset, and is_dst fields from pg_timezone_names.

• Using PostgreSQL AT TIME ZONE operator and pg_timezone_name, we convert UTC event timestamp to local event timestamp completing the process, e.g.

`timestamp '2021-07-05 00:59:12' at time zone 'America/Denver' → 2021-07-04 18:59:12+00’`

Transformation Query SQL:

``````SELECT event_id, latitude, longitude, abbrev,
utc_offset,is_dst, eventdatetime,
((eventdatetime::timestamp WITH TIME ZONE AT TIME ZONE abbrev)::timestamp WITH TIME ZONE)
AS eventdatetime_local
FROM public.events
JOIN timezone_with_oceans ON ST_Intersects(ST_Point(longitude, latitude, 4326) , geom)
JOIN pg_timezone_names ON tzid = name;
``````

## Closing Thoughts

PostgreSQL and PostGIS allow you to easily and dynamically solve timezone transformation. I hope this blog was helpful, and we at Crunchy Data wish you happy learning.

## December 04, 2022

### Boston GIS (Regina Obe, Leo Hsu)

#### PostGIS Day 2022 Videos are out and some more highlights

Elizabeth Chistensen already gave a succinct summary of some of the PostGIS Day 2022 presentations, which you can see here.

There were many case study presentations which involved use of PostGIS, QGIS, OpenStreetMap, and pgRouting (and other extensions that extend PostGIS) as well as many "How to" videos. There were also talks on "How PostGIS is made". I'll highlight some of these, which overlap with Elizabeth's list but different angle of view.

Continue reading "PostGIS Day 2022 Videos are out and some more highlights"

## December 02, 2022

### Crunchy Data

#### PostGIS Day 2022

Crunchy Data hosted the 4th annual PostGIS Day on November 17, 2022. PostGIS Day always comes a day after GIS Day which occurs annually on the 3rd Wednesday of November.

We had speakers from 10 different countries and attendees from more than 70 countries.

PostGIS is the most popular spatial relational database worldwide with:

• An extensive catalog of spatial functions
• Rich ecosystem of in-db extensions for routing, event management, external database linkage, point clouds, rasters and more
• Numerous web frameworks for serving tiles and vector images
• Easy integration with QGIS, the most popular open-source mapping software

## PostGIS is in our Everyday Lives

What blew me away this year is how PostGIS touches so many parts of our everyday lives. From high speed internet, to water and flood management, there’s PostGIS supporting major parts of the infrastructure we depend on every day.

### Water, Soil, Floods, and Environment

We heard a talk about the GISwater project, which is open source asset management for water and wastewater tracking. Helping humanity, bringing the world water resources, using the open source tools. What could be better? This project just gives everyone all the feels.

In a similar way, my ears really perked up with the Vizzuality talk on using PostGIS for conversation efforts, tracking food supply chains, and sustainability.

We were really lucky to have two federal agencies join us as well. USDA talked about using PostGIS as the engine behind the Dynamic Soils Hub project, which tracked a number of key soil changes like erosion and farmland loss across America. A developer who works alongside the National Flood Insurance Program talked about some of his work to collect several massive open data sets and merge them into a set of building outlines for the entire country.

I also really appreciated several talks that went a little bit outside the box to talk about how PostGIS can be used in a more people focused context.  We had a nice talk about using PostGIS/QGIS to look at how you can explore infrastructure and how it impacts social inequity and health outcomes. There was also a great talk about routing people's movements outside roads, like on foot or inside buildings. Another talk went outside the box, or should I say inside the goal box, with a talk on how geospatial data is being collected and used in hockey and sports.

### Developing nations

We were really fortunate to get a couple talks on how developers in Africa are using PostGIS. We had a talk on using QField, QGIS, and PostGIS for urban planning in Nigeria. A developer from Tanzania showed an isochrone map plugin he wrote for QGIS. Isochrone maps show distance from a certain point with time considered.

We had two talks from Canada, where they’re aiming to have the entire country covered by high speed internet by 2030. One of the talks described building a design system in QGIS and PostGIS and even training internal designers to use these tools. We had a second talk about the telecom industry, this one more focused on getting things up to scale for backing a really large enterprise.

## Tooling Showcases

PosGIS Day also had a really broad set of talks about ecosystem tools. Seeing these all in total, you realize how big the world of PostGIS is and how many people are building large scale applications out of geospatial data.

### PostGIS at scale

• pgSTAC, an asset catalog used by the Microsoft Planetary Computer and NASA projects
• pg_eventserv, an extension for doing real time updates and web sockets

### Ecosystem projects

Thanks to everyone who participated this year! If you missed it, plan to join us next year. Mark your calendar, we’ll open a call for papers in September of 2023. Videos are linked above or available as a full playlist on Crunchy Data’s YouTube channel.

### PostGIS at Crunchy Data

PostGIS is a supported extension with Crunchy Certified Postgres and is packaged with our products on traditional infrastructure, Kubernetes, and our fully managed Crunchy Bridge. Crunchy Data supports a variety of enterprise clients using PostGIS in production. If you’re using PostGIS and are looking for a host or supported distribution, contact us, we’d love to chat.

## November 19, 2022

### Boston GIS (Regina Obe, Leo Hsu)

#### Post GIS Day 2022 Celebrations

I didn't think last two years Post GIS Day conferences could be topped, but I was wrong. This year's was absolutely fabulous and better than all the others. We had two key conferences going on this year (yesterday). The first was the Cruncy Data PostGIS Day 2022 12 hour marathon of nothing but PostGIS related talks. Many thanks to Elizabeth Christensen and Paul Ramsey for putting it all together and for 12 hrs. The PostGIS day celebrations climaxed with many of us playing capture the flag on Paul's new shiny invention until we managed to crash it.

Continue reading "Post GIS Day 2022 Celebrations"

## November 13, 2022

### PostGIS Development

#### PostGIS 3.3.2, 3.2.4, 3.1.8, 3.0.8 Patch Releases

The PostGIS development team is pleased to provide security, bug fixes and performance enhancements 3.3.2, 3.2.4, 3.1.8 and 3.0.8 for the 3.3, 3.2, 3.1, and 3.0 stable branches.

### PostGIS Development

#### PostGIS 3.3.2, 3.2.4, 3.1.8, 3.0.8 Patch Releases

The PostGIS development team is pleased to provide security, bug fixes and performance enhancements 3.3.2, 3.2.4, 3.1.8 and 3.0.8 for the 3.3, 3.2, 3.1, and 3.0 stable branches.

## November 12, 2022

### PostGIS Development

#### PostGIS 2.5.9 EOL

The PostGIS Team is pleased to release PostGIS 2.5.9! This is the end-of-life release of the 2.5 branch.

This release is a bug fix release, addressing issues found in the previous 2.5 releases.

### PostGIS Development

#### PostGIS 2.5.9 EOL

The PostGIS Team is pleased to release PostGIS 2.5.9! This is the end-of-life release of the 2.5 branch.

This release is a bug fix release, addressing issues found in the previous 2.5 releases.

## October 24, 2022

### Crunchy Data

#### Moving Objects and Geofencing with Postgres & PostGIS

In a recent post, we introduced pg_eventserv and the real-time web notifications from database actions.

In this post, we will dive into a practical use case: displaying state, calculating events, and tracking historical location for a set of moving objects.

This demonstration uses pg_eventserv for eventing, and pg_featureserv for external web API, and OpenLayers as the map API, to build a small example application that shows off the common features of moving objects systems.

Try it out!

## Features

Moving objects applications can be very complex or very simple, but they usually include a few common baseline features:

• A real-time view of the state of the objects.
• Live notifications when objects enter and leave a set of "geofences".
• Querying the history of the system, to see where objects have been, and to summarize their state (eg, "truck 513 spent 50% of its time in the yard").

## Tables

The model has three tables:

• objects where the current location of the objects is stored.
• geofences where the geofences are stored.
• objects_history where the complete set of all object locations are stored.

## Architecture (in brief)

From the outside, the system has the following architecture:

Changes to objects are communicated in via a web API backed by pg_featureserv, those changes fires a bunch of triggers that generate events that pg_eventserv pushes out to listening clients via WebSockets.

## Architecture (in detail)

• The user interface generates object movements, via the arrow buttons for each object. This is in lieu of a real "moving object" fleet in the real world generating timestamped GPS tracks.

• Every movement click on the UI fires a call to a web API, which is just a function published via pg_featureserv, `object_move(object_id, direction)`.

postgisftw.object_move(object_id, direction)
``````CREATE OR REPLACE FUNCTION postgisftw.object_move(
move_id integer, direction text)
RETURNS TABLE(id integer, geog geography)
AS \$\$
DECLARE
xoff real = 0.0;
yoff real = 0.0;
step real = 2.0;
BEGIN

yoff := CASE
WHEN direction = 'up' THEN 1 * step
WHEN direction = 'down' THEN -1 * step
ELSE 0.0 END;

xoff := CASE
WHEN direction = 'left' THEN -1 * step
WHEN direction = 'right' THEN 1 * step
ELSE 0.0 END;

RETURN QUERY UPDATE moving.objects mo
SET geog = ST_Translate(mo.geog::geometry, xoff, yoff)::geography,
ts = now()
WHERE mo.id = move_id
RETURNING mo.id, mo.geog;

END;
\$\$
LANGUAGE 'plpgsql' VOLATILE;
``````
• The `object_move(object_id, direction)` function just converts the "direction" parameter into a movement vector, and `UPDATES` the relevant row of the `objects` table.

• The change to the `objects` table fires off the `objects_geofence()` trigger, which calculates the fences the object is now in.

objects_geofence()
``````CREATE FUNCTION objects_geofence() RETURNS trigger AS \$\$
DECLARE
fences_new integer[];
BEGIN
-- Add the current geofence state to the input
-- tuple every time.
SELECT coalesce(array_agg(id), ARRAY[]::integer[])
INTO fences_new
FROM moving.geofences
WHERE ST_Intersects(geofences.geog, new.geog);

RAISE DEBUG 'fences_new %', fences_new;
-- Ensure geofence state gets saved
NEW.fences := fences_new;
RETURN NEW;
END;

\$\$
LANGUAGE 'plpgsql';
``````
• The change to the `objects` table then fires off the `objects_update()` trigger, which:
• Compares the current set of geofences to the previous set, and thus detects any enter/leave events.
• Adds the new location of the object to the `objects_history` tracking table.
• Composes the new location and any geofence events into a JSON object and puts it into the "objects" `NOTIFY` queue using `pg_notify()`.
objects_update()
``````CREATE FUNCTION objects_update() RETURNS trigger AS
\$\$

DECLARE
channel text := 'objects';
fences_old integer[];
fences_entered integer[];
fences_left integer[];
events_json jsonb;
location_json jsonb;
BEGIN
-- Place a copy of the value into the history table
INSERT INTO moving.objects_history (id, geog, ts, props)
VALUES (NEW.id, NEW.geog, NEW.ts, NEW.props);

-- Clean up any nulls
fences_old := coalesce(OLD.fences, ARRAY[]::integer[]);
RAISE DEBUG 'fences_old %', fences_old;

-- Compare to previous fences state
fences_entered = NEW.fences - fences_old;
fences_left = fences_old - NEW.fences;

RAISE DEBUG 'fences_entered %', fences_entered;
RAISE DEBUG 'fences_left %', fences_left;

-- Form geofence events into JSON for notify payload
WITH r AS (
SELECT 'entered' AS action,
g.id AS geofence_id,
g.label AS geofence_label
FROM moving.geofences g
WHERE g.id = ANY(fences_entered)
UNION
SELECT 'left' AS action,
g.id AS geofence_id,
g.label AS geofence_label
FROM moving.geofences g
WHERE g.id = ANY(fences_left)
)
SELECT json_agg(row_to_json(r))
INTO events_json
FROM r;

SELECT json_build_object(
'type', 'objectchange',
'object_id', NEW.id,
'events', events_json,
'location', json_build_object(
'longitude', ST_X(NEW.geog::geometry),
'latitude', ST_Y(NEW.geog::geometry)),
'ts', NEW.ts,
'color', NEW.color,
'props', NEW.props)

-- Send the payload out on the channel
PERFORM (
);

RETURN NEW;
END;

\$\$
LANGUAGE 'plpgsql';
``````
• pg_eventserv picks the event off the `NOTIFY` queue and pushes it out to all listening clients over WebSockets.
• The user interface recieves the JSON payload, parses it, and applies the new location to the appropriate object. If there is a enter/leave event on a geofence, the UI also changes the geofence outline color appropriately.

Phew! That's a lot!

• Side note, the `geofences` table also has a trigger, `layer_change()` that catches insert/update/delete events and publishes a JSON notification with `pg_notify()`. This is also published by pg_eventserv and when the UI receives it, it simply forces a full re-load of geofence data.
layer_change()
``````CREATE FUNCTION layer_change() RETURNS trigger AS
\$\$

DECLARE
layer_change_json json;
channel text := 'objects';
BEGIN
-- Tell the client what layer changed and how
SELECT json_build_object(
'type', 'layerchange',
'layer', TG_TABLE_NAME::text,
'change', TG_OP)
INTO layer_change_json;

RAISE DEBUG 'layer_change %', layer_change_json;
PERFORM (
SELECT pg_notify(channel, layer_change_json::text)
);
RETURN NEW;
END;

\$\$
LANGUAGE 'plpgsql';
``````

OK, all done.

## Trying It Out Yourself

All the code and instructions are available in the moving objects example of `pg_eventserv`.

## Conclusion

• Moving objects are a classic case of "system state stored in the database".
• PostgreSQL provides the LISTEN/NOTIFY system to update clients about real-time changes.
• The pg_eventserv service allows you to push LISTEN/NOTIFY events further out to any WebSockets client and generate a moving object map.
• Because the state is managed in the database, storing the historical state of the system is trivially easy. \$\$

## October 17, 2022

### Boston GIS (Regina Obe, Leo Hsu)

#### Using pg_upgrade to upgrade PostgreSQL 9.6 PostGIS 2.4 to PostgreSQL 15 3.3 on Yum

PostgreSQL 15 came out just last week. To celebrate the arrival of PostgreSQL 15, I will revisit the number one problem people have with PostGIS, how to upgrade your PostGIS enabled cluster, without installing an old PostGIS version.

In a previous article Using pg upgrade to upgrade PostGIS without installing older version I demonstrated a trick for upgrading to a newer PostgreSQL instance from PostGIS 2.2 - 2.whatever without having to install the older version of PostGIS in your new PostgreSQL service. This is a revisit of that article, but with considerations for upgrading from a 2 series to a 3 series.

Fear not. Going from PostGIS 2+ to PostGIS 3+ can still be done without installing the old PostGIS 2+ in your new cluster. Additionally once you are on PostGIS 3.1+, you should never have to do symlink or copy hacks to upgrade say PostGIS 3.1 to PostGIS 3.4. At least not for the same major version.

If per chance you are on PostGIS 1 and trying to jump all the way to PostGIS 3+, then sadly you need to do a pg_dump of your old database/cluster and pg_restore into your new cluster.

## October 13, 2022

### Martin Davis (Lin.ear th.inking)

#### Relational Properties of DE-9IM spatial predicates

There is an elegant mathematical theory of binary relations.  Homogeneous relations are an important subclass of binary relations in which both domains are the same.  A homogeneous relation R is a subset of all ordered pairs (x,y) with x and y elements of the domain.  This can be thought of as a boolean-valued function R(x,y), which is true if the pair has the relationship and false if not.

The restricted structure of homogeneous relations allows describing them by various properties, including:

Reflexive - R(x, x) is always true

Irreflexive - R(x, x) is never true

Symmetric - if R(x, y), then R(y, x)

Antisymmetric - if R(x, y) and R(y, x), then x = y

Transitive - if R(x, y) and R(y, z), then R(x, z)

The Dimensionally Extended 9-Intersection Model (DE-9IM) represents the topological relationship between two geometries.

Various useful subsets of spatial relationships are specified by named spatial predicates.  These are the basis for spatial querying in many spatial systems including the JTS Topology Suite, GEOS and PostGIS.

The spatial predicates are homogeneous binary relations over the Geometry domain  They can thus be categorized in terms of the relational properties above.  The following table shows the properties of the standard predicates:

#### Notes

• Contains and Within are not Transitive because of the quirk that "Polygons do not contain their Boundary" (explained in this post).  A counterexample is a Polygon that contains a LineString lying in its boundary and interior, with the LineString containing a Point that lies in the Polygon boundary.  The Polygon does not contain the Point.  So Contains(poly, line) = true and Contains(line, pt) = true, but Contains(poly, pt) = false.   The predicates Covers and CoveredBy do not have this idiosyncrasy, and thus are transitive.

## September 30, 2022

### Anita Graser (Underdark)

#### Detecting close encounters using MobilityDB 1.0

It’s been a while since we last talked about MobilityDB in 2019 and 2020. Since then, the project has come a long way. It joined OSGeo as a community project and formed a first PSC, including the project founders Mahmoud Sakr and Esteban Zimányi as well as Vicky Vergara (of pgRouting fame) and yours truly.

This post is a quick teaser tutorial from zero to computing closest points of approach (CPAs) between trajectories using MobilityDB.

## Setting up MobilityDB with Docker

The easiest way to get started with MobilityDB is to use the ready-made Docker container provided by the project. I’m using Docker and WSL (Windows Subsystem Linux on Windows 10) here. Installing WLS/Docker is out of scope of this post. Please refer to the official documentation for your operating system.

Once Docker is ready, we can pull the official container and fire it up:

```docker pull mobilitydb/mobilitydb
docker volume create mobilitydb_data
docker run --name "mobilitydb" -d -p 25432:5432 -v mobilitydb_data:/var/lib/postgresql mobilitydb/mobilitydb
psql -h localhost -p 25432 -d mobilitydb -U docker
```

Currently, the container provides PostGIS 3.2 and MobilityDB 1.0:

Once the container is running, we can already connect to it from QGIS. This is my preferred way to load data into MobilityDB because we can simply drag-and-drop any timestamped point layer into the database:

For this post, I’m using an AIS data sample in the region of Gothenburg, Sweden.

After loading this data into a new table called ais, it is necessary to remove duplicate and convert timestamps:

```CREATE TABLE AISInputFiltered AS
SELECT DISTINCT ON("MMSI","Timestamp") *
FROM ais;

ALTER TABLE AISInputFiltered ADD COLUMN t timestamp;
UPDATE AISInputFiltered SET t = "Timestamp"::timestamp;
```

Afterwards, we can create the MobilityDB trajectories:

```CREATE TABLE Ships AS
SELECT "MMSI" mmsi,
tgeompoint_seq(array_agg(tgeompoint_inst(Geom, t) ORDER BY t)) AS Trip,
tfloat_seq(array_agg(tfloat_inst("SOG", t) ORDER BY t) FILTER (WHERE "SOG" IS NOT NULL) ) AS SOG,
tfloat_seq(array_agg(tfloat_inst("COG", t) ORDER BY t) FILTER (WHERE "COG" IS NOT NULL) ) AS COG
FROM AISInputFiltered
GROUP BY "MMSI";

ALTER TABLE Ships ADD COLUMN Traj geometry;
UPDATE Ships SET Traj = trajectory(Trip);
```

Once this is done, we can load the resulting Ships layer and the trajectories will be loaded as lines:

## Computing closest points of approach

To compute the closest point of approach between two moving objects, MobilityDB provides a shortestLine function. To be correct, this function computes the line connecting the nearest approach point between the two tgeompoint_seq. In addition, we can use the time-weighted average function twavg to compute representative average movement speeds and eliminate stationary or very slowly moving objects:

```SELECT S1.MMSI mmsi1, S2.MMSI mmsi2,
shortestLine(S1.trip, S2.trip) Approach,
ST_Length(shortestLine(S1.trip, S2.trip)) distance
FROM Ships S1, Ships S2
WHERE S1.MMSI > S2.MMSI AND
twavg(S1.SOG) > 1 AND twavg(S2.SOG) > 1 AND
dwithin(S1.trip, S2.trip, 0.003)
```

In the QGIS Browser panel, we can right-click the MobilityDB connection to bring up an SQL input using Execute SQL:

The resulting query layer shows where moving objects get close to each other:

To better see what’s going on, we’ll look at individual CPAs:

## Having a closer look with the Temporal Controller

Since our filtered AIS layer has proper timestamps, we can animate it using the Temporal Controller. This enables us to replay the movement and see what was going on in a certain time frame.

I let the animation run and stopped it once I spotted a close encounter. Looking at the AIS points and the shortest line, we can see that MobilityDB computed the CPAs along the trajectories:

A more targeted way to investigate a specific CPA is to use the Temporal Controllers’ fixed temporal range mode to jump to a specific time frame. This is helpful if we already know the time frame we are interested in. For the CPA use case, this means that we can look up the timestamp of a nearby AIS position and set up the Temporal Controller accordingly:

## More

I hope you enjoyed this quick dive into MobilityDB. For more details, including talks by the project founders, check out the project website.

This post is part of a series. Read more about movement data in GIS.