How to get land_area polygons with osm2pgsql?

Question

I want to get the polygons of several administrative land areas. For example https://wiki.openstreetmap.org/wiki/DE:Grenze#Innerstaatliche_Grenzen contains the relation numbers for the land areas of Germany and of five german states that differ from boundary=administrative and are tagged as land_area=administrative. I get entries in osm_polygon for that six boundary relations 51477, 62718, 62782, 28322, 62771 and 51529, but entries for the land_area relations 62781, 2833343, 451087, 62774, 454192 and 62775 are missing.

For osm2pgsql I copied empty.style as fds.style and added a line

way         land_area               text    polygon,nocolumn

(but I notice empty.style has no entry for "boundary" keys).

Data is for Germany and some surrounding areas, downloaded from geofabrik.de and merged with osmosis.

This is the command line I tried for osm2pgsql:

C:\Daten2\osm2pgsql-bin\osm2pgsql.exe --host ubuntu.forplan.local
--username postgres --password --slim --cache 30000 --proj 32632
--hstore --hstore-add-index --keep-coastlines --multi-geometry --prefix osm 
--style C:\Daten2\osm2pgsql-bin\fds.style
C:\Daten2\osm-2019-11-06.pbf

Answer 1

You can add the instruction for boundaries

node,way   boundary     text         linear

But since you have used the --hstore option, you could query it:

SELECT tags from public.planet_osm_polygon 
WHERE tags @> '"boundary"=>"administrative"'::hstore;

Answer 2

I asked osm2pgsql maintainers for support. Maybe this tagging should be changed. Maybe this relations shouldn't be used at all, as land areas can be constructed from boundary=administrative and coastlines (and I desparately miss a land area relation for the Netherlands).

For the meantime I use the attached query to get the polygons. Note that this has to be adjusted: It uses table name prefix "osm" instead of "planet_osm", it contains patches for some osm data errors at the time of my export, it constructs geometries in EPSG:32632 and it doesn't respect member=subarea. Additionally your planet_osm_polygon should have a hstore index; query takes under 20 seconds in my case.

WITH s1 AS (SELECT
    id rid,
    members,
    hstore(tags) tags
    FROM osm_rels
    WHERE hstore(tags) -> 'type' = 'land_area'),
s2 AS (SELECT
    s1.rid,
    m.mid,
    m.member m1,
    lead(m.member) over (partition by s1.rid order by m.mid) m2
    FROM s1, unnest(s1.members) with ordinality m(member, mid)),
s3 AS (SELECT
    rid,
    row_number() over (partition by rid order by mid) mid,
    cast(split_part(m1, 'w', 2) as bigint) wid
    FROM s2
    WHERE m2 in ('inner', 'outer') or
          m1 = 'w639948393'), -- patch data error
s4 AS (SELECT
    s3.rid,
    s3.mid,
    s3.wid,
    w.nodes
    FROM s3 LEFT JOIN osm_ways w ON w.id = s3.wid
UNION SELECT -- patch data error
    62781,
    10000 + row_number() over (order by id),
    id,
    nodes
    FROM osm_ways
    WHERE id in (742213135, 742213136)),
s5 AS (SELECT
    s4.rid, s4.mid,
    wn.wnid, wn.nid
    FROM s4, unnest(s4.nodes) with ordinality wn(nid, wnid)),
s6 AS (SELECT
    s5.rid, s5.mid, s5.wnid, s5.nid,
    st_transform(st_setsrid(st_point(
        n.lon / 10000000.0, n.lat / 10000000.0), 4326), 32632) geom
    FROM s5 LEFT JOIN osm_nodes n ON n.id = s5.nid),
s7 AS (SELECT
    rid,
    st_makeline(geom ORDER BY wnid) geom
    FROM s6
    GROUP BY rid, mid
UNION SELECT -- patch data error
    5409700,
    st_makeline(
    (SELECT DISTINCT geom FROM s6 WHERE nid = 3605382213),
    (SELECT DISTINCT geom FROM s6 WHERE nid = 4760866421)
    )),
s8 AS (SELECT
    rid,
    st_multi(st_buildarea(st_collect(geom))) geom
    FROM s7
    GROUP BY rid)
SELECT
    s1.rid,
    s1.tags,
    s8.geom
    FROM s1 LEFT JOIN s8 ON s8.rid = s1.rid
    ORDER BY s1.rid;