"Catchment areas" with PostgreSQL and PostGIS

Recently a colleague in our sales department asked me for a way to partition an area of interest spatially. He wanted to approximate customer potential and optimize our sales strategies respective trips.
Furthermore he wanted the resulting regions to be created around international airports first, and then intersected by potential customer locations, in order to support a basic ranking. Well, I thought this would be a good opportunity to experiment with lesser-known PostGIS functions ????.

Table of Contents

In our case, customer locations could be proxied by major cities.
It's also worth mentioning that I only used freely available datasets for this demo.

Here is the basic structure of the process:

Dataset import
Layer setup
Setup areas around airports
Intersect areas with major cities

Dataset import

Our sales team will kick off their trips from international airports in Germany, so we require locations and classifications of airports from Germany. Airports can be extracted from OpenStreetMap datasets filtering osm_points by Tag:aeroway=aerodrome. Alternatively, any other data source containing the most recent airport locations and classifications is welcome. For our demo, I decided to give a freely available csv from “OurAirports” a chance ????.

In addition, I used OpenStreetMap datasets to gather both spatial and attributive data for cities located in Germany. To do so, I downloaded the most recent dataset covering Germany from Geofabrik and subsequently triggered the PostGIS import utilizing osm2pgsql. Please checkout my blogpost on this topic to get more details on how to accomplish this easily.

At the end of this article, you will find data structures for airports, osm_points and osm_polygons (country border). These are the foundation for our analysis.

Layer setup

Let’s start with cities and filter osm_points by place to retrieve only the major cities and towns.

SELECT planet_osm_point.name,
       planet_osm_point.way,
       planet_osm_point.tags -> 'population' :: text AS population,
       planet_osm_point.place
FROM   osmtile_germany.planet_osm_point
WHERE  planet_osm_point.place = ANY ( array [ 'city' :: text, 'town' :: text] );

SELECT planet_osm_point.name,

planet_osm_point.way,

planet_osm_point.tags -> 'population' :: text AS population,

planet_osm_point.place

FROM osmtile_germany.planet_osm_point

WHERE planet_osm_point.place = ANY ( array [ 'city' :: text, 'town' :: text] );

Since we are only interested in international airports within Germany, let’s get rid of those tiny airports and heliports which are not relevant for our simple analysis.

WITH border AS
(
       SELECT way
       FROM osmtile_germany.planet_osm_polygon
       WHERE admin_level = '2'
       AND boundary = 'administrative')
SELECT geom AS geom,
       NAME,
       type
FROM   osmtile_germany.airports,
       border
WHERE  type = 'large_airport'
AND    St_intersects(St_transform(airports.geom, 3857), border.way)

WITH border AS

(

SELECT way

FROM osmtile_germany.planet_osm_polygon

WHERE admin_level = '2'

AND boundary = 'administrative')

SELECT geom AS geom,

NAME,

type

FROM osmtile_germany.airports,

border

WHERE type = 'large_airport'

AND St_intersects(St_transform(airports.geom, 3857), border.way)

Check out Figures 1 and 2 below to get an impression of what our core layers look like.

Areas around airports

So far, so good. It is time to start with our analysis by generating catchment areas around international airports first.
Instead of creating Voronoi polygons (Voronoi-Diagrams) using our preferred GIS-client, let us instead investigate how we can accomplish this task using only PostGIS’ functions.

WITH border AS
(
       SELECT way
       FROM osmtile_germany.planet_osm_polygon
       WHERE admin_level = '2'
       AND boundary = 'administrative'), 
airports AS
(
       SELECT St_collect(St_transform(geom, 3857)) AS geom
       FROM osmtile_germany.airports,
       border
       WHERE type = 'large_airport'
       AND St_intersects(St_transform(airports.geom, 3857), border.way))
SELECT (St_dump(St_voronoipolygons(geom))).geom
FROM   airports

WITH border AS

(

SELECT way

FROM osmtile_germany.planet_osm_polygon

WHERE admin_level = '2'

AND boundary = 'administrative'),

airports AS

(

SELECT St_collect(St_transform(geom, 3857)) AS geom

FROM osmtile_germany.airports,

border

WHERE type = 'large_airport'

AND St_intersects(St_transform(airports.geom, 3857), border.way))

SELECT (St_dump(St_voronoipolygons(geom))).geom

FROM airports

Intersection of catchment areas with potential customers

So how do we proceed? Let’s count the major cities by area as a proxy for potential customer counts, in order to rank the most interesting areas. As a reminder - major cities are only used here as a proxy for more detailed customer locations. They should be enriched by further datasets in order to more realistically approximate customer potential. Below, you see the final query which outputs major city count by polygon. Figure 4 represents the resulting Voronoi polygons extended by their major city counts. Query results and visualisation might be used as building blocks for optimizing customer acquisition in the future.

WITH border AS
(
       SELECT way
       FROM osmtile_germany.planet_osm_polygon
       WHERE admin_level = '2'
       AND boundary = 'administrative'), 
voronoi_polys AS
(
       SELECT (St_dump(St_voronoipolygons(St_collect(St_transform(geom, 3857))))).geom geom
       FROM osmtile_germany.airports,
            border
       WHERE type = 'large_airport'
       AND St_intersects(St_transform(airports.geom, 3857), border.way))
SELECT voronoi_polys.geom,
       a.NAME,
       b.majorcitycount
FROM   voronoi_polys,
       lateral
       (
              SELECT NAME
              FROM   osmtile_germany.airports,
                     border
              WHERE St_intersects(voronoi_polys.geom, St_transform(osmtile_germany.airports.geom, 3857))
              AND St_intersects(border.way, St_transform(osmtile_germany.airports.geom, 3857))
              AND osmtile_germany.airports.type = 'large_airport' ) AS a,
       lateral
       (
              SELECT count(*) majorcitycount
              FROM osmtile_germany.planet_osm_point
              WHERE place IN ('city')
              AND St_intersects(voronoi_polys.geom, osmtile_germany.planet_osm_point.way) ) AS b
ORDER BY 3 DESC ;

WITH border AS

(

SELECT way

FROM osmtile_germany.planet_osm_polygon

WHERE admin_level = '2'

AND boundary = 'administrative'),

voronoi_polys AS

(

SELECT (St_dump(St_voronoipolygons(St_collect(St_transform(geom, 3857))))).geom geom

FROM osmtile_germany.airports,

border

WHERE type = 'large_airport'

AND St_intersects(St_transform(airports.geom, 3857), border.way))

SELECT voronoi_polys.geom,

a.NAME,

b.majorcitycount

FROM voronoi_polys,

lateral

(

SELECT NAME

FROM osmtile_germany.airports,

border

WHERE St_intersects(voronoi_polys.geom, St_transform(osmtile_germany.airports.geom, 3857))

AND St_intersects(border.way, St_transform(osmtile_germany.airports.geom, 3857))

AND osmtile_germany.airports.type = 'large_airport' ) AS a,

lateral

(

SELECT count(*) majorcitycount

FROM osmtile_germany.planet_osm_point

WHERE place IN ('city')

AND St_intersects(voronoi_polys.geom, osmtile_germany.planet_osm_point.way) ) AS b

ORDER BY 3 DESC ;

Final thoughts and outlook

Today we started with some basic spatial analytics and took the opportunity to experiment with less known PostGIS functions.
Hopefully you enjoyed this session, and even more importantly, you’re now eager to play around with PostGIS and its rich feature set. It's just waiting to tackle your real-world challenges.

Annex

CREATE TABLE IF NOT EXISTS osmtile_germany.planet_osm_point (
    osm_id bigint,
    name text,
    place text,
    tags hstore,
    way geometry(point, 3857)
);

CREATE INDEX IF NOT EXISTS planet_osm_point_way_idx ON osmtile_germany.planet_osm_point USING gist (way);

CREATE INDEX IF NOT EXISTS idx_osm_point_place ON osmtile_germany.planet_osm_point (place);

CREATE TABLE IF NOT EXISTS osmtile_germany.airports (
    gid serial NOT NULL CONSTRAINT airports_pkey PRIMARY KEY,
    type varchar(254),
    name varchar(254),
    geom geometry(point, 4326)
);

CREATE INDEX IF NOT EXISTS airports_geom_idx ON osmtile_germany.airports USING gist (geom);

CREATE INDEX IF NOT EXISTS airports_geom_3857_idx ON osmtile_germany.airports USING gist (St_transform (geom, 3857));

CREATE INDEX IF NOT EXISTS airports_type_index ON osmtile_germany.airports (type);

CREATE TABLE IF NOT EXISTS osmtile_germany.planet_osm_polygon (
    osm_id bigint,
    admin_level text,
    boundary text,
    name text,
    way geometry(Geometry, 3857)
);

CREATE INDEX IF NOT EXISTS planet_osm_polygon_way_idx ON osmtile_germany.planet_osm_polygon USING gist (way);

CREATE INDEX IF NOT EXISTS planet_osm_polygon_osm_id_idx ON osmtile_germany.planet_osm_polygon (osm_id);

CREATE INDEX IF NOT EXISTS planet_osm_polygon_admin_level_boundary_index ON osmtile_germany.planet_osm_polygon (admin_level, boundary);

CREATE TABLE IF NOT EXISTS osmtile_germany.planet_osm_point (

osm_id bigint,

name text,

place text,

tags hstore,

way geometry(point, 3857)

);

CREATE INDEX IF NOT EXISTS planet_osm_point_way_idx ON osmtile_germany.planet_osm_point USING gist (way);

CREATE INDEX IF NOT EXISTS idx_osm_point_place ON osmtile_germany.planet_osm_point (place);

CREATE TABLE IF NOT EXISTS osmtile_germany.airports (

gid serial NOT NULL CONSTRAINT airports_pkey PRIMARY KEY,

type varchar(254),

name varchar(254),

geom geometry(point, 4326)

);

CREATE INDEX IF NOT EXISTS airports_geom_idx ON osmtile_germany.airports USING gist (geom);

CREATE INDEX IF NOT EXISTS airports_geom_3857_idx ON osmtile_germany.airports USING gist (St_transform (geom, 3857));

CREATE INDEX IF NOT EXISTS airports_type_index ON osmtile_germany.airports (type);

CREATE TABLE IF NOT EXISTS osmtile_germany.planet_osm_polygon (

osm_id bigint,

admin_level text,

boundary text,

name text,

way geometry(Geometry, 3857)

);

CREATE INDEX IF NOT EXISTS planet_osm_polygon_way_idx ON osmtile_germany.planet_osm_polygon USING gist (way);

CREATE INDEX IF NOT EXISTS planet_osm_polygon_osm_id_idx ON osmtile_germany.planet_osm_polygon (osm_id);

CREATE INDEX IF NOT EXISTS planet_osm_polygon_admin_level_boundary_index ON osmtile_germany.planet_osm_polygon (admin_level, boundary);

"Catchment areas" with PostgreSQL and PostGIS

Dataset import

Layer setup

Areas around airports

Intersection of catchment areas with potential customers

Final thoughts and outlook

Annex

Leave a Reply Cancel reply

Florian Nadler

Blog Tags

NEWSLETTER

Articles by our PostgreSQL Experts