“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 😊.

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:

1. Dataset import
2. Layer setup
3. Setup areas around airports
4. 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] ); 

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) 

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

Figure 1 Major cities and international airports

Figure 2 Major cities, towns and international airports

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 

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 ; 

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);