A sales colleague needed a way to partition an area of interest spatially, in order to approximate customer potential and prioritize sales trips. The regions were to be built around international airports and then intersected with potential customer locations to support a basic ranking. The exercise also served as an opportunity to experiment with lesser-known PostGIS functions. Customer locations are proxied by major cities, and only freely available datasets are used.

The process has four stages: importing the datasets, preparing the layers, building areas around airports, and intersecting those areas with major cities.

Importing the data

Because the sales team departs from international airports in Germany, airport locations and classifications are required for that country. Airports can be pulled out of OpenStreetMap data by filtering osm_points on Tag:aeroway=aerodrome, or any other source with current airport locations and classifications will do. For this demo, a freely available CSV from "OurAirports" was used.

Spatial and attributive data for cities in Germany came from OpenStreetMap: the most recent dataset covering Germany was downloaded from Geofabrik and imported into PostGIS with osm2pgsql. The airport, osm_points and osm_polygons (country border) data structures form the foundation of the analysis.

Preparing the layers

Cities are obtained by filtering osm_points on place, which returns only major cities and towns.

1

2

3

4

5

6

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

For the airports, small airports and heliports are removed, since only international airports within Germany matter for this analysis.

1

2

3

4

5

6

7

8

9

10

11

12

13

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)

The resulting core layers are shown below.

catchment areas postgis

Figure 1 Major cities and international airports


Figure 2 Major cities, towns and international airports


Catchment areas around airports

With the layers in place, catchment areas can be generated around the international airports. Rather than creating Voronoi polygons in a preferred GIS client, the task is handled entirely with PostGIS functions.

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

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

Intersecting catchment areas with cities

The next step ranks the most interesting areas by counting major cities per area as a proxy for potential customer counts. Major cities stand in for more detailed customer locations here and should be enriched with further datasets to approximate customer potential more realistically. The query below outputs major city count by polygon; the figure shows the resulting Voronoi polygons extended by those counts. Both the results and the visualization can act as building blocks for optimizing customer acquisition later on.

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

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 ;

Closing notes

The session began with basic spatial analytics and used the opportunity to experiment with less-known PostGIS functions, whose feature set is available for real-world challenges.

Annex

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

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