Extracting Points from Polygons using osm2pgsql and PostGIS

last updated 2022-09-24

A common use case for OpenStreetMap is extracting POIs (points of interest). This is easy for points that are modelled as nodes. However, some features of interest actually appear in the shape of polygons (and sometimes even lines). A typical example is a building (polygon) that is tagged as supermarket, but in our case we are only interested in a single point and not in the shape of the geometry.

A simple approach is to compute a centroid of a polygon. This works well for most common cases, but fails for polygons whose centroid lies outside of its geometry, like for “U”-shaped polygons (see image below). A more advanced approach is to use a function that searches a point on the polygon’s surface. An even more advanced approach is to use the center point of a maximum inscribed circle. All of these options can be triggered with osm2pgsql - a tool for importing OSM data into a PostGIS database.

Adding a centroid can be done directly in the Lua script that defines the OSM data import. It has the function centroid() that computes the centroid of a geometry as seen here:

function osm2pgsql.process_way(object)
    if object.is_closed and has_area_tags(object.tags) and (object.tags.amenity or object.tags.shop) then
        print(inspect(object))

        local row = {
            name = object.tags.name,
            tags = object.tags,
            centroid = object:as_multipolygon():centroid(),
            geom = object:as_multipolygon()
        }

        pois:insert(row)
    else
        return
    end
end

The other two options are done by PostGIS during the import, but are also defined in the Lua script. They make use of generated columns and are configured in the table definition. See the respective section in the osm2pgsql docs:

local pois = osm2pgsql.define_table({
    name = 'pois',
    ids = { type = 'any', type_column = 'osm_type', id_column = 'osm_id' },
    columns = {
        { column = 'name' },
        { column = 'tags', type = 'jsonb' },
        { column = 'centroid', type = 'point', projection = srid },
        { column = 'geom', type = 'geometry', projection = srid, not_null = true },
        {
            column = 'point_on_surface',
            create_only = true,
            sql_type = 'geometry(point, 4326) GENERATED ALWAYS AS (ST_PointOnSurface(geom)) STORED'
        },
        {
            column = 'point_inscribed_circle_center',
            create_only = true,
            sql_type = 'geometry(point, 4326) GENERATED ALWAYS AS ((ST_MaximumInscribedCircle(geom)).center) STORED'
        }
    }
})

These example Lua scripts use these approach(es). Both are released under Public Domain;

• Multiple options
• Point on surface (only)

They can be run like:

osm2pgsql \
   -d {{database-name}}\
   --output=flex \
   --style={{path/to/style.lua}} \
   {{path/to/osm-file.osm.pbf}}

This adds the original geometry and the computed points to the resulting table. PostGIS can handle this, but many file formats like GeoPackage, GeoJSON or Shapefile can not have multiple geometry columns. Luckily GDAL/OGR the tool for converting geodata between different formats can select a geometry column for further processing. All required columns including the geometry column have to be specified using the -sql flag:

ogr2ogr \
  -sql "SELECT {{column_1}}, {{column_2}}, {{column_3}}, {{geom_name}} FROM {{table_name}}" \
  {{path/to/output.gpkg}} \
  PG:dbname={{database-name}}

Optionally the geometry column can be renamed during the export with {{geom_name}} AS {{new_name}}.