Work-in-progress: Add geometry functions in Lua

This is intended to show how geometry processing in Lua could work. Only
some example functions are implemented and not for all geometry types,
error detection is limited etc.

Geometry creation functions:
* as_point() (for nodes)
* as_linestring() (for ways)
* as_polygon() (for ways)
* as_multipolygon() (for relations)

Geometry functions:
* area() (for polygons)
* centroid() (for polygons)
* simplify() (for linestrings)

See the file flex-config/custom-geometries.lua for an example that shows
how to use these.

Author: joto

.github/actions/win-install/action.yml

@@ -4,7 +4,7 @@ runs:
   using: composite
   steps:
     - name: Install packages
-      run: vcpkg install bzip2:x64-windows expat:x64-windows zlib:x64-windows proj4:x64-windows boost-system:x64-windows boost-filesystem:x64-windows boost-property-tree:x64-windows lua:x64-windows libpq:x64-windows
+      run: vcpkg install bzip2:x64-windows expat:x64-windows zlib:x64-windows proj4:x64-windows boost-geometry:x64-windows boost-system:x64-windows boost-filesystem:x64-windows boost-property-tree:x64-windows lua:x64-windows libpq:x64-windows
       shell: bash
     - name: Install psycopg2
       run: python -m pip install psycopg2

flex-config/custom-geometries.lua

@@ -0,0 +1,165 @@
+-- This config example file is released into the Public Domain.
+
+-- This is Lua config showing the handling of custom geometries.
+--
+-- This is "work in progress", changes are possible.
+--
+-- NOTE that all the geometry types have a leading "c" for "custom", this is
+--      different than the geometry types before.
+
+local tables = {}
+
+-- This table will get all nodes and areas with an "addr:street" tag, for
+-- areas we'll calculate the centroid.
+tables.addr = osm2pgsql.define_table({
+    name = 'addr',
+    ids = { type = 'any', id_column = 'osm_id', type_column = 'osm_type' },
+    columns = {
+        { column = 'tags', type = 'jsonb' },
+        { column = 'geom', type = 'cpoint', projection = 4326 },
+    }
+})
+
+tables.ways = osm2pgsql.define_way_table('ways', {
+    { column = 'tags', type = 'jsonb' },
+    { column = 'geom', type = 'clinestring', projection = 4326 },
+    { column = 'poly', type = 'cpolygon', projection = 4326 },
+    { column = 'geom3857', type = 'clinestring', projection = 3857 },
+})
+
+tables.major_roads = osm2pgsql.define_way_table('major_roads', {
+    { column = 'tags', type = 'jsonb' },
+    { column = 'geom', type = 'clinestring' },
+    { column = 'sgeom', type = 'clinestring' }
+})
+
+tables.polygons = osm2pgsql.define_area_table('polygons', {
+    { column = 'tags', type = 'jsonb' },
+    { column = 'geom', type = 'cgeometry' },
+    { column = 'area', type = 'real' }
+})
+
+-- Helper function to remove some of the tags we usually are not interested in.
+-- Returns true if there are no tags left.
+function clean_tags(tags)
+    tags.odbl = nil
+    tags.created_by = nil
+    tags.source = nil
+    tags['source:ref'] = nil
+
+    return next(tags) == nil
+end
+
+-- Helper function that looks at the tags and decides if this is possibly
+-- an area.
+function has_area_tags(tags)
+    if tags.area == 'yes' then
+        return true
+    end
+    if tags.area == 'no' then
+        return false
+    end
+
+    return tags.aeroway
+        or tags.amenity
+        or tags.building
+        or tags.harbour
+        or tags.historic
+        or tags.landuse
+        or tags.leisure
+        or tags.man_made
+        or tags.military
+        or tags.natural
+        or tags.office
+        or tags.place
+        or tags.power
+        or tags.public_transport
+        or tags.shop
+        or tags.sport
+        or tags.tourism
+        or tags.water
+        or tags.waterway
+        or tags.wetland
+        or tags['abandoned:aeroway']
+        or tags['abandoned:amenity']
+        or tags['abandoned:building']
+        or tags['abandoned:landuse']
+        or tags['abandoned:power']
+        or tags['area:highway']
+end
+
+function osm2pgsql.process_node(object)
+    if clean_tags(object.tags) then
+        return
+    end
+
+    if object.tags['addr:street'] then
+        tables.addr:add_row({
+           tags = object.tags,
+           geom = object:as_point()
+        })
+    end
+end
+
+function osm2pgsql.process_way(object)
+    if clean_tags(object.tags) then
+        return
+    end
+
+    if object.is_closed and object.tags['addr:street'] then
+        tables.addr:add_row({
+            tags = object.tags,
+            geom = object:as_polygon():centroid()
+        })
+    end
+
+    if object.tags.highway == 'motorway' or
+       object.tags.highway == 'trunk' or
+       object.tags.highway == 'primary' then
+        tables.major_roads:add_row({
+            tags = object.tags,
+            geom = object:as_linestring(),
+            sgeom = object:as_linestring():simplify(100),
+        })
+    end
+
+    -- A closed way that also has the right tags for an area is a polygon.
+    if object.is_closed and has_area_tags(object.tags) then
+        local a = object:as_polygon():area()
+        if a < 0.0000001 then
+            tables.polygons:add_row({
+                tags = object.tags,
+                geom = object:as_polygon(),
+                area = object:as_polygon():area()
+            })
+        end
+    else
+        tables.ways:add_row({
+            tags = object.tags,
+            geom = object:as_linestring(),
+            poly = object:as_polygon(),
+            geom3857 = object:as_linestring()
+        })
+    end
+end
+
+function osm2pgsql.process_relation(object)
+    if clean_tags(object.tags) then
+        return
+    end
+
+    local type = object:grab_tag('type')
+
+    if type == 'multipolygon' then
+        local g = object:as_multipolygon()
+        local a = g:area()
+        if a < 0.0000001 then
+            tables.polygons:add_row({
+                tags = object.tags,
+                geom = g,
+                area = g:area()
+            })
+        end
+    end
+end
+

src/flex-table-column.cpp

@@ -23,7 +23,7 @@ struct column_type_lookup
     table_column_type type;
 };
 
-static std::array<column_type_lookup, 25> const column_types = {
+static std::array<column_type_lookup, 32> const column_types = {
     {{"text", table_column_type::text},
      {"boolean", table_column_type::boolean},
      {"bool", table_column_type::boolean},
@@ -39,6 +39,13 @@ static std::array<column_type_lookup, 25> const column_types = {
      {"json", table_column_type::json},
      {"jsonb", table_column_type::jsonb},
      {"direction", table_column_type::direction},
+     {"cgeometry", table_column_type::custom_geometry},
+     {"cpoint", table_column_type::custom_point},
+     {"clinestring", table_column_type::custom_linestring},
+     {"cpolygon", table_column_type::custom_polygon},
+     {"cmultipoint", table_column_type::custom_multipoint},
+     {"cmultilinestring", table_column_type::custom_multilinestring},
+     {"cmultipolygon", table_column_type::custom_multipolygon},
      {"geometry", table_column_type::geometry},
      {"point", table_column_type::point},
      {"linestring", table_column_type::linestring},
@@ -140,6 +147,20 @@ std::string flex_table_column_t::sql_type_name() const
         return "jsonb";
     case table_column_type::direction:
         return "int2";
+    case table_column_type::custom_geometry:
+        return "Geometry(GEOMETRY, {})"_format(m_srid);
+    case table_column_type::custom_point:
+        return "Geometry(POINT, {})"_format(m_srid);
+    case table_column_type::custom_linestring:
+        return "Geometry(LINESTRING, {})"_format(m_srid);
+    case table_column_type::custom_polygon:
+        return "Geometry(POLYGON, {})"_format(m_srid);
+    case table_column_type::custom_multipoint:
+        return "Geometry(MULTIPOINT, {})"_format(m_srid);
+    case table_column_type::custom_multilinestring:
+        return "Geometry(MULTILINESTRING, {})"_format(m_srid);
+    case table_column_type::custom_multipolygon:
+        return "Geometry(MULTIPOLYGON, {})"_format(m_srid);
     case table_column_type::geometry:
         return "Geometry(GEOMETRY, {})"_format(m_srid);
     case table_column_type::point:

src/flex-table-column.hpp

@@ -32,6 +32,14 @@ enum class table_column_type : uint8_t
 
     direction,
 
+    custom_geometry,
+    custom_point,
+    custom_linestring,
+    custom_polygon,
+    custom_multipoint,
+    custom_multilinestring,
+    custom_multipolygon,
+
     geometry,
     point,
     linestring,
@@ -84,6 +92,12 @@ class flex_table_column_t
                (m_type <= table_column_type::multipolygon);
     }
 
+    bool is_custom_geometry_column() const noexcept
+    {
+        return (m_type >= table_column_type::custom_geometry) &&
+               (m_type <= table_column_type::custom_multipolygon);
+    }
+
     /**
      * Do we need an ST_IsValid() check in the database for this geometry
      * column? If the SRID is 4326 the geometry validity is already assured

src/geom-from-osm.cpp

@@ -17,6 +17,13 @@
 
 namespace geom {
 
+void create_point(geometry_t *geom, osmium::Node const &node)
+{
+    auto &point = geom->set<point_t>();
+    point.set_x(node.location().lon());
+    point.set_y(node.location().lat());
+}
+
 geometry_t create_point(osmium::Node const &node)
 {
     return geometry_t{point_t{node.location()}};
@@ -36,18 +43,22 @@ static void fill_point_list(point_list_t *list,
     }
 }
 
-geometry_t create_linestring(osmium::Way const &way)
+void create_linestring(geometry_t *geom, osmium::Way const &way)
 {
-    geometry_t geom{linestring_t{}};
-    auto &line = geom.get<linestring_t>();
+    auto &line = geom->set<linestring_t>();
 
     fill_point_list(&line, way.nodes());
 
     // Return nullgeom_t if the line geometry is invalid
     if (line.size() <= 1U) {
-        geom.reset();
+        geom->reset();
     }
+}
 
+geometry_t create_linestring(osmium::Way const &way)
+{
+    geometry_t geom{};
+    create_linestring(&geom, way);
     return geom;
 }
 

src/geom-from-osm.hpp

@@ -31,6 +31,8 @@ namespace geom {
  */
 geometry_t create_point(osmium::Node const &node);
 
+void create_point(geometry_t *geom, osmium::Node const &node);
+
 /**
  * Create a linestring geometry from a way. Nodes without location are ignored.
  * Consecutive nodes with the same location will only end up once in the
@@ -44,6 +46,8 @@ geometry_t create_point(osmium::Node const &node);
  */
 geometry_t create_linestring(osmium::Way const &way);
 
+void create_linestring(geometry_t *geom, osmium::Way const &way);
+
 /**
  * Create a polygon geometry from a way.
  *

src/geom-functions.cpp

@@ -7,6 +7,7 @@
  * For a full list of authors see the git log.
  */
 
+#include "geom-boost-adaptor.hpp"
 #include "geom-functions.hpp"
 
 #include <algorithm>
@@ -519,4 +520,48 @@ geometry_t line_merge(geometry_t geom)
     return output;
 }
 
+namespace {
+
+class centroid_helper
+{
+public:
+    centroid_helper(point_t *output) : m_output(output) {}
+
+    void operator()(geom::nullgeom_t const & /*geom*/) const {}
+
+    // XXX not implemented yet
+    void operator()(geom::collection_t const & /*geom*/) const {}
+
+    template <typename GEOM>
+    void operator()(GEOM const &input) const
+    {
+        boost::geometry::centroid(input, *m_output);
+    }
+
+private:
+    point_t *m_output;
+};
+
+} // anonymous namespace
+
+geometry_t centroid(geometry_t const &geom)
+{
+    geom::geometry_t center{point_t{}, geom.srid()};
+    geom.visit(centroid_helper{&center.get<point_t>()});
+    return center;
+}
+
+geometry_t simplify(geometry_t const &geom, double tolerance)
+{
+    geom::geometry_t simplified{linestring_t{}, geom.srid()};
+    if (!geom.is_linestring()) {
+        simplified.reset();
+        return simplified;
+    }
+
+    boost::geometry::simplify(geom.get<linestring_t>(), simplified.get<linestring_t>(), tolerance);
+
+    return simplified;
+}
+
 } // namespace geom

src/geom-functions.hpp

@@ -86,6 +86,10 @@ std::vector<geometry_t> split_multi(geometry_t geom, bool split_multi = true);
  */
 geometry_t line_merge(geometry_t geom);
 
+geometry_t centroid(geometry_t const &geom);
+
+geometry_t simplify(geometry_t const &geom, double tolerance);
+
 } // namespace geom
 
 #endif // OSM2PGSQL_GEOM_FUNCTIONS_HPP

src/init.lua

@@ -179,7 +179,7 @@ local inner_metatable = {
 }
 
 object_metatable = {
-    __index =  {
+    __index = {
         grab_tag = function(data, tag)
             if not tag then
                 error("Missing tag key", 2)
@@ -193,3 +193,7 @@ object_metatable = {
 
 setmetatable(object_metatable.__index, inner_metatable)
 
+-- This will be the metatable for the geometry objects
+geom_metatable = {
+    __index = {}
+}