Geographic Information Systems Asked by Redoute on January 1, 2021
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:Daten2osm2pgsql-binosm2pgsql.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:Daten2osm2pgsql-binfds.style
C:Daten2osm-2019-11-06.pbf
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;
Answered by JGH on January 1, 2021
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;
Answered by Redoute on January 1, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP